From 868cb606ce8384312f0e5723243de18e1e7dbdfa Mon Sep 17 00:00:00 2001 From: Ronan Fleming Date: Mon, 14 Oct 2024 19:39:31 +0100 Subject: [PATCH 1/3] refactored parameters gurobi and mosek --- initCobraToolbox.m | 53 +-- src/analysis/FBA/optimizeCbModel.m | 126 +++-- .../wholeBody/PSCMToolbox/optimizeWBModel.m | 2 +- src/base/install/updateCobraToolbox.m | 2 +- src/base/solvers/buildOptProblemFromModel.m | 187 ++++---- .../cplex/buildCplexProblemFromCOBRAStruct.m | 2 +- .../cplex/setCplexParametersForProblem.m | 52 +- .../entropicFBA/entropicFluxBalanceAnalysis.m | 58 ++- .../solvers/entropicFBA/mosekParamSetEFBA.m | 125 ----- .../entropicFBA/processConcConstraints.m | 11 +- .../entropicFBA/processFluxConstraints.m | 13 +- src/base/solvers/entropicFBA/solveCobraEP.m | 446 +++++++----------- .../solvers/getSetSolver/changeCobraSolver.m | 40 +- src/base/solvers/gurobi/setGurobiParam.m | 342 ++++++++++++++ src/base/solvers/init/configEnvVars.m | 2 +- .../{entropicFBA => mosek}/mosekParamStrip.m | 0 .../solvers/{msk => mosek}/parseMskResult.m | 15 +- src/base/solvers/mosek/setMosekParam.m | 373 +++++++++++++++ src/base/solvers/param/getCobraSolverParams.m | 101 ++-- .../getCobraSolverParamsOptionsForType.m | 33 +- src/base/solvers/param/mergeCobraParams.m | 11 + .../solvers/param/parseSolverParameters.m | 63 ++- src/base/solvers/solveCobraLP.m | 199 ++------ src/base/solvers/solveCobraQP.m | 143 ++---- .../createToyModelForLifting.m | 0 .../testLiftModel.m | 0 .../base/testSolvers/runLPvariousSolvers.m | 4 +- .../base/testSolvers/testSolveCobraLP.m | 4 +- 28 files changed, 1445 insertions(+), 962 deletions(-) delete mode 100644 src/base/solvers/entropicFBA/mosekParamSetEFBA.m create mode 100644 src/base/solvers/gurobi/setGurobiParam.m rename src/base/solvers/{entropicFBA => mosek}/mosekParamStrip.m (100%) rename src/base/solvers/{msk => mosek}/parseMskResult.m (97%) create mode 100644 src/base/solvers/mosek/setMosekParam.m create mode 100644 src/base/solvers/param/mergeCobraParams.m rename test/verifiedTests/base/{testSolvers => testLifting}/createToyModelForLifting.m (100%) rename test/verifiedTests/base/{testSolvers => testLifting}/testLiftModel.m (100%) diff --git a/initCobraToolbox.m b/initCobraToolbox.m index 3404c44ef5..ca323154b3 100644 --- a/initCobraToolbox.m +++ b/initCobraToolbox.m @@ -16,11 +16,8 @@ function initCobraToolbox(updateToolbox) % % initCobraToolbox % -or- -% changeCobraSolver('gurobi'); -% changeCobraSolver('gurobi', 'MILP'); -% changeCobraSolver('tomlab_cplex', 'QP'); -% changeCobraSolver('tomlab_cplex', 'MIQP'); -% changeCbMapOutput('svg'); +% changeCobraSolver('gurobi','all'); +% changeCobraSolver('mosek', 'CLP'); % % Maintained by Ronan M.T. Fleming @@ -362,18 +359,7 @@ function initCobraToolbox(updateToolbox) % requires the solver compatibility to be re-read at each initialisation clear isCompatible -%These default tolerances are based on the default values for the Gurobi LP -%solver. Do not change them without first consulting with other developers. -%https://www.gurobi.com/documentation/9.0/refman/parameters.html -% (primal) feasibility tolerance -changeCobraSolverParams('LP', 'feasTol', 1e-6); -% (dual) optimality tolerance -changeCobraSolverParams('LP', 'optTol', 1e-6); -% (primal) feasibility tolerance -changeCobraSolverParams('EP', 'feasTol', 1e-8); -% (dual) optimality tolerance -changeCobraSolverParams('EP', 'optTol', 1e-12); % Check that SBML toolbox is installed and accessible if ~exist('TranslateSBML', 'file') @@ -409,10 +395,8 @@ function initCobraToolbox(updateToolbox) %Define a set of "use first" solvers, other supported solvers will also be added to the struct. %This allows to assign them in any order but keep the most commonly used ones on top of the struct. SOLVERS = struct('gurobi',struct(),... - 'ibm_cplex',struct(),... - 'tomlab_cplex',struct(),... - 'glpk',struct(),... 'mosek',struct(),... + 'glpk',struct(),... 'matlab',struct()); % active support - supported solvers @@ -423,17 +407,17 @@ function initCobraToolbox(updateToolbox) SOLVERS.quadMinos.type = {'LP'}; SOLVERS.dqqMinos.type = {'LP','QP'}; SOLVERS.matlab.type = {'LP', 'NLP'}; -% active support of cplex interfaces - supported solvers -SOLVERS.cplex_direct.type = {'LP', 'MILP', 'QP'}; -SOLVERS.ibm_cplex.type = {'LP', 'MILP', 'QP', 'MIQP'}; -SOLVERS.cplexlp.type = {'LP'}; -SOLVERS.tomlab_cplex.type = {'LP', 'MILP', 'QP', 'MIQP'}; + % passive support - solver interfaces SOLVERS.qpng.type = {'QP'}; SOLVERS.tomlab_snopt.type = {'NLP'}; % legacy solvers +SOLVERS.cplex_direct.type = {'LP', 'MILP', 'QP'}; +SOLVERS.ibm_cplex.type = {'LP', 'MILP', 'QP', 'MIQP'}; +SOLVERS.cplexlp.type = {'LP'}; +SOLVERS.tomlab_cplex.type = {'LP', 'MILP', 'QP', 'MIQP'}; %SOLVERS.gurobi_mex.type = {'LP', 'MILP', 'QP', 'MIQP'}; %SOLVERS.lindo_old.type = {'LP'}; %SOLVERS.lindo_legacy.type = {'LP'}; @@ -441,7 +425,6 @@ function initCobraToolbox(updateToolbox) %SOLVERS.opti.type = {'LP', 'MILP', 'QP', 'MIQP', 'NLP'}; % definition of category of solvers with active support - SOLVERS.dqqMinos.categ = 'active'; SOLVERS.glpk.categ = 'active'; SOLVERS.gurobi.categ = 'active'; @@ -450,11 +433,6 @@ function initCobraToolbox(updateToolbox) SOLVERS.pdco.categ = 'active'; SOLVERS.quadMinos.categ = 'active'; -SOLVERS.cplex_direct.categ = 'active'; -SOLVERS.ibm_cplex.categ = 'active'; -SOLVERS.cplexlp.categ = 'active'; -SOLVERS.tomlab_cplex.categ = 'active'; - % definition of category of solvers with passive support SOLVERS.qpng.categ = 'passive'; SOLVERS.tomlab_snopt.categ = 'passive'; @@ -464,6 +442,10 @@ function initCobraToolbox(updateToolbox) %SOLVERS.lindo_old.categ = 'legacy'; %SOLVERS.lindo_legacy.categ = 'legacy'; SOLVERS.lp_solve.categ = 'legacy'; +SOLVERS.cplex_direct.categ = 'legacy'; +SOLVERS.ibm_cplex.categ = 'legacy'; +SOLVERS.cplexlp.categ = 'legacy'; +SOLVERS.tomlab_cplex.categ = 'legacy'; %SOLVERS.opti.categ = 'legacy'; % definition of categories of solvers @@ -484,17 +466,19 @@ function initCobraToolbox(updateToolbox) end end -% check the installation of the solver +supportedSolversNames = setdiff(supportedSolversNames,{'cplex_direct','ibm_cplex','cplexlp','tomlab_cplex'}); + +% check the installation of the solver - except cplex for i = 1:length(supportedSolversNames) if 0 %set to 1 to debug a new solver disp(supportedSolversNames{i}) - if strcmp(supportedSolversNames{i},'quadMinos') - pause(0.1) + if strcmp(supportedSolversNames{i},'quadMinos') + disp(supportedSolversNames{i}) end end %We will validate all solvers in init. After this, all solvers are %checked, whether they actually work and the SOLVERS field is set. - [solverOK,solverInstalled] = changeCobraSolver(supportedSolversNames{i},SOLVERS.(supportedSolversNames{i}).type{1},0, 2); + [solverOK,solverInstalled] = changeCobraSolver(supportedSolversNames{i},SOLVERS.(supportedSolversNames{i}).type{1},0, 1); if strcmp(supportedSolversNames{i},'gurobi') && 0%use fordebugging disp(supportedSolversNames{i}); end @@ -667,7 +651,6 @@ function initCobraToolbox(updateToolbox) % use Gurobi (if installed) as the default solver for LP, QP and MILP problems changeCobraSolver('gurobi', 'ALL', 0); -%changeCobraSolver('ibm_cplex', 'QP', 0); %until problem with gurobi QP sorted % check if a new update exists if installedGit && ENV_VARS.printLevel && status_curl == 0 && contains(result_curl, ' 200') && updateToolbox diff --git a/src/analysis/FBA/optimizeCbModel.m b/src/analysis/FBA/optimizeCbModel.m index 731bcdb1a9..2e5be448d2 100644 --- a/src/analysis/FBA/optimizeCbModel.m +++ b/src/analysis/FBA/optimizeCbModel.m @@ -61,6 +61,13 @@ % * ctrs `k x 1` Cell Array of Strings giving IDs of the coupling constraints % * dsense - `k x 1` character array with entries in {L,E,G} % +% * `.evars` : evars x 1 Column Cell Array of Strings IDs of the additional variables +% * `.E` : n x evars The additional Variable Matrix +% * `.evarub`: evars x 1 The upper bounds of the variables from E; +% * `.evarlb`: evars x 1 The lower bounds of the variables from E; +% * `.evarc` : evars x 1 The objective coefficients of the variables from E; +% * `.D` : k x evars The matrix coupling additional Constraints (form C), with additional Variables (from E); +% % * g0 - `n x 1` weights on zero norm, where positive is minimisation, negative is maximisation, zero is neither. % * g1 - `n x 1` weights on one norm, where positive is minimisation, negative is maximisation, zero is neither. % * g2 - `n x 1` weights on two norm @@ -152,11 +159,11 @@ % % OUTPUT: % solution: solution object: -% +% * obj - Primal objective value (sum of f0,f1,f2 terms, ignoring NaN) % * f - Linear objective value (from LP problem) % * f0 - Zero-norm objective value -% * f1 - One-norm objective value -% * f2 - Two-norm objective value +% * f1 - Linear part of objective value (c'*v or ||v||_1) +% * f2 - Quadratic part of objective value % * v - Reaction rates (Optimal primal variable, legacy FBAsolution.x) % * y - Dual to the matrix inequality constraints (Shadow prices) % * w - Dual to the box constraints (Reduced costs) @@ -222,18 +229,41 @@ % Figure out linear objective sense -if exist('osenseStr', 'var') +if exist('osenseStr', 'var') if isempty(osenseStr) model.osenseStr = 'max'; else - model.osenseStr = osenseStr; - end -else - if isfield(model, 'osenseStr') - model.osenseStr = model.osenseStr; - else - model.osenseStr = 'max'; + %second argument may be the parameter structure + if isstruct(osenseStr) + if exist('param','var') + error('osenseStr is a structure and param structure is also present') + else + %second argument is be the parameter structure + param = osenseStr; + if isfield(param,'osenseStr') && ~isfield(model, 'osenseStr') + model.osenseStr = param.osenseStr; + else + model.osenseStr = 'max'; + end + if isfield(param,'minNorm') + minNorm=param.minNorm; + else + minNorm=[]; + end + if isfield(param,'allowLoops') + allowLoops=param.allowLoops; + else + allowLoops=1; + end + end + end end + % % override if osenseStr already in the model + % if isfield(model, 'osenseStr') + % model.osenseStr = model.osenseStr; + % else + % model.osenseStr = 'max'; + % end end if ~exist('param','var') @@ -246,7 +276,7 @@ if minNorm == true minNorm = 1e-6; else - minNorm = 0; + minNorm = []; end end @@ -367,31 +397,31 @@ end %by default, do linear optimisation unless not required -doLinearOptimisation = 1; +doLinearOptimisationFirst = 1; %if there is no linear objective, do a QP if all(model.c==0) - doLinearOptimisation = 0; + doLinearOptimisationFirst = 0; end % if there is no linear objective and no quadratic objective, do an LP if isempty(minNorm) - doLinearOptimisation = 1; + doLinearOptimisationFirst = 1; end % If there is are linear and quadrative objectives but the bounds on the % corresponding reaction are fixed, then there is no need to solve an LP % first, so do a QP if all( (model.lb == model.ub & model.c~=0) == (model.c~=0)) && ~isempty(minNorm) - doLinearOptimisation = 0; + doLinearOptimisationFirst = 0; end % If this is a quadratically regularised LP, go straight to QP % TODO This is a hack of the param.minNorm to direct solution to QRLP or QRQP if isfield(param,'solveWBMmethod') if any(strcmp(param.solveWBMmethod,{'QRLP','QRQP'})) - doLinearOptimisation = 0; - param.minNormWBM = minNorm; - minNorm = param.solveWBMmethod; + doLinearOptimisationFirst = 0; + param.minNormWBM = param.minNorm; + param.minNorm = param.solveWBMmethod; else param.solveWBMmethod = []; end @@ -401,13 +431,15 @@ % build the optimization problem -optProblem = buildOptProblemFromModel(model,verify); +optProblem = buildOptProblemFromModel(model,verify,param); % save the original size of the problem [~,nTotalVars] = size(optProblem.A); % nTotalVars needed even if optProblem not used for an LP + + %% t1 = clock; -if doLinearOptimisation +if doLinearOptimisationFirst if allowLoops && ~strcmp(minNorm,'optimizeCardinality') clear model end @@ -444,7 +476,7 @@ %only run if minNorm is not empty, and either there is no linear objective %or there is a linear objective and the LP problem solved to optimality -if (doLinearOptimisation==0 && ~isempty(minNorm)) || (doLinearOptimisation==1 && solution.stat==1 && ~isempty(minNorm)) +if (doLinearOptimisationFirst==0 && ~isempty(minNorm)) || (doLinearOptimisationFirst==1 && solution.stat==1 && ~isempty(minNorm)) if strcmp(minNorm, 'optimizeCardinality') % DC programming for solving the cardinality optimization problem @@ -519,7 +551,7 @@ % * .thetaMultiplier - at each iteration: theta = theta*thetaMultiplier % * .eta - Smallest value considered non-zero (Default value feasTol*1000) - if doLinearOptimisation + if doLinearOptimisationFirst optProblem2 = optProblem; optProblem2.A = [optProblem.A ; optProblem.c']; optProblem2.b = [optProblem.b ; objectiveLP]; @@ -555,7 +587,7 @@ % lb <= v <= ub % Define the constraints structure - if doLinearOptimisation + if doLinearOptimisationFirst optProblem2.A = [optProblem.A ; optProblem.c']; optProblem2.b = [optProblem.b ; objectiveLP]; optProblem2.csense = [optProblem.csense;'E']; @@ -676,8 +708,8 @@ solution = solveCobraLP(optProblem2, param); elseif strcmp(minNorm, 'QRLP') - - buildOptProblemFromModel_param.method='QRLP'; + buildOptProblemFromModel_param = param; + buildOptProblemFromModel_param.minNorm = param.minNormWBM; optProblem = buildOptProblemFromModel(model, 0, buildOptProblemFromModel_param); solutionQRLP = solveCobraQP(optProblem,param); @@ -698,8 +730,8 @@ solution.q = -solution.q; % lb + p <= x <= ub + q elseif strcmp(minNorm, 'QRQP') - - buildOptProblemFromModel_param.method='QRQP'; + buildOptProblemFromModel_param = param; + buildOptProblemFromModel_param.minNorm = param.minNormWBM; optProblem = buildOptProblemFromModel(model, 0, buildOptProblemFromModel_param); solution = solveCobraQP(optProblem); @@ -737,7 +769,7 @@ end % quadratic minimization of the norm. - if doLinearOptimisation + if doLinearOptimisationFirst % set previous linear optimum as constraint. optProblem2 = optProblem; optProblem2.A = [optProblem.A;optProblem.c']; @@ -757,14 +789,13 @@ solution = solveCobraMIQP(MIQPproblem); end else - optProblem.F = spdiags(minNorm,0,nTotalVars,nTotalVars); if allowLoops %quadratic optimization will get rid of the loops unless you are maximizing a flux which is %part of a loop. By definition, exchange reactions are not part of these loops, more %properly called stoichiometrically balanced cycles. - solution = solveCobraQP(optProblem); + solution = solveCobraQP(optProblem,param); else %this is slow, but more useful than minimizing the Euclidean norm if one is trying to %maximize the flux through a reaction in a loop. e.g. in flux variablity analysis @@ -835,23 +866,35 @@ %the value of the second part of the objective depends on the norm switch minNorm case 'empty' + solution.f0 = NaN; solution.f1 = optProblem.c'*solution.full(1:nTotalVars,1); + solution.f2 = NaN; case 'zero' %zero norm solution.f0 = sum(abs(solution.full(1:nTotalVars,1)) > feasTol); + solution.f1 = NaN; + solution.f2 = NaN; case 'one' + solution.f0 = NaN; %one norm solution.f1 = sum(abs(solution.full(1:nTotalVars,1))); + solution.f2 = NaN; case 'two' + solution.f0 = NaN; if isfield(optProblem,'c') solution.f1 = optProblem.c'*solution.full(1:nTotalVars,1); if isfield(solution,'objLinear') solution = rmfield(solution,'objLinear'); end else - solution.f1 = 0; + solution.f1 = NaN; + end + if isfield(optProblem,'F') + solution.f2 = 0.5*solution.full'*optProblem.F*solution.full; + else + disp(param) + warning('optProblem.F missing') end - solution.f2 = 0.5*solution.full'*optProblem.F*solution.full; if isfield(solution,'objQuadratic') solution.f2 = solution.objQuadratic; solution = rmfield(solution,'objQuadratic'); @@ -859,12 +902,17 @@ otherwise if exist('LPproblem2','var') if isfield(optProblem2,'F') - solution.f0 = 0; + solution.f0 = NaN; solution.f1 = optProblem.c'*solution.full(1:nTotalVars,1); solution.f2 = 0.5*solution.full'*optProblem2.F*solution.full; end + else + solution.f0 = NaN; + solution.f1 = NaN; + solution.f2 = NaN; end end + solution.obj = sum([solution.f0,solution.f1,solution.f2],'omitnan'); %primal optimal variables solution.v = solution.full(1:nRxns); @@ -922,7 +970,7 @@ solution.time = etime(clock, t1); - fieldOrder = {'f';'f0';'f1';'f2';'v';'y';'w';'s';'solver';'method';'stat';'origStat';'time';'basis';'vars_v';'vars_w';'ctrs_y';'ctrs_s';'x';'full';'obj';'rcost';'dual';'slack'}; + fieldOrder = {'f';'f0';'f1';'f2';'v';'y';'w';'s';'solver';'lpmethod';'qpmethod';'stat';'origStat';'time';'basis';'vars_v';'vars_w';'ctrs_y';'ctrs_s';'x';'full';'obj';'rcost';'dual';'slack'}; % reorder fields for better readability currentfields = fieldnames(solution); presentfields = ismember(fieldOrder,currentfields); @@ -932,6 +980,7 @@ if 0 %return NaN of correct dimensions if problem does not solve properly solution.f = NaN; + solution.v = NaN*ones(nRxns,1); solution.y = NaN*ones(nMets,1); solution.w = NaN*ones(nRxns,1); @@ -947,7 +996,11 @@ else %return empty fields if problem does not solve properly (backward %compatible) + solution.obj = NaN; solution.f = NaN; + solution.f0 = NaN; + solution.f1 = NaN; + solution.f2 = NaN; solution.v = []; solution.y = []; solution.w = []; @@ -968,9 +1021,6 @@ if 1 %this may not be very backward compatible %remove fields coming from solveCobraLP/QP but not part of the specification %of the output from optimizeCbModel - if isfield(solution,'obj') - solution = rmfield(solution,'obj'); - end if isfield(solution,'full') solution = rmfield(solution,'full'); end diff --git a/src/analysis/wholeBody/PSCMToolbox/optimizeWBModel.m b/src/analysis/wholeBody/PSCMToolbox/optimizeWBModel.m index e9e788ae7d..107f4445c7 100644 --- a/src/analysis/wholeBody/PSCMToolbox/optimizeWBModel.m +++ b/src/analysis/wholeBody/PSCMToolbox/optimizeWBModel.m @@ -177,7 +177,7 @@ % is especially well suited for models with a wide range of coefficients in the constraint matrix rows or columns. % Settings 1 and 3 are not as directly connected to any specific model characteristics, so experimentation with both % settings may be needed to assess performance impact. - param.scaleFlag=0; + param.ScaleFlag=0; case 'ibm_cplex' % https://www.ibm.com/docs/en/icos/12.10.0?topic=infeasibility-coping-ill-conditioned-problem-handling-unscaled-infeasibilities diff --git a/src/base/install/updateCobraToolbox.m b/src/base/install/updateCobraToolbox.m index 5af2667281..cce06c1cb4 100644 --- a/src/base/install/updateCobraToolbox.m +++ b/src/base/install/updateCobraToolbox.m @@ -179,7 +179,7 @@ function updateCobraToolbox(fetchAndCheckOnly) %disp([codeBaseDir filesep forkNames{j}]) %disp([toolboxDir filesep dirNames{j}]) fprintf('%s\n',['> Adding path to ' codeBaseDir filesep forkNames{j} ', and removing path to ' toolboxDir filesep dirNames{j}]); - rmpath([toolboxDir filesep dirNames{j}]) + rmpath(genpath([toolboxDir filesep dirNames{j}])) addpath(genpath([codeBaseDir filesep forkNames{j}])) end end diff --git a/src/base/solvers/buildOptProblemFromModel.m b/src/base/solvers/buildOptProblemFromModel.m index 9b88a68365..5d47cece35 100644 --- a/src/base/solvers/buildOptProblemFromModel.m +++ b/src/base/solvers/buildOptProblemFromModel.m @@ -155,13 +155,6 @@ model.c = zeros(nRxn,1); end -if ~isfield(param, 'method') - if isfield(model,'F') - param.method = 'QP'; - else - param.method = 'LP'; - end -end % case 'LP' if ~modelC && ~modelE @@ -195,86 +188,122 @@ end end -switch param.method - case 'LP' - %nothing to do - done above already - case 'QP' - if modelE - optProblem.F = spdiags(zeros(size(optProblem.A,2),1),0,size(optProblem.A,2),size(optProblem.A,2)); - %assume that the remainder of the variables are not being quadratically - %minimised - optProblem.F(1:size(model.F,1),1:size(model.F,1)) = model.F; - else - optProblem.F = model.F; - end - case 'QRLP' - [m,n]=size(optProblem.A); - optProblem.A = [... - % v r p z - optProblem.A, speye(m,m), sparse(m,n), sparse(m,n); - speye(n,n), sparse(n,m), -speye(n,n), -speye(n,n)]; % x - p = z; +if isfield(model,'F') + if modelE + optProblem.F = spdiags(zeros(size(optProblem.A,2),1),0,size(optProblem.A,2),size(optProblem.A,2)); + %assume that the remainder of the variables are not being quadratically + %minimised + optProblem.F(1:size(model.F,1),1:size(model.F,1)) = model.F; + else + optProblem.F = model.F; + end +end - optProblem.b = [optProblem.b;sparse(n,1)]; - optProblem.lb = [-inf(2*n+m,1); model.lb]; - optProblem.ub = [ inf(2*n+m,1); model.ub]; - optProblem.csense = [optProblem.csense;repmat('E',n,1)]; - optProblem.c = [model.c; sparse(m+2*n,1)]; - optProblem.F = speye(size(optProblem.A,2),size(optProblem.A,2))*1e-16;%small amound of regularistion makes matrix positive definite numerically - optProblem.F(n+1:2*n+m,n+1:2*n+m) = speye(n+m,n+m)*max(abs(model.c))*100; %regularisation must dominate linear objective - %dimensions needed to extract non-regularised part of solution - optProblem.m = m; - optProblem.n = n; +if isfield(param,'solveWBMmethod') && ~isempty(param.solveWBMmethod) + switch param.solveWBMmethod + case 'LP' + %nothing to do - done above already + case 'QP' + %nothing to do - done above already + case 'QRLP' + [m,n]=size(optProblem.A); + optProblem.A = [... + % v r p z + optProblem.A, speye(m,m), sparse(m,n), sparse(m,n); + speye(n,n), sparse(n,m), -speye(n,n), -speye(n,n)]; % x - p = z; - if 0 %use to test if regularisation enough - F2 = optProblem.F; - %This line modifies the diagonal elements of F2 by setting them to zero. - F2(1:size(optProblem.F,1):end)=0; - if all(all(F2)) == 0 - %only nonzeros in QPproblem.F are on the diagonal - try - %try cholesky decomposition - B = chol(optProblem.F); - catch - optProblem.F = optProblem.F + diag((diag(optProblem.F)==0)*1e-16); - end - try - B = chol(optProblem.F); - catch - error('QPproblem.F only has non-zeros along the main diagnoal and is still not positive semidefinite after adding 1e-16') + optProblem.b = [optProblem.b;sparse(n,1)]; + optProblem.lb = [-inf(2*n+m,1); model.lb]; + optProblem.ub = [ inf(2*n+m,1); model.ub]; + optProblem.csense = [optProblem.csense;repmat('E',n,1)]; + optProblem.c = [model.c; sparse(m+2*n,1)]; + optProblem.F = speye(size(optProblem.A,2),size(optProblem.A,2))*1e-16;%small amound of regularistion makes matrix positive definite numerically + optProblem.F(n+1:2*n+m,n+1:2*n+m) = speye(n+m,n+m)*max(abs(model.c))*100; %regularisation must dominate linear objective + %dimensions needed to extract non-regularised part of solution + optProblem.m = m; + optProblem.n = n; + + if 0 %use to test if regularisation enough + F2 = optProblem.F; + %This line modifies the diagonal elements of F2 by setting them to zero. + F2(1:size(optProblem.F,1):end)=0; + if all(all(F2)) == 0 + %only nonzeros in QPproblem.F are on the diagonal + try + %try cholesky decomposition + B = chol(optProblem.F); + catch + optProblem.F = optProblem.F + diag((diag(optProblem.F)==0)*1e-16); + end + try + B = chol(optProblem.F); + catch + error('QPproblem.F only has non-zeros along the main diagnoal and is still not positive semidefinite after adding 1e-16') + end end end - end - case 'QRQP' - if modelE - optProblem.F = spdiags(zeros(size(optProblem.A,2),1),0,size(optProblem.A,2),size(optProblem.A,2)); - %assume that the remainder of the variables are not being quadratically - %minimised - optProblem.F(1:size(model.F,1),1:size(model.F,1)) = model.F; - else - optProblem.F = model.F; - end + case 'QRQP' + if modelE + optProblem.F = spdiags(zeros(size(optProblem.A,2),1),0,size(optProblem.A,2),size(optProblem.A,2)); + %assume that the remainder of the variables are not being quadratically + %minimised + optProblem.F(1:size(model.F,1),1:size(model.F,1)) = model.F; + else + optProblem.F = model.F; + end - [m,n]=size(optProblem.A); - optProblem.A = [... - optProblem.A, speye(m,m), sparse(m,n), sparse(m,n); % A*x + r <=> b - speye(n,n), sparse(n,m), -speye(n,n), -speye(n,n)]; % x - p = z; + [m,n]=size(optProblem.A); + optProblem.A = [... + optProblem.A, speye(m,m), sparse(m,n), sparse(m,n); % A*x + r <=> b + speye(n,n), sparse(n,m), -speye(n,n), -speye(n,n)]; % x - p = z; - optProblem.b = [optProblem.b;sparse(n,1)]; - optProblem.lb = [-inf(2*n+m,1); model.lb]; % lb <= x - p ----> lb <= z - optProblem.ub = [ inf(2*n+m,1); model.ub]; % x - p <= ub ----> z <= ub - optProblem.csense = model.csense; - optProblem.c = [model.c; sparse(m+2*n,1)]; - optProblem.F = sparse(size(optProblem.A,2),size(optProblem.A,2)); - optProblem.F(1:n,1:n) = model.F; - optProblem.F(n+1:2*n+m,n+1:2*n+m) = speye(n+m,n+m)*max(abs(model.c))*100; %regularisation must dominate linear and quadratic objective - %dimensions needed to extract non-regularised part of solution - optProblem.m = m; - optProblem.n = n; + optProblem.b = [optProblem.b;sparse(n,1)]; + optProblem.lb = [-inf(2*n+m,1); model.lb]; % lb <= x - p ----> lb <= z + optProblem.ub = [ inf(2*n+m,1); model.ub]; % x - p <= ub ----> z <= ub + optProblem.csense = model.csense; + optProblem.c = [model.c; sparse(m+2*n,1)]; + optProblem.F = sparse(size(optProblem.A,2),size(optProblem.A,2)); + optProblem.F(1:n,1:n) = model.F; + optProblem.F(n+1:2*n+m,n+1:2*n+m) = speye(n+m,n+m)*max(abs(model.c))*100; %regularisation must dominate linear and quadratic objective + %dimensions needed to extract non-regularised part of solution + optProblem.m = m; + optProblem.n = n; - otherwise - error('param.method unrecognised') + otherwise + error('param.method unrecognised') + end end - [~,optProblem.osense] = getObjectiveSense(model); + +if isfield(param,'debug') && param.debug + switch param.solver + case 'mosek' + % names + % This structure is used to store all the names of individual items in the optimization problem such as the constraints and the variables. + % + % Fields + % name (string) – contains the problem name. + % + % obj (string) – contains the name of the objective. + % + % con (cell) – a cell array where names.con{i} contains the name of the + % -th constraint. + % + % var (cell) – a cell array where names.var{j} contains the name of the + % -th variable. + optProblem.names.name='optimizeCbModel'; + optProblem.names.obj=model.rxns{model.c~=0}; + if isfield(model,'ctrs') + optProblem.names.con=[model.mets;model.ctrs]; + else + optProblem.names.con=model.mets; + end + if isfield(model,'evars') + optProblem.names.con=[model.mets;model.ctrs]; + else + optProblem.names.con=model.mets; + end + end +end \ No newline at end of file diff --git a/src/base/solvers/cplex/buildCplexProblemFromCOBRAStruct.m b/src/base/solvers/cplex/buildCplexProblemFromCOBRAStruct.m index 8891b9bac5..61f1bd766e 100644 --- a/src/base/solvers/cplex/buildCplexProblemFromCOBRAStruct.m +++ b/src/base/solvers/cplex/buildCplexProblemFromCOBRAStruct.m @@ -75,7 +75,7 @@ %always set the problem to minimise, but change the linear objective sign cplexProblem.Model.sense = 'minimize'; if isfield(Problem,'c') - if isfield(Problem,'osense') && isfield(Problem,'c') + if isfield(Problem,'osense') cplexProblem.Model.obj = Problem.osense*Problem.c; else cplexProblem.Model.obj = Problem.c; diff --git a/src/base/solvers/cplex/setCplexParametersForProblem.m b/src/base/solvers/cplex/setCplexParametersForProblem.m index e4398ab27d..f53b7f0ca4 100644 --- a/src/base/solvers/cplex/setCplexParametersForProblem.m +++ b/src/base/solvers/cplex/setCplexParametersForProblem.m @@ -19,6 +19,12 @@ %set the default parameters so we can see what they are cplexProblem.setDefault; +%TODO ? add these +% valDef.DATACHECK = 1; +% valDef.DEPIND = 1; +% valDef.checkNaN = 0; +% valDef.warning = 0; + % set the printLevel to the cobra Parameters cplexProblem.Param.output.writelevel.Cur = problemTypeParams.printLevel; cplexProblem.Param.barrier.display.Cur = problemTypeParams.printLevel; @@ -67,12 +73,56 @@ end % set tolerances + +% simplex.tolerances.feasibility +% Specifies the feasibility tolerance, that is, the degree to which values of the basic variables +% calculated by the simplex method may violate their bounds. CPLEX® does not use this tolerance to +% relax the variable bounds nor to relax right hand side values. This parameter specifies an +% allowable violation. Feasibility influences the selection of an optimal basis and can be reset +% to a higher value when a problem is having difficulty maintaining feasibility during optimization. +% You can also lower this tolerance after finding an optimal solution if there is any doubt +% that the solution is truly optimal. If the feasibility tolerance is set too low, CPLEX may falsely +% conclude that a problem is infeasible. If you encounter reports of infeasibility during Phase II of +% the optimization, a small adjustment in the feasibility tolerance may improve performance. +% Values +% Any number from 1e-9 to 1e-1; default: 1e-06. cplexProblem.Param.simplex.tolerances.feasibility.Cur = problemTypeParams.feasTol; -cplexProblem.Param.simplex.tolerances.optimality.Cur = problemTypeParams.optTol; + +% network.tolerances.feasibility +% Specifies feasibility tolerance for network primal optimization. The feasibility tolerance specifies +% the degree to which the flow value of a model may violate its bounds. This tolerance influences +% the selection of an optimal basis and can be reset to a higher value when a problem is having +% difficulty maintaining feasibility during optimization. You may also wish to lower this tolerance +% after finding an optimal solution if there is any doubt that the solution is truly optimal. +% If the feasibility tolerance is set too low, CPLEX may falsely conclude that a problem is infeasible. +% If you encounter reports of infeasibility during Phase II of the optimization, a small adjustment +% in the feasibility tolerance may improve performance. +% Values +% Any number from 1e-11 to 1e-1; default: 1e-6. cplexProblem.Param.network.tolerances.feasibility.Cur = problemTypeParams.feasTol; + +% Influences the reduced-cost tolerance for optimality. This parameter governs +% how closely CPLEX must approach the theoretically optimal solution. +% The simplex algorithm halts when it has found a basic feasible solution with +% all reduced costs nonnegative. CPLEX uses this optimality tolerance to make +% the decision of whether or not a given reduced cost should be considered nonnegative. +% CPLEX considers "nonnegative" a negative reduced cost having absolute value less +% than the optimality tolerance. For example, if your optimality tolerance is set +% to 1e-6, then CPLEX considers a reduced cost of -1e-9 as nonnegative for +% the purpose of deciding whether the solution is optimal. +% Values +% Any number from 1e-9 to 1e-1; default: 1e-06. +cplexProblem.Param.simplex.tolerances.optimality.Cur = problemTypeParams.optTol; + +% network.tolerances.optimality +% Specifies the optimality tolerance for network optimization; that is, +% the amount a reduced cost may violate the criterion for an optimal solution. +% Values +% Any number from 1e-11 to 1e-1; default: 1e-6. cplexProblem.Param.network.tolerances.optimality.Cur = problemTypeParams.optTol; + %https://www.ibm.com/support/knowledgecenter/SSSA5P_12.7.0/ilog.odms.cplex.help/CPLEX/Parameters/topics/BarEpComp.html %Sets the tolerance on complementarity for convergence. The barrier algorithm terminates with an optimal solution if the relative complementarity is smaller than this value. %Changing this tolerance to a smaller value may result in greater numerical precision of the solution, but also increases the chance of failure to converge in the algorithm and consequently may result in no solution at all. Therefore, caution is advised in deviating from the default setting. diff --git a/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m b/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m index 4e74b983b6..c867a17a3b 100644 --- a/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m +++ b/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m @@ -219,6 +219,36 @@ %% processing for fluxes [vl,vu,vel,veu,vfl,vfu,vrl,vru,ci,ce,cf,cr,g] = processFluxConstraints(model,param); +if param.debug + modelProcessed = model; + modelProcessed.lb(model.SConsistentRxnBool)=vl; + modelProcessed.lb(~model.SConsistentRxnBool)=vel; + modelProcessed.ub(model.SConsistentRxnBool)=vu; + modelProcessed.ub(~model.SConsistentRxnBool)=veu; + solutionLP = optimizeCbModel(modelProcessed); + + switch solutionLP.stat + case 0 + solution = solutionLP; + message = ['solveCobraEP: LP part of EPproblem is infeasible according to solveCobraLP with ' param.solver '.']; + warning(message) + + return + case 2 + solution = solutionLP; + message = ['solveCobraEP: LP part of EPproblem is unbounded according to solveCobraLP with ' param.solver '.']; + warning(message) + + return + case 1 + message =['solveCobraEP: LP part of EPproblem is feasible according to solveCobraLP with ' param.solver '.']; + fprintf('%s\n',message) + otherwise + error('inconclusive solveCobraLP') + end + messages = cellstr(message); + +end %% optionally processing for concentrations %processConcConstraints if contains(lower(param.entropicFBAMethod),'conc') @@ -1096,25 +1126,11 @@ expConeBool = EPproblem.d~=0; nExpCone = nnz(expConeBool); - % - if 1 - mosekParam=param; - mosekParam.printLevel=param.printLevel-1; - solution = solveCobraEP(EPproblem,mosekParam); - else - [verify,method,printLevel,debug,feasTol,optTol,solver,param] =... - getCobraSolverParams('EP',getCobraSolverParamsOptionsForType('EP'),param); - - solution = solveCobraEP(EPproblem,... - 'verify',verify,... - 'method',method,... - 'printLevel',printLevel,... - 'debug',debug,... - 'feasTol',feasTol,... - 'optTol',optTol,... - 'solver',solver,... - param); - end + + solveCobraEPparam=param; + solveCobraEPparam.printLevel=solveCobraEPparam.printLevel-1; + solution = solveCobraEP(EPproblem,solveCobraEPparam); + switch solution.stat case 1 @@ -1606,7 +1622,9 @@ solution.messages = []; end otherwise - solution_optimizeCbModel = optimizeCbModel(model); + %solution = optimizeCbModel(model, osenseStr, minNorm, allowLoops, param) + param.debug=1; + solution_optimizeCbModel = optimizeCbModel(model,'min',[],1,param); switch solution_optimizeCbModel.stat case 0 message = 'entropicFluxBalanceAnalysis: EPproblem is not feasible, because LP part of model is not feasible according to optimizeCbModel.'; diff --git a/src/base/solvers/entropicFBA/mosekParamSetEFBA.m b/src/base/solvers/entropicFBA/mosekParamSetEFBA.m deleted file mode 100644 index c435a5708e..0000000000 --- a/src/base/solvers/entropicFBA/mosekParamSetEFBA.m +++ /dev/null @@ -1,125 +0,0 @@ -function param=mosekParamSetEFBA(param) -%creates a structure of pertinent user defined options for MOSEK -%OUTPUT -%param parameter structure to be passed to the MOSEK solver - -%MSK_IPAR_LOG_PRESOLVE -% Description:Controls amount of output printed by the presolve procedure. A higher level implies that more information is logged. -% Possible Values:Any number between 0 and +inf. -% Default value:1 -param.MSK_IPAR_LOG_PRESOLVE =1; - -%MSK_IPAR_INTPNT_SCALING -% Controls how the problem is scaled before the interior-point optimizer is used. -% Possible Values: -% MSK_SCALING_NONE -% No scaling is performed. -% MSK_SCALING_MODERATE -% A conservative scaling is performed. -% MSK_SCALING_AGGRESSIVE -% A very aggressive scaling is performed. -% MSK_SCALING_FREE -% The optimizer chooses the scaling heuristic. -% Default value: -% MSK_SCALING_FREE -param.MSK_IPAR_INTPNT_SCALING ='MSK_SCALING_FREE'; -%param.MSK_IPAR_INTPNT_SCALING ='MSK_SCALING_NONE'; - -% MSK_IPAR_INTPNT_REGULARIZATION_USE -% Description:Controls whether regularization is allowed. -% Possible Values: MSK_ON Switch the option on. -% MSK_OFF Switch the option off. -% Default value: MSK_ON -param.MSK_IPAR_INTPNT_REGULARIZATION_USE='MSK_ON'; - -%%%%%%%%%%%%% NONLINEAR TERMINATION CRITERIA%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% MSK_DPAR_INTPNT_CO_TOL_DFEAS -%Dual feasibility tolerance used by the interior-point optimizer for conic problems. -%Default:1.0e-8 -%Accepted: [0.0; 1.0] -param.MSK_DPAR_INTPNT_CO_TOL_DFEAS = 1e-10; - -%MSK_DPAR_INTPNT_CO_TOL_PFEAS -%Primal feasibility tolerance used by the interior-point optimizer for conic problems. -%Default: 1.0e-8 -% Accepted: [0.0; 1.0] -param.MSK_DPAR_INTPNT_CO_TOL_PFEAS = 1.0e-10; %was 1e-11 may be too aggressive -RF - -%MSK_DPAR_INTPNT_CO_TOL_REL_GAP -%Relative gap termination tolerance used by the interior-point optimizer for conic problems. -%Default:1.0e-8 -%Accepted: [0.0; 1.0] -param.MSK_DPAR_INTPNT_CO_TOL_REL_GAP = 1.0e-9; %was 1e-11 may be too aggressive -RF - -%useful for ensuring dual feasibility is as good as primal - -% MSK_IPAR_INTPNT_MAX_ITERATIONS -% Controls the maximum number of iterations allowed in the interior-point optimizer. -% Possible Values:Any number between 0 and +inf. -% Default value: 400 -param.MSK_IPAR_INTPNT_MAX_ITERATIONS=400; - - -%%%%%%%%%%%%% NONLINEAR SOLVER INTEGER PARAM %%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% MSK_IPAR_BI_IGNORE_MAX_ITER -% If the parameter MSK_IPAR_INTPNT_BASIS has the value MSK_BI_NO_ERROR and -% the interior-point optimizer has terminated due to maximum number of -% iterations, then basis identification is performed if this parameter has -% the value MSK_ON. -% Possible Values: -% MSK_ON Switch the option on. -% MSK_OFF Switch the option off. -% Default value: -% MSK_OFF -param.MSK_IPAR_BI_IGNORE_MAX_ITER='MSK_OFF'; - - -%%%%%%%%%%% Solution Approach -% MSK_IPAR_INTPNT_SOLVE_FORM -% Controls whether the primal or the dual problem is solved. -% Possible Values: -% MSK_SOLVE_PRIMAL -% The optimizer should solve the primal problem. -% MSK_SOLVE_DUAL -% The optimizer should solve the dual problem. -% MSK_SOLVE_FREE -% The optimizer is free to solve either the primal or the dual problem. -% Default value:MSK_SOLVE_FREE -param.MSK_IPAR_INTPNT_SOLVE_FORM='MSK_SOLVE_FREE'; -%param.MSK_IPAR_INTPNT_SOLVE_FORM='MSK_SOLVE_PRIMAL'; - -%%%%%%% Infeasibility -% MSK_DPAR_INTPNT_TOL_INFEAS -% Controls when the optimizer declares the model primal or dual infeasible. -% A small number means the optimizer gets more conservative about declaring the model infeasible. -% Possible Values:Any number between 0.0 and 1.0. -% Default value: 1.0e-8 -param.MSK_DPAR_INTPNT_TOL_INFEAS=1e-10; - -%%%%%%%%%%%%%%%%%%%%%%OUTPUT%%%%%%%%%%%%%%%%%% -%MSK_IPAR_LOG_INTPNT -% Controls amount of output printed printed by the interior-point optimizer. -%A higher level implies that more information is logged. -% Possible Values: Any number between 0 and +inf. -% Default value: 4 -param.MSK_IPAR_LOG_INTPNT=5; - -%infesibility report -% MSK_IPAR_INFEAS_REPORT_AUTO -%Controls the amount of information presented in an infeasibility report. -% Possible Values: -% MSK_ON -% Switch the option on. -% MSK_OFF -% Switch the option off. -% Default value: -% MSK_OFF -param.MSK_IPAR_INFEAS_REPORT_AUTO='MSK_OFF'; - -% MSK_IPAR_INFEAS_REPORT_LEVEL -% Controls the amount of information presented in an infeasibility report. Higher values imply more information. -% Possible Values:Any number between 0 and +inf. -% Default value: 1 -%Higher values imply more information. -param.MSK_IPAR_INFEAS_REPORT_LEVEL=100; diff --git a/src/base/solvers/entropicFBA/processConcConstraints.m b/src/base/solvers/entropicFBA/processConcConstraints.m index 57caf4398c..bef526d96e 100644 --- a/src/base/solvers/entropicFBA/processConcConstraints.m +++ b/src/base/solvers/entropicFBA/processConcConstraints.m @@ -27,8 +27,6 @@ % model.dcu: m x 1 real valued upper bound on difference between final and initial initial molecular concentrations (default inf) % model.gasConstant: scalar gas constant (default 8.31446261815324 J K^-1 mol^-1) % model.temperature: scalar temperature (default 310.15 Kelvin) -% -% param.method: 'fluxConc' % param.maxConc: (1e4) maximim micromolar concentration allowed % param.maxConc: (1e-4) minimum micromolar concentration allowed % param.externalNetFluxBounds: ('original') = @@ -87,8 +85,13 @@ nMetabolitesPerRxn = sum(model.S~=0,1)'; bool = nMetabolitesPerRxn>1 & ~model.SConsistentRxnBool; if any(bool) - warning('Exchange reactions involving more than one metabolite, check bounds on x - x0') - disp(model.rxns(bool)) + warning([ int2str(nnz(bool)) ' stoichiometrically inconsistent reactions involving more than one metabolite, check bounds on x - x0']) + if nnz(bool)>10 + ind=find(bool); + disp(model.rxns(ind(1:10))) + else + disp(model.rxns(bool)) + end end if any(~model.SConsistentRxnBool) diff --git a/src/base/solvers/entropicFBA/processFluxConstraints.m b/src/base/solvers/entropicFBA/processFluxConstraints.m index 65a61fb22d..c0944f5764 100644 --- a/src/base/solvers/entropicFBA/processFluxConstraints.m +++ b/src/base/solvers/entropicFBA/processFluxConstraints.m @@ -26,7 +26,7 @@ % param.printLevel: % param.solver: {'pdco',('mosek')} % param.debug: {(0),1} 1 = run in debug mode -% param.method: {('fluxes'),'fluxesConcentrations'} maximise entropy of fluxes (default) or also concentrations +% param.entropicFBAMethod: {('fluxes'),'fluxesConcentrations'} maximise entropy of fluxes (default) or also concentrations % param.maxUnidirectionalFlux: maximum unidirectional flux (1e5 by default) % param.minUnidirectionalFlux: minimum unidirectional flux (zero by default) % param.internalNetFluxBounds: ('original') = use model.lb and model.ub to set the direction and magnitude of internal net flux bounds @@ -74,8 +74,13 @@ param.externalNetFluxBounds='original'; end -if ~isfield(param,'method') - param.method='fluxes'; +if ~isfield(param,'entropicFBAMethod') + if isfield(param,'method') && contains(param.method,'flux') + param.entropicFBAMethod=param.method; + param=rmfield(param,'method'); + else + param.entropicFBAMethod='fluxes'; + end end %find the maximal set of metabolites and reactions that are stoichiometrically consistent @@ -353,7 +358,7 @@ end if ~isfield(model,'g') || isempty(model.g) - if isequal(param.method,'fluxes') + if isequal(param.entropicFBAMethod,'fluxes') model.g='one'; else model.g='two'; diff --git a/src/base/solvers/entropicFBA/solveCobraEP.m b/src/base/solvers/entropicFBA/solveCobraEP.m index 0f008254d1..1a091ce0b1 100644 --- a/src/base/solvers/entropicFBA/solveCobraEP.m +++ b/src/base/solvers/entropicFBA/solveCobraEP.m @@ -1,4 +1,4 @@ -function solution = solveCobraEP(EPproblem, varargin) +function sol = solveCobraEP(EPproblem, varargin) % Solves the following optimisation problem: % % minimize osense*(c.*d)'x + d.*x'(log(x) -1) + (1/2)*x'*Q*x @@ -30,7 +30,7 @@ % % % USAGE: -% solution = solveCobraEP(EPproblem, varargin) +% sol = solveCobraEP(EPproblem, varargin) % % INPUT: % EPproblem: Structure containing the following fields describing the EP problem to be solved @@ -92,37 +92,37 @@ % optTol: Optimality tolerance % % OUTPUT: -% solution: Structure containing the following fields describing a LP solution: -% * .obj: Objective value +% sol: Structure containing the following fields describing a LP sol: +% *.obj: Objective value % *.objLinear osense*c'*x; % *.objEntropy d.*x'*(log(x) -1); % *.objQuadratic (1/2)*x'*Q*x; -% * .full: Primal solution vector +% * .full: Primal sol vector % * .slack: bl = A*x + s = bu -% * .rcost: Reduced costs, dual solution to :math:`lb <= x <= ub` -% * .dual: dual solution to constraints :math: `A*x ('E' | 'G' | 'L') b` +% * .rcost: Reduced costs, dual sol to :math:`lb <= x <= ub` +% * .dual: dual sol to constraints :math: `A*x ('E' | 'G' | 'L') b` % % * .solver: Solver used to solve EP problem % * .stat: Solver status in standardized form % * 0 - Infeasible problem -% * 1 - Optimal solution -% * 2 - Unbounded solution -% * 3 - Almost optimal solution +% * 1 - Optimal sol +% * 2 - Unbounded sol +% * 3 - Almost optimal sol % * -1 - Some other problem (timelimit, numerical problem etc) % * .origStat: Original status returned by the specific solver % * .origStatText: Original status text returned by the specific solver % * .time: Solve time in seconds % % OPTIONAL OUTPUT (from conic optimisation with mosek): -% solution.auxPrimal: auxiliary primal variable -% solution.auxRcost: dual to auxiliary primal variable -% solution.coneF: affine constraint matrix -% solution.coneDual: dual to affine constraints -% solution.dualNorm: dual to the probability normalisation constraint +% sol.auxPrimal: auxiliary primal variable +% sol.auxRcost: dual to auxiliary primal variable +% sol.coneF: affine constraint matrix +% sol.coneDual: dual to affine constraints +% sol.dualNorm: dual to the probability normalisation constraint % % OPTIONAL OUTPUT (from optimisation with pdco): -% solution.d1: primal regularisation parameter, see pdco.m -% solution.d2: dual regularisation parameter, see pdco.m +% sol.d1: primal regularisation parameter, see pdco.m +% sol.d2: dual regularisation parameter, see pdco.m % % EXAMPLE: % @@ -130,15 +130,14 @@ % % Author(s): Ronan M.T. Fleming, 2021 -[problemTypeParams, solverOnlyParams] = parseSolverParameters('EP', varargin{:}); -if ~isfield(problemTypeParams,'debug') - problemTypeParams.debug = 1; -end +[problemTypeParams, solverParams] = parseSolverParameters('EP', varargin{:}); +param = mergeCobraParams(solverParams,problemTypeParams); +clear problemTypeParams solverParams -% Remove outer function specific parameters to avoid crashing solver interfaces -% Default EP parameters are removed within solveCobraEP, so are not removed here -solverOnlyParams = mosekParamStrip(solverOnlyParams); +if ~isfield(param,'debug') + param.debug = 1; +end if any(EPproblem.lb>EPproblem.ub) error('EPproblem.lb>EPproblem.ub'); @@ -153,7 +152,7 @@ end end -if isequal(problemTypeParams.solver,'mosek') +if isequal(param.solver,'mosek') if ~(isfield(EPproblem,'blc') || isfield(EPproblem,'blc')) % blc <= A*x <= buc EPproblem.blc = EPproblem.b; @@ -167,53 +166,38 @@ end %% if in debug mode, test to see if the LP part of the problem is feasible -if problemTypeParams.debug - switch problemTypeParams.solver +if param.debug + switch param.solver case 'pdco' solutionLP2 = solveCobraLP(EPproblem); - if problemTypeParams.printLevel>2 + if param.printLevel>2 disp(solutionLP2) end case 'mosek' - %https://docs.mosek.com/8.1/toolbox/solving-linear.html - if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_PFEAS') - solverOnlyParams.MSK_DPAR_INTPNT_TOL_PFEAS=problemTypeParams.feasTol; - end - if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_DFEAS.') - solverOnlyParams.MSK_DPAR_INTPNT_TOL_DFEAS=problemTypeParams.feasTol; - end - - %remove any fields with names that do not begin with 'MSK_' - solverOnlyParams = mosekParamStrip(solverOnlyParams); + [cmd,mosekParam] = setMosekParam(param); - [res] = msklpopt(EPproblem.c,EPproblem.A,EPproblem.blc,EPproblem.buc,EPproblem.lb,EPproblem.ub,solverOnlyParams,'minimize'); + [res] = msklpopt(EPproblem.c,EPproblem.A,EPproblem.blc,EPproblem.buc,EPproblem.lb,EPproblem.ub,mosekParam,'minimize'); - %If the feasibility tolerance is changed by the solverParams - %struct, this needs to be forwarded to the cobra Params for the - %final consistency test! - if isfield(problemTypeParams,'MSK_DPAR_INTPNT_TOL_PFEAS') - solverOnlyParams.feasTol = solverOnlyParams.MSK_DPAR_INTPNT_TOL_PFEAS; - end - %parse mosek result structure - [solutionLP2.stat,solutionLP2.origStat,x,y,yl,yu,z,zl,zu,s,k,bas,pobjval,dobjval] = parseMskResult(res,EPproblem,solverOnlyParams,problemTypeParams.printLevel); + %[ stat, origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res) + [solutionLP2.stat,solutionLP2.origStat,x,y,yl,yu,z,zl,zu,k, basis,pobjval,dobjval] = parseMskResult(res); switch solutionLP2.stat case 0 - solution = solutionLP2; - message = ['solveCobraEP: LP part of EPproblem is infeasible according to solveCobraLP with ' problemTypeParams.solver '.']; + sol = solutionLP2; + message = ['solveCobraEP: LP part of EPproblem is infeasible according to solveCobraLP with ' param.solver '.']; warning(message) return case 2 - solution = solutionLP2; - message = ['solveCobraEP: LP part of EPproblem is unbounded according to solveCobraLP with ' problemTypeParams.solver '.']; + sol = solutionLP2; + message = ['solveCobraEP: LP part of EPproblem is unbounded according to solveCobraLP with ' param.solver '.']; warning(message) return case 1 - message =['solveCobraEP: LP part of EPproblem is feasible according to solveCobraLP with ' problemTypeParams.solver '.']; + message =['solveCobraEP: LP part of EPproblem is feasible according to solveCobraLP with ' param.solver '.']; fprintf('%s\n',message) otherwise error('inconclusive solveCobraLP') @@ -277,7 +261,7 @@ [mlt,nlt]=size(A); -switch problemTypeParams.solver +switch param.solver case 'pdco' % solves optimization problems of the form % @@ -334,38 +318,38 @@ end end - if isfield(solverOnlyParams,'d1') - d1 = solverOnlyParams.d1; + if isfield(param,'d1') + d1 = param.d1; else d1 = 1e-4; end - if isfield(solverOnlyParams,'d2') - d2 = solverOnlyParams.d2; + if isfield(param,'d2') + d2 = param.d2; else d2 = 1e-4; end - if isfield(solverOnlyParams,'x0') - x0 = solverOnlyParams.x0; + if isfield(param,'x0') + x0 = param.x0; else x0 = ones(size(Aeq,2),1); end - if isfield(solverOnlyParams,'y0') - y0 = solverOnlyParams.y0; + if isfield(param,'y0') + y0 = param.y0; else y0 = ones(size(Aeq,1),1); end - if isfield(solverOnlyParams,'z0') - z0 = solverOnlyParams.z0; + if isfield(param,'z0') + z0 = param.z0; else z0 = ones(size(Aeq,2),1); end - if isfield(solverOnlyParams,'xsize') - xsize = solverOnlyParams.xsize; + if isfield(param,'xsize') + xsize = param.xsize; else xsize = 1; end - if isfield(solverOnlyParams,'zsize') - zsize = solverOnlyParams.zsize; + if isfield(param,'zsize') + zsize = param.zsize; else zsize = 1; end @@ -374,8 +358,8 @@ options = pdcoSet; %options.mu0 = 1; %very small only for entropy function options.mu0 = 0; %pdco chooses its own - options.FeaTol = problemTypeParams.feasTol; - options.OptTol = problemTypeParams.optTol; + options.FeaTol = param.feasTol; + options.OptTol = param.optTol; % If getting linesearch failures, slacken tolerances % i.e. Linesearch failed (nf too big) %options.FeaTol = 1e-6; %%Ecoli core working at 1e-7 @@ -443,7 +427,7 @@ else objHandle = @(x) entropyObj(x,ceq,deq); end - options.Print = problemTypeParams.printLevel-1; + options.Print = param.printLevel-1; saveAndDebug=0; if saveAndDebug @@ -457,7 +441,7 @@ logx(deq~=0) = reallog(x(deq~=0)); % error if negative grad = ceq + deq.*logx; - if problemTypeParams.printLevel > 2 || problemTypeParams.debug + if param.printLevel > 2 || param.debug % determine the residuals fprintf('\n%s\n','KKT with pdco signs:') fprintf('%8.2g %s\n',norm(Aeq*x - beq,inf), '|| Aeq*x - beq ||_inf'); @@ -473,15 +457,16 @@ return end - solution.time = toc; - - % inform = 0 if a solution is found; + sol.time = toc; + sol.epmethod = options.Method; + + % inform = 0 if a sol is found; % = 1 if too many iterations were required; % = 2 if the linesearch failed too often; % = 3 if the step lengths became too small; % = 4 if Cholesky said ADDA was not positive definite. if (inform == 0) - solution.stat = 1; + sol.stat = 1; if ~any(csense == 'L' | csense == 'G') @@ -495,71 +480,73 @@ y(csense == 'G') = -y(csense == 'G'); end - solution.slack = slack; - solution.full = x(1:nlt,1); - solution.dual = -y; - solution.rcost = -z(1:nlt,1); - solution.origStat = inform; + sol.slack = slack; + sol.full = x(1:nlt,1); + sol.dual = -y; + sol.rcost = -z(1:nlt,1); + sol.origStat = inform; %objective logx = zeros(size(A,2),1); - logx(d~=0) = reallog(solution.full(d~=0)); % error if negative + logx(d~=0) = reallog(sol.full(d~=0)); % error if negative if isfield(EPproblem,'Q') - solution.obj = c'*solution.full + (d.*solution.full)'*logx + (1/2)*solution.full'*EPproblem.Q*solution.full; - grad = c + d.*logx + EPproblem.Q*solution.full; + sol.obj = c'*sol.full + (d.*sol.full)'*logx + (1/2)*sol.full'*EPproblem.Q*sol.full; + grad = c + d.*logx + EPproblem.Q*sol.full; else - solution.obj = c'*solution.full + (d.*solution.full)'*logx; + sol.obj = c'*sol.full + (d.*sol.full)'*logx; grad = c + d.*logx; end Aty = -A'*y; - if problemTypeParams.printLevel > 2 || problemTypeParams.debug + if param.printLevel > 2 || param.debug fprintf('\n%s\n','KKT with Rockafellar signs:') - fprintf('%8.2g %s\n',norm(A*solution.full + solution.slack - b,inf), '|| A*x + s - b ||_inf'); - fprintf('%8.2g %s\n',norm(A*solution.full + solution.slack - b - (d2^2)*solution.dual,inf), '|| A*x + s - b - (d2^2)*y ||_inf'); - res2 = grad + A'*solution.dual + solution.rcost; + fprintf('%8.2g %s\n',norm(A*sol.full + sol.slack - b,inf), '|| A*x + s - b ||_inf'); + fprintf('%8.2g %s\n',norm(A*sol.full + sol.slack - b - (d2^2)*sol.dual,inf), '|| A*x + s - b - (d2^2)*y ||_inf'); + res2 = grad + A'*sol.dual + sol.rcost; fprintf('%8.2g %s\n',norm(res2,inf), '|| grad + A''*y + z ||_inf'); - fprintf('%8.2g %s\n',norm(grad + A'*solution.dual + solution.rcost - (d1^2)*solution.full,inf), '|| grad + A''*y + z - (d1^2)*x ||_inf'); - if problemTypeParams.debug - %res2 = grad + Aty + solution.rcost; - res2 = grad + A'*solution.dual + solution.rcost; - solution.T = table(res2,c,d.*logx, -A'*y, solution.rcost,... + fprintf('%8.2g %s\n',norm(grad + A'*sol.dual + sol.rcost - (d1^2)*sol.full,inf), '|| grad + A''*y + z - (d1^2)*x ||_inf'); + if param.debug + %res2 = grad + Aty + sol.rcost; + res2 = grad + A'*sol.dual + sol.rcost; + sol.T = table(res2,c,d.*logx, -A'*y, sol.rcost,... 'VariableNames',{'total','c','dlogx','Aty','z'}); end if any(~isfinite(res2)) warning('Infinite variables in dual optimality condition') - solution.Tinf=solution.T(~isfinite(solution.T.total),:); - ind = find(~isfinite(solution.T.total)); - solution.Tinf.ind = ind; + sol.Tinf=sol.T(~isfinite(sol.T.total),:); + ind = find(~isfinite(sol.T.total)); + sol.Tinf.ind = ind; end end - solution.d1=d1; - solution.d2=d2; + sol.d1=d1; + sol.d2=d2; elseif (inform == 1 || inform == 2 || inform == 3) - solution.stat = 0; - solution.obj = NaN; + sol.stat = 0; + sol.obj = NaN; else - solution.stat = -1; - solution.obj = NaN; + sol.stat = -1; + sol.obj = NaN; end - solution.origStat = inform; + sol.origStat = inform; %update parameters for testing optimality criterion - problemTypeParams.feasTol = options.FeaTol; - problemTypeParams.optTol = options.OptTol; - % * .full: Full LP solution vector + param.feasTol = options.FeaTol; + param.optTol = options.OptTol; + % * .full: Full LP sol vector % * .obj: Objective value - % * .rcost: Reduced costs, dual solution to :math:`lb <= v <= ub` - % * .dual: dual solution to `A*v ('E' | 'G' | 'L') b` + % * .rcost: Reduced costs, dual sol to :math:`lb <= v <= ub` + % * .dual: dual sol to `A*v ('E' | 'G' | 'L') b` % * .solver: Solver used to solve LP problem % * .algorithm: Algorithm used by solver to solve LP problem % * .stat: Solver status in standardized form case 'mosek' %% + [cmd, mosekParam] = setMosekParam(param); + % https://docs.mosek.com/modeling-cookbook/expo.html % https://docs.mosek.com/modeling-cookbook/qcqo.html#conic-reformulation % min (d.*x)'*(log(x./y) + c) + (1/2)*x'*Q*x @@ -669,8 +656,8 @@ % Specify conic part of the problem % https://docs.mosek.com/9.2/toolbox/data-types.html#cones - if problemTypeParams.printLevel>1 || problemTypeParams.debug - [~, res] = mosekopt('symbcon'); + if param.printLevel>1 || param.debug + [~, res] = mosekopt('symbcon',mosekParam); else [~, res] = mosekopt('symbcon echo(0)'); end @@ -809,57 +796,7 @@ mosekopt('write(problem.opf)',prob) end - %set default mosek parameters for this type of problem - paramMosek=mosekParamSetEFBA; - - if ~isfield(solverOnlyParams,'MSK_DPAR_INTPNT_CO_TOL_PFEAS') - if isfield(solverOnlyParams,'MSK_DPAR_INTPNT_CO_TOL_PFEAS') - paramMosek.MSK_DPAR_INTPNT_CO_TOL_PFEAS = solverOnlyParams.feasTol; - else - paramMosek.MSK_DPAR_INTPNT_CO_TOL_PFEAS = problemTypeParams.feasTol; - end - end - if ~isfield(solverOnlyParams,'MSK_DPAR_INTPNT_CO_TOL_DFEAS') - if isfield(solverOnlyParams,'MSK_DPAR_INTPNT_CO_TOL_DFEAS') - paramMosek.MSK_DPAR_INTPNT_CO_TOL_DFEAS = solverOnlyParams.optTol; - else - paramMosek.MSK_DPAR_INTPNT_CO_TOL_DFEAS = problemTypeParams.optTol; - end - end - - % only set the print level if not already set via solverParams structure - if ~isfield(solverOnlyParams, 'MSK_IPAR_LOG') - switch problemTypeParams.printLevel - case 0 - echolev = 0; - case 1 - echolev = 3; - case 2 - paramMosek.MSK_IPAR_LOG_INTPNT = 1; - paramMosek.MSK_IPAR_LOG_SIM = 1; - echolev = 3; - otherwise - echolev = 0; - end - if echolev == 0 - paramMosek.MSK_IPAR_LOG = 0; - cmd = ['minimize echo(' int2str(echolev) ')']; - else - cmd = 'minimize'; - end - - end - %overide if in debug mode - if problemTypeParams.debug - cmd = 'minimize'; - end - - if problemTypeParams.debug && 0 - probBeforeMosekopt = prob; - save('probBeforeMosekopt','probBeforeMosekopt'); - end - %param = updateStructData(param,solverParams); %call mosek exponential cone solver tic; @@ -867,17 +804,21 @@ %default [~,res]=mosekopt('minimize',prob); else - [~,res]=mosekopt(cmd,prob,paramMosek); + [~,res]=mosekopt(cmd,prob,mosekParam); end - solution.time = toc; - + sol.time = toc; + if isfield(param,'epmethod') + sol.epmethod = param.epmethod; + else + sol.epmethod = 'FREE'; + end + %parse mosek result structure - %[stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res,solverParams,printLevel) - [stat,origStat,x,y,yl,yu,z,zl,zu,k,bas,pobjval,dobjval] = parseMskResult(res,solverOnlyParams,problemTypeParams.printLevel); - %[stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res,solverOnlyParams,printLevel) + [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res); - solution.stat = stat; - solution.origStat = origStat; + sol.stat = stat; + sol.origStat = origStat; + sol.obj = pobjval; switch stat case 1 %check for zeros in variables within entropy functions @@ -895,20 +836,20 @@ sbl = prob.a*x - prob.blc; sbu = prob.buc - prob.a*x; s = sbu - sbl; %TODO -double check this - if problemTypeParams.printLevel>1 + if param.printLevel>1 fprintf('%8.2g %s\n',min(sbl), ' min(sbl) = min(A*x - bl), (should be positive)'); fprintf('%8.2g %s\n',min(sbu), ' min(sbu) = min(bu - A*x), (should be positive)'); end - if problemTypeParams.printLevel > 1 + if param.printLevel > 1 % Problem definition here: https://docs.mosek.com/9.2/toolbox/prob-def-affine-conic.html fprintf('%s\n','Optimality conditions (numerical)') - % Guide to interpreting the solution summary: https://docs.mosek.com/9.2/toolbox/debugging-log.html#continuous-problem + % Guide to interpreting the sol summary: https://docs.mosek.com/9.2/toolbox/debugging-log.html#continuous-problem fprintf('%8.2g %s\n',norm(prob.a(prob.blc==prob.buc,:)*x - prob.blc(prob.blc==prob.buc),inf), '|| A*x - b ||_inf'); val = norm(prob.c - prob.a'*y - z - prob.f'*k,inf); fprintf('%8.2g %s\n',val, '|| c - A''*y - z - F''*k ||_inf'); - if val>1e-6 || problemTypeParams.debug - solution.T0 = table(prob.c - prob.a'*y - z - prob.f'*k,prob.c, prob.a'*y, z,prob.f'*k,'VariableNames',{'tot','c','Aty','z','Ftdoty'}); + if val>1e-6 || param.debug + sol.T0 = table(prob.c - prob.a'*y - z - prob.f'*k,prob.c, prob.a'*y, z,prob.f'*k,'VariableNames',{'tot','c','Aty','z','Ftdoty'}); end %fprintf('%8.2g %s\n',norm(prob.c - prob.f'*s,inf), '|| c - F''k ||_inf'); fprintf('%8.2g %s\n',norm(-y + res.sol.itr.slc - res.sol.itr.suc,inf), '|| -y + res.sol.itr.slc - res.sol.itr.suc ||_inf'); @@ -926,14 +867,14 @@ %check with the original order of the affine cone constraints val = norm(prob.c - prob.a'*y - z - F'*y_K,inf); - if problemTypeParams.printLevel > 1 + if param.printLevel > 1 fprintf('%8.2g %s\n',val, '|| c - A''*y - z - F''*y_K ||_inf'); end - if val>1e-6 || problemTypeParams.debug - solution.T = table(prob.c - prob.a'*y - z - F'*y_K,prob.c, prob.a'*y, z,prob.f'*k,F'*y_K,'VariableNames',{'tot','c','Aty','z','Ftdoty','Fty_K'}); + if val>1e-6 || param.debug + sol.T = table(prob.c - prob.a'*y - z - F'*y_K,prob.c, prob.a'*y, z,prob.f'*k,F'*y_K,'VariableNames',{'tot','c','Aty','z','Ftdoty','Fty_K'}); end - if problemTypeParams.printLevel > 1 + if param.printLevel > 1 x1 = F(1:nCone,:)*x; x2 = F(nCone+1:2*nCone,:)*x; x3 = F(2*nCone+1:3*nCone,:)*x; @@ -973,40 +914,40 @@ fprintf('%7.2g\t%s\n',min(y1_K(1:nExpCone) + y3_K(1:nExpCone).*exp(y2_K(1:nExpCone)./y3_K(1:nExpCone))/exp(1)), 'min(y1_k + y3_k.*exp(y2_K./y3_K)/exp(1)) >= 0'); end - solution.full = x(1:size(A,2)); + sol.full = x(1:size(A,2)); %switch to Rockafellar signs - solution.dual = -y(1:size(A,1)); - solution.dualNorm = -y(size(A,1)+1:size(A,1)+p); - solution.rcost = -z(1:size(A,2)+p); - solution.slack = s; + sol.dual = -y(1:size(A,1)); + sol.dualNorm = -y(size(A,1)+1:size(A,1)+p); + sol.rcost = -z(1:size(A,2)+p); + sol.slack = s; %need to zero out the NaN due to log(0) for some variables - logSolutionFull = real(log(solution.full)); + logSolutionFull = real(log(sol.full)); logSolutionFull(~isfinite(logSolutionFull))=0; if isfield(EPproblem,'Q') - solution.obj = EPproblem.c'*solution.full + (EPproblem.d.*solution.full)'*(logSolutionFull -1) + (1/2)*solution.full'*EPproblem.Q*solution.full; - solution.objLinear = EPproblem.c'*solution.full; - solution.objEntropy = -(EPproblem.d.*solution.full)'*(logSolutionFull -1); - solution.objQuadratic = (1/2)*solution.full'*EPproblem.Q*solution.full; + sol.obj = EPproblem.c'*sol.full + (EPproblem.d.*sol.full)'*(logSolutionFull -1) + (1/2)*sol.full'*EPproblem.Q*sol.full; + sol.objLinear = EPproblem.c'*sol.full; + sol.objEntropy = -(EPproblem.d.*sol.full)'*(logSolutionFull -1); + sol.objQuadratic = (1/2)*sol.full'*EPproblem.Q*sol.full; else - solution.obj = EPproblem.c'*solution.full + (EPproblem.d.*solution.full)'*(logSolutionFull -1); - solution.objLinear = EPproblem.c'*solution.full; - solution.objEntropy = -(EPproblem.d.*solution.full)'*(logSolutionFull -1); - solution.objQuadratic = 0; + sol.obj = EPproblem.c'*sol.full + (EPproblem.d.*sol.full)'*(logSolutionFull -1); + sol.objLinear = EPproblem.c'*sol.full; + sol.objEntropy = -(EPproblem.d.*sol.full)'*(logSolutionFull -1); + sol.objQuadratic = 0; end - posRcost = solution.rcost>0; - negRcost = solution.rcost<0; + posRcost = sol.rcost>0; + negRcost = sol.rcost<0; blx = prob.blx(1:size(A,2)); bux = prob.bux(1:size(A,2)); - solution.lagRcost = sum(solution.rcost(negRcost)'*blx(negRcost) + solution.rcost(posRcost)'*bux(posRcost)); + sol.lagRcost = sum(sol.rcost(negRcost)'*blx(negRcost) + sol.rcost(posRcost)'*bux(posRcost)); %pass back the F matrix to check biochemical optimality criteria - solution.coneF = F; - solution.auxPrimal = x(size(A,2)+p+1:end); - solution.auxRcost = -z(size(A,2)+p+1:end); - solution.coneDual = -y_K; + sol.coneF = F; + sol.auxPrimal = x(size(A,2)+p+1:end); + sol.auxRcost = -z(size(A,2)+p+1:end); + sol.coneDual = -y_K; % variable to determine the residual 1 b = prob.blc; @@ -1021,146 +962,111 @@ k = NaN*ones(size(prob.f,1),1); end otherwise - error([problemTypeParams.solver ' is an unrecognised solver']) + error([param.solver ' is an unrecognised solver']) end -switch solution.stat +sol.solver=param.solver; + +switch sol.stat case 0 - switch problemTypeParams.solver + switch param.solver case 'pdco' %infeasible, debug the situtation - disp(solution.origStat) - %solution.origStat: 'PRIMAL_INFEASIBLE_CER' + disp(sol.origStat) + %sol.origStat: 'PRIMAL_INFEASIBLE_CER' solutionLP = solveCobraLP(EPproblem); statLP=solutionLP.stat; case 'mosek' - %https://docs.mosek.com/8.1/toolbox/solving-linear.html - if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_PFEAS') - solverOnlyParams.MSK_DPAR_INTPNT_TOL_PFEAS=problemTypeParams.feasTol; - end - if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_DFEAS.') - solverOnlyParams.MSK_DPAR_INTPNT_TOL_DFEAS=problemTypeParams.feasTol; - end - %If the feasibility tolerance is changed by the solverParams - %struct, this needs to be forwarded to the cobra Params for the - %final consistency test! - if isfield(problemTypeParams,'MSK_DPAR_INTPNT_TOL_PFEAS') - solverOnlyParams.feasTol = solverOnlyParams.MSK_DPAR_INTPNT_TOL_PFEAS; - end + [cmd,mosekParam] = setMosekParam(param); - - % only set the print level if not already set via solverParams structure - if ~isfield(solverOnlyParams, 'MSK_IPAR_LOG') - switch problemTypeParams.printLevel - case 0 - echolev = 0; - case 1 - echolev = 3; - case 2 - solverOnlyParams.MSK_IPAR_LOG_INTPNT = 1; - solverOnlyParams.MSK_IPAR_LOG_SIM = 1; - echolev = 3; - otherwise - echolev = 0; - end - end - if echolev == 0 - solverOnlyParams.MSK_IPAR_LOG = 0; - cmd = ['minimize echo(' int2str(echolev) ')']; - else - cmd = 'minimize'; - end - [res] = msklpopt(EPproblem.c,EPproblem.A,EPproblem.blc,EPproblem.buc,EPproblem.lb,EPproblem.ub,solverOnlyParams,cmd); - - %[stat,origStat,x,y,yl,yu,z,zl,zu,s,k,bas,pobjval,dobjval] = parseMskResult(res,prob,problemTypeParams.printLevel); - [statLP,origStat,x,y,yl,yu,z,zl,zu,k,bas,pobjval,dobjval] = parseMskResult(res,solverOnlyParams,problemTypeParams.printLevel); - - + [res] = msklpopt(EPproblem.c,EPproblem.A,EPproblem.blc,EPproblem.buc,EPproblem.lb,EPproblem.ub,mosekParam,cmd); + + [statLP,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res); end switch statLP case 1 - message =['solveCobraEP: EPproblem with ' problemTypeParams.solver ' is infeasible, but corresponding LPproblem is feasible according to solveCobraLP with ' problemTypeParams.solver]; + message =['solveCobraEP: EPproblem with ' param.solver ' is infeasible, but corresponding LPproblem is feasible according to solveCobraLP with ' param.solver]; warning(message) otherwise - message = ['solveCobraEP: EPproblem with ' problemTypeParams.solver ' is infeasible, because corresponding LPproblem is infeasible according to solveCobraLP with ' problemTypeParams.solver]; + message = ['solveCobraEP: EPproblem with ' param.solver ' is infeasible, because corresponding LPproblem is infeasible according to solveCobraLP with ' param.solver]; warning(message) end if exist('messages','var') - if isfield(solution,'messages') - solution.messages = [messages;solution.messages;message]; + if isfield(sol,'messages') + sol.messages = [messages;sol.messages;message]; else - solution.messages = [messages;message]; + sol.messages = [messages;message]; end else - solution.messages = cellstr(message); + sol.messages = cellstr(message); end case 1 % check the optimality conditions for various solvers - if ~isempty(solution.slack) && ~isempty(solution.full) + if ~isempty(sol.slack) && ~isempty(sol.full) % determine the residual 1 - switch problemTypeParams.solver + switch param.solver case 'pdco' feasTol = 1e-3; - res1 = A*solution.full + solution.slack - b; + res1 = A*sol.full + sol.slack - b; res1(~isfinite(res1))=0; case 'mosek' - feasTol = problemTypeParams.feasTol * 1e2; - res1 = A(blc==buc,:)*solution.full - blc(blc==buc); + feasTol = param.feasTol * 1e2; + res1 = A(blc==buc,:)*sol.full - blc(blc==buc); end tmp1 = norm(res1, inf); % evaluate the optimality condition 1 if tmp1 > feasTol - if strcmp(problemTypeParams.solver,'pdco') - res1b = norm(A*solution.full + solution.slack - b + (d2^2)*y,inf); + if strcmp(param.solver,'pdco') + res1b = norm(A*sol.full + sol.slack - b + (d2^2)*y,inf); tmp1b = norm(res1b, inf); if tmp1b > feasTol displayError = 1; else displayError = 0; - warning(['[' problemTypeParams.solver '] Primal optimality condition in solveCobraEP only approximately satisfied, residual = ' num2str(tmp1) ', regularised residual = ' num2str(tmp1b) ', while problem feasTol = ' num2str(feasTol) '. origStat = ' solution.origStat]) + warning(['[' param.solver '] Primal optimality condition in solveCobraEP only approximately satisfied, residual = ' num2str(tmp1) ', regularised residual = ' num2str(tmp1b) ', while problem feasTol = ' num2str(feasTol) '. origStat = ' sol.origStat]) end else %TODO - debug why solver reporting optimal but unscaled seems less so. displayError = 0; end if displayError - %disp(solution.origStat) - fprintf('%s\n',['[' problemTypeParams.solver '] Primal optimality condition in solveCobraEP not satisfied, residual = ' num2str(tmp1) ', while problem feasTol = ' num2str(feasTol) '. origStat = ' solution.origStat]) + %disp(sol.origStat) + fprintf('%s\n',['[' param.solver '] Primal optimality condition in solveCobraEP not satisfied, residual = ' num2str(tmp1) ', while problem feasTol = ' num2str(feasTol) '. origStat = ' sol.origStat]) end else - if problemTypeParams.printLevel > 0 - fprintf(['\n > [' problemTypeParams.solver '] Primal optimality condition in solveCobraEP satisfied.']); + if param.printLevel > 0 + fprintf(['\n > [' param.solver '] Primal optimality condition in solveCobraEP satisfied.']); end end end %gradient may differ depending on the solver - res2 = grad + Aty + solution.rcost; + res2 = grad + Aty + sol.rcost; tmp2 = norm(res2, inf); if 0 - optTol = problemTypeParams.optTol * 1e2; + optTol = param.optTol * 1e2; else optTol = 5e-5; end % evaluate the optimality condition 2 if tmp2 > optTol - disp(solution.origStat) - if ~(length(A)==1 && strcmp(problemTypeParams.solver,'pdco')) %todo, why does pdco choke on small A? - warning(['[' problemTypeParams.solver '] Dual optimality condition in solveCobraEP not satisfied, residual = ' num2str(tmp2) ', while problem optTol = ' num2str(optTol)]) + disp(sol.origStat) + if ~(length(A)==1 && strcmp(param.solver,'pdco')) %todo, why does pdco choke on small A? + warning(['[' param.solver '] Dual optimality condition in solveCobraEP not satisfied, residual = ' num2str(tmp2) ', while problem optTol = ' num2str(optTol)]) end else - if problemTypeParams.printLevel > 0 - fprintf(['\n > [' problemTypeParams.solver '] Dual optimality condition in solveCobraEP satisfied.\n']); + if param.printLevel > 0 + fprintf(['\n > [' param.solver '] Dual optimality condition in solveCobraEP satisfied.\n']); end end end -solution.solver=problemTypeParams.solver; + end diff --git a/src/base/solvers/getSetSolver/changeCobraSolver.m b/src/base/solvers/getSetSolver/changeCobraSolver.m index 0dfac287dc..19bbc7c041 100644 --- a/src/base/solvers/getSetSolver/changeCobraSolver.m +++ b/src/base/solvers/getSetSolver/changeCobraSolver.m @@ -1,4 +1,4 @@ -function [solverOK, solverInstalled] = changeCobraSolver(solverName, solverType, printLevel, validationLevel) +function [solverOK, solverInstalled] = changeCobraSolver(solverName, problemType, printLevel, validationLevel) % Changes the Cobra Toolbox optimization solver(s) % % USAGE: @@ -7,10 +7,10 @@ % % INPUTS: % solverName: Solver name -% solverType: Solver type ('LP' by default) +% problemType: Problem type ('LP' by default) % (a) One of the following: `LP` `MILP`, `QP`, `MIQP` 'EP', 'CLP' % (b) 'all' attempts to change all applicable solvers to solverName. This is purely a shorthand convenience. -% (c) Cell array of solverTypes, e.g. {'LP','QP'} +% (c) Cell array of problemTypes, e.g. {'LP','QP'} % printLevel: verbose level % % * if `0`, warnings and errors are silenced @@ -21,7 +21,7 @@ % % * `-1`: assign only the global variable. Do not assign any path. % * `0`: adjust solver paths but don't validate the solver -% * `1`: validate but remove outputs (default) +% * `1`: validate but remove outputs, silent (default) % * `2`: validate and keep any outputs % % OUTPUT: @@ -207,7 +207,7 @@ solverInstalled = true; if validationLevel == -1 - switch solverType + switch problemType case 'LP' CBT_LP_SOLVER = solverName; case 'QP' @@ -299,9 +299,9 @@ end if nargin < 2 - solverType = 'LP'; + problemType = 'LP'; else - solverType = upper(solverType); + problemType = upper(problemType); end % print an error message if the solver is not supported @@ -333,12 +333,12 @@ end % Attempt to set the user provided solver for all optimization problem types -if strcmpi(solverType, 'all') +if strcmpi(problemType, 'all') solvedProblems = SOLVERS.(solverName).type; for i = 1:length(solvedProblems) %this looks to be self referential (calling changeCobraSolver in %changeCobraSolver) - [solverOK,solverInstalled] = changeCobraSolver(solverName, solvedProblems{i}, printLevel); + [solverOK,solverInstalled] = changeCobraSolver(solverName, solvedProblems{i}, printLevel,validationLevel); if printLevel > 0 fprintf([' > changeCobraSolver: Solver for ', solvedProblems{i}, ' problems has been set to ', solverName, '.\n']); end @@ -360,24 +360,24 @@ % check if the given solver is able to solve the given problem type. solverOK = false; -if ~contains(solverType, OPT_PROB_TYPES) +if ~contains(problemType, OPT_PROB_TYPES) %This is not done during init, so at this point, the solver is already %checked for installation solverInstalled = SOLVERS.(solverName).installed; if printLevel > 0 - error('changeCobraSolver: %s problems cannot be solved in The COBRA Toolbox', solverType); + error('changeCobraSolver: %s problems cannot be solved in The COBRA Toolbox', problemType); else return end end % check if the given solver is able to solve the given problem type. -if ~contains(solverType, SOLVERS.(solverName).type) +if ~contains(problemType, SOLVERS.(solverName).type) %This is not done during init, so at this point, the solver is already %checked for installation solverInstalled = SOLVERS.(solverName).installed; if printLevel > 0 - error('Solver %s cannot solve %s problems', solverName, solverType); + error('Solver %s cannot solve %s problems', solverName, problemType); else return end @@ -528,10 +528,10 @@ if validationLevel > 0 cwarn = warning; warning('off'); - eval(['oldval = CBT_', solverType, '_SOLVER;']); - eval(['CBT_', solverType, '_SOLVER = solverName;']); + eval(['oldval = CBT_', problemType, '_SOLVER;']); + eval(['CBT_', problemType, '_SOLVER = solverName;']); % validate with a simple problem. - if strcmp(solverName,'mosek') && strcmp(solverType,'CLP') || strcmp(solverType,'all') + if strcmp(solverName,'mosek') && strcmp(problemType,'CLP') || strcmp(problemType,'all') problem = struct('A',[0 1],'b',0,'c',[1;1],'osense',-1,'lb',[0;0],'ub',[0;0],'csense','E','vartype',['C';'I'],'x0',[0;0]); else problem = struct('A',[0 1],'b',0,'c',[1;1],'osense',-1,'F',speye(2),'lb',[0;0],'ub',[0;0],'csense','E','vartype',['C';'I'],'x0',[0;0]); @@ -540,9 +540,9 @@ %This is the code that actually tests if a solver is working if validationLevel>1 %display progress - eval(['solveCobra' solverType '(problem,''printLevel'',3);']); + eval(['solveCobra' problemType '(problem,''printLevel'',3);']); else - eval(['solveCobra' solverType '(problem,''printLevel'',0);']); + eval(['solveCobra' problemType '(problem,''printLevel'',0);']); end catch ME %This is the code that describes what went wrong if a call to a solver does not work @@ -553,12 +553,12 @@ rethrow(ME) end solverOK = false; - eval(['CBT_', solverType, '_SOLVER = oldval;']); + eval(['CBT_', problemType, '_SOLVER = oldval;']); end warning(cwarn) else % if unvalidated, simply set the solver without testing. - eval(['CBT_', solverType, '_SOLVER = solverName;']); + eval(['CBT_', problemType, '_SOLVER = solverName;']); end else switch solverName diff --git a/src/base/solvers/gurobi/setGurobiParam.m b/src/base/solvers/gurobi/setGurobiParam.m new file mode 100644 index 0000000000..f5f47f0e94 --- /dev/null +++ b/src/base/solvers/gurobi/setGurobiParam.m @@ -0,0 +1,342 @@ +function gurobiParam = setGurobiParam(param) +% The params struct contains Gurobi parameters. A full list may be +% found on the Parameter page of the reference manual: +% https://www.gurobi.com/documentation/current/refman/parameter_descriptions.html +% Parameters must be in TimeLimit not timelimit or timeLimit format, see +% below for full list that are eligible to be passed to the solver. + +if isfield(param,'timelimit') + param.TimeLimit=param.timelimit; +end + +if isfield(param,'multiscale') && param.multiscale==1 + param.ScaleFlag=0; +end + +%backward compatibility +if isfield(param,'method') + if isempty(param.method) + param = rmfield(param,'method'); + else + if ~isfield(param,[lower(param.problemType) 'method']) + param.([lower(param.problemType) 'method'])=param.method; + end + end +end + +% https://www.gurobi.com/documentation/current/refman/method.html +% params.method gives the method used to solve continuous models +% -1=automatic, +% 0=primal simplex, +% 1=dual simplex, +% 2=barrier, +% 3=concurrent, +% 4=deterministic concurrent +% i.e. params.method = 1; % use dual simplex method +if isfield(param,'lpmethod') + %gurobiAlgorithms = {'AUTOMATIC','PRIMAL','DUAL','BARRIER','CONCURRENT','CONCURRENT_DETERMINISTIC'}; + % -1=automatic, + % 0=primal simplex, + % 1=dual simplex, + % 2=barrier, + % 3=concurrent, + % 4=deterministic concurrent + switch param.lpmethod + case 'AUTOMATIC' + param.Method = -1; + case 'PRIMAL' + param.Method = 0; + case 'DUAL' + param.Method = 1; + case 'BARRIER' + param.Method = 2; + case 'CONCURRENT' + param.Method = 3; + case 'DETERMINISTIC_CONCURRENT' + param.Method = 4; + otherwise + %https://www.gurobi.com/documentation/current/refman/method.html + %Concurrent methods aren't available for QP and QCP. + warning([param.lpmethod ' is an unrecognised param.qpmethod for gurobi']) + end + param = rmfield(param,'lpmethod'); +end + +% param.qpmethod gives the qpmethod used to solve continuous models +% -1=automatic, +% 0=primal simplex, +% 1=dual simplex, +% 2=barrier, +% 3=concurrent, +% 4=deterministic concurrent +% i.e. param.qpmethod = 1; % use dual simplex method +if isfield(param,'qpmethod') + %gurobi QP algorithms = {'AUTOMATIC','PRIMAL','DUAL','BARRIER'}; + % -1=automatic, + % 0=primal simplex, + % 1=dual simplex, + % 2=barrier, + switch param.qpmethod + case 'AUTOMATIC' + param.Method = -1; + case 'PRIMAL' + param.Method = 0; + case 'DUAL' + param.Method = 1; + case 'BARRIER' + param.Method = 2; + otherwise + %https://www.gurobi.com/documentation/current/refman/method.html + %Concurrent methods aren't available for QP and QCP. + warning([param.qpmethod ' is an unrecognised param.qpmethod for gurobi']) + end + param = rmfield(param,'qpmethod'); +end + +if ~isfield(param,'OutputFlag') + switch param.printLevel + case 0 + param.OutputFlag = 0; + case 1 + param.OutputFlag = 0; + otherwise + % silent + param.OutputFlag = 0; + end +end + +if ~isfield(param,'DisplayInterval') + switch param.printLevel + case 0 + param.DisplayInterval = 1; + case 1 + param.DisplayInterval = 1; + otherwise + % silent + param.DisplayInterval = 1; + end +end + +if ~isfield(param,'FeasibilityTol') + % Primal feasibility tolerance + % Type: double + % Default value: 1e-6 + % Minimum value: 1e-9 + % Maximum value: 1e-2 + % All constraints must be satisfied to a tolerance of FeasibilityTol. + % Tightening this tolerance can produce smaller constraint violations, but for + % numerically challenging models it can sometimes lead to much larger iteration counts. + param.FeasibilityTol = param.feasTol; +end +if ~isfield(param,'OptimalityTol') + % Dual feasibility tolerance + % Type: double + % Default value: 1e-6 + % Minimum value: 1e-9 + % Maximum value: 1e-2 + % Reduced costs must all be smaller than OptimalityTol in the improving direction in order for a model to be declared optimal. + param.OptimalityTol = param.optTol; +end + + +% Permitted parameter fields +permittedFields = {... +'AggFill'... +'Aggregate'... +'BarConvTol'... +'BarCorrectors'... +'BarHomogeneous'... +'BarIterLimit'... +'BarOrder'... +'BarQCPConvTol'... +'BestBdStop'... +'BestObjStop'... +'BQPCuts'... +'BranchDir'... +'CliqueCuts'... +'CloudAccessID'... +'CloudHost'... +'CloudSecretKey'... +'CloudPool'... +'ComputeServer'... +'ConcurrentJobs'... +'ConcurrentMethod'... +'ConcurrentMIP'... +'ConcurrentSettings'... +'CoverCuts'... +'Crossover'... +'CrossoverBasis'... +'CSAPIAccessID'... +'CSAPISecret'... +'CSAppName'... +'CSAuthToken'... +'CSBatchMode'... +'CSClientLog'... +'CSGroup'... +'CSIdleTimeout'... +'CSManager'... +'CSPriority'... +'CSQueueTimeout'... +'CSRouter'... +'CSTLSInsecure'... +'CutAggPasses'... +'Cutoff'... +'CutPasses'... +'Cuts'... +'DegenMoves'... +'Disconnected'... +'DisplayInterval'... +'DistributedMIPJobs'... +'DualReductions'... +'FeasibilityTol'... +'FeasRelaxBigM'... +'FlowCoverCuts'... +'FlowPathCuts'... +'FuncPieceError'... +'FuncPieceLength'... +'FuncPieceRatio'... +'FuncPieces'... +'FuncMaxVal'... +'FuncNonlinear'... +'GomoryPasses'... +'GUBCoverCuts'... +'Heuristics'... +'IgnoreNames'... +'IISMethod'... +'ImpliedCuts'... +'ImproveStartGap'... +'ImproveStartNodes'... +'ImproveStartTime'... +'InfProofCuts'... +'InfUnbdInfo'... +'InputFile'... +'IntegralityFocus'... +'IntFeasTol'... +'IterationLimit'... +'JobID'... +'JSONSolDetail'... +'LazyConstraints'... +'LicenseID'... +'LiftProjectCuts'... +'LPWarmStart'... +'LogFile'... +'LogToConsole'... +'MarkowitzTol'... +'MemLimit'... +'Method'... +'MinRelNodes'... +'MIPFocus'... +'MIPGap'... +'MIPGapAbs'... +'MIPSepCuts'... +'MIQCPMethod'... +'MIRCuts'... +'MixingCuts'... +'ModKCuts'... +'MultiObjMethod'... +'MultiObjPre'... +'MultiObjSettings'... +'NetworkAlg'... +'NetworkCuts'... +'NLPHeur'... +'NodefileDir'... +'NodefileStart'... +'NodeLimit'... +'NodeMethod'... +'NonConvex'... +'NoRelHeurTime'... +'NoRelHeurWork'... +'NormAdjust'... +'NumericFocus'... +'OBBT'... +'ObjNumber'... +'ObjScale'... +'OptimalityTol'... +'OutputFlag'... +'PartitionPlace'... +'PerturbValue'... +'PoolGap'... +'PoolGapAbs'... +'PoolSearchMode'... +'PoolSolutions'... +'PreCrush'... +'PreDepRow'... +'PreDual'... +'PreMIQCPForm'... +'PrePasses'... +'PreQLinearize'... +'Presolve'... +'PreSOS1BigM'... +'PreSOS1Encoding'... +'PreSOS2BigM'... +'PreSOS2Encoding'... +'PreSparsify'... +'ProjImpliedCuts'... +'PSDCuts'... +'PSDTol'... +'PumpPasses'... +'QCPDual'... +'Quad'... +'Record'... +'ResultFile'... +'RINS'... +'RelaxLiftCuts'... +'RLTCuts'... +'ScaleFlag'... +'ScenarioNumber'... +'Seed'... +'ServerPassword'... +'ServerTimeout'... +'Sifting'... +'SiftMethod'... +'SimplexPricing'... +'SoftMemLimit'... +'SolutionLimit'... +'SolutionTarget'... +'SolFiles'... +'SolutionNumber'... +'StartNodeLimit'... +'StartNumber'... +'StrongCGCuts'... +'SubMIPCuts'... +'SubMIPNodes'... +'Symmetry'... +'Threads'... +'TimeLimit'... +'TokenServer'... +'TSPort'... +'TuneBaseSettings'... +'TuneCleanup'... +'TuneCriterion'... +'TuneDynamicJobs'... +'TuneJobs'... +'TuneMetric'... +'TuneOutput'... +'TuneResults'... +'TuneTargetMIPGap'... +'TuneTargetTime'... +'TuneTimeLimit'... +'TuneTrials'... +'TuneUseFilename'... +'UpdateMode'... +'UserName'... +'VarBranch'... +'WLSAccessID'... +'WLSSecret'... +'WLSToken'... +'WLSTokenDuration'... +'WLSTokenRefresh'... +'WorkerPassword'... +'WorkerPool'... +'WorkLimit'... +'ZeroHalfCuts'... +'ZeroObjNodes'}; + +fields = fieldnames(param); +fieldsToRemove = setdiff(fields, permittedFields); +%remove fields that should not be present +gurobiParam = rmfield(param, fieldsToRemove); + +if isfield(gurobiParam,'logFile') && isempty(gurobiParam.logFile) + gurobiParam = rmfield(gurobiParam,'logFile'); +end diff --git a/src/base/solvers/init/configEnvVars.m b/src/base/solvers/init/configEnvVars.m index 89c40b6c92..5d9eb639a2 100644 --- a/src/base/solvers/init/configEnvVars.m +++ b/src/base/solvers/init/configEnvVars.m @@ -191,7 +191,7 @@ end % Method 1: check if the solver is already on the MATLAB path - isOnPath = ~isempty(strfind(lower(path), lower(possibleDir))); + isOnPath = contains(lower(path), lower(possibleDir)); % find the index of the most recently added solver path tmp = path; diff --git a/src/base/solvers/entropicFBA/mosekParamStrip.m b/src/base/solvers/mosek/mosekParamStrip.m similarity index 100% rename from src/base/solvers/entropicFBA/mosekParamStrip.m rename to src/base/solvers/mosek/mosekParamStrip.m diff --git a/src/base/solvers/msk/parseMskResult.m b/src/base/solvers/mosek/parseMskResult.m similarity index 97% rename from src/base/solvers/msk/parseMskResult.m rename to src/base/solvers/mosek/parseMskResult.m index 0a61e1081b..8282b0cd1e 100644 --- a/src/base/solvers/msk/parseMskResult.m +++ b/src/base/solvers/mosek/parseMskResult.m @@ -1,4 +1,4 @@ -function [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res,solverOnlyParams,printLevel) +function [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res) %parse the res structure returned from mosek % INPUTS: % res: mosek results structure returned by mosekopt @@ -54,12 +54,6 @@ pobjval =[]; dobjval =[]; -if ~exist('printLevel','var') - printLevel = 0; -end -if ~exist('solverOnlyParams','var') - solverOnlyParams = struct(); -end % prosta (string) – Problem status (prosta). % prosta @@ -126,7 +120,7 @@ % The primal solution is integer optimal. % https://docs.mosek.com/latest/toolbox/accessing-solution.html -accessSolution=[]; +accessSolution=''; if isfield(res, 'sol') if isfield(res.sol,'itr') && isfield(res.sol,'bas') if any(strcmp(res.sol.bas.solsta,{'OPTIMAL','MSK_SOL_STA_OPTIMAL','MSK_SOL_STA_NEAR_OPTIMAL'})) && any(strcmp(res.sol.itr.solsta,{'UNKNOWN'})) @@ -217,6 +211,8 @@ otherwise accessSolution = 'dontAccess'; end + otherwise + accessSolution = 'dontAccess'; end if strcmp(accessSolution,'dontAccess') @@ -241,3 +237,6 @@ end end end + + +% https://themosekblog.blogspot.com/2014/06/what-if-solver-stall.html \ No newline at end of file diff --git a/src/base/solvers/mosek/setMosekParam.m b/src/base/solvers/mosek/setMosekParam.m new file mode 100644 index 0000000000..d867594482 --- /dev/null +++ b/src/base/solvers/mosek/setMosekParam.m @@ -0,0 +1,373 @@ +function [cmd,mosekParam] = setMosekParam(param) +% set mosek parameters from param fields +% strip any non mosek compatible fields from param and return it as +% mosekParam + + +%tests if solver correctly interfaced and licence running +if param.printLevel>1 || param.debug + [~, res] = mosekopt('symbcon'); +else + [~, res] = mosekopt('symbcon echo(0)'); +end + +% only set the print level if not already set via param structure +if ~isfield(param, 'MSK_IPAR_LOG') + % Controls the amount of log information. + % The value 0 implies that all log information is suppressed. + % A higher level implies that more information is logged. + switch param.printLevel + case 0 + echolev = 0; + case 1 + echolev = 3; + case 2 + param.MSK_IPAR_WRITE_DATA_PARAM='MSK_ON'; + param.MSK_IPAR_LOG_INTPNT = 1; + param.MSK_IPAR_LOG_SIM = 1; + %MSK_IPAR_LOG_PRESOLVE + % Description:Controls amount of output printed by the presolve procedure. A higher level implies that more information is logged. + % Possible Values:Any number between 0 and +inf. + % Default value:1 + param.MSK_IPAR_LOG_PRESOLVE=10; + + %MSK_IPAR_LOG_INTPNT + % Controls amount of output printed printed by the interior-point optimizer. + %A higher level implies that more information is logged. + % Possible Values: Any number between 0 and +inf. + % Default value: 4 + if ~isfield(param,'MSK_IPAR_LOG_INTPNT') + param.MSK_IPAR_LOG_INTPNT=5; + end + + %infesibility report + % MSK_IPAR_INFEAS_REPORT_AUTO + %Controls the amount of information presented in an infeasibility report. + % Possible Values: + % MSK_ON + % Switch the option on. + % MSK_OFF + % Switch the option off. + % Default value: + % MSK_OFF + if ~isfield(param,'MSK_IPAR_INFEAS_REPORT_AUTO') + param.MSK_IPAR_INFEAS_REPORT_AUTO='MSK_ON'; + end + + % MSK_IPAR_INFEAS_REPORT_LEVEL + % Controls the amount of information presented in an infeasibility report. Higher values imply more information. + % Possible Values:Any number between 0 and +inf. + % Default value: 1 + %Higher values imply more information. + if ~isfield(param,'MSK_IPAR_INFEAS_REPORT_LEVEL') + param.MSK_IPAR_INFEAS_REPORT_LEVEL=1; + end + + echolev = 3; + otherwise + echolev = 0; + end + if echolev == 0 && ~param.debug + param.MSK_IPAR_LOG = 0; + cmd = ['minimize echo(' int2str(echolev) ')']; + else + cmd = 'minimize'; + end +end + + +if ~isfield(param, 'MSK_DPAR_OPTIMIZER_MAX_TIME') && isfield(param,'timelimit') + % MSK_DPAR_OPTIMIZER_MAX_TIME + % Maximum amount of time the optimizer is allowed to spent on the optimization. A negative number means infinity. + % Default + % -1.0 + % Accepted + % [-inf; +inf] + % Example + % param.MSK_DPAR_OPTIMIZER_MAX_TIME = -1.0 + % Groups + % Termination criteria + param.MSK_DPAR_OPTIMIZER_MAX_TIME = param.timelimit; +end + + + +if ~isfield(param, 'MSK_DPAR_INTPNT_TOL_PFEAS') + % Primal feasibility tolerance used by the interior-point optimizer for linear problems. + % Default + % 1.0e-8 + % Accepted + % [0.0; 1.0] + % Example + % param.MSK_DPAR_INTPNT_TOL_PFEAS = 1.0e-8 + % Groups + % Interior-point method, Termination criteria + param.MSK_DPAR_INTPNT_TOL_PFEAS=param.feasTol; +end + + +if ~isfield(param,'MSK_DPAR_INTPNT_QO_TOL_PFEAS') + % Primal feasibility tolerance used by the interior-point optimizer for quadratic problems. + % Default + % 1.0e-8 + % Accepted + % [0.0; 1.0] + % Example + % param.MSK_DPAR_INTPNT_QO_TOL_PFEAS = 1.0e-8 + % See also + % MSK_DPAR_INTPNT_QO_TOL_NEAR_REL + % Groups + % Interior-point method, Termination criteria + param.MSK_DPAR_INTPNT_QO_TOL_PFEAS=param.feasTol; +end + +if ~isfield(param, 'MSK_DPAR_INTPNT_CO_TOL_PFEAS') + % Primal feasibility tolerance used by the interior-point optimizer for conic problems. + % Default + % 1.0e-8 + % Accepted + % [0.0; 1.0] + % Example + % param.MSK_DPAR_INTPNT_CO_TOL_PFEAS = 1.0e-8 + % See also + % MSK_DPAR_INTPNT_CO_TOL_NEAR_REL + % Groups + % Interior-point method, Termination criteria, Conic interior-point method + param.MSK_DPAR_INTPNT_CO_TOL_PFEAS=param.feasTol; +end + +if ~isfield(param, 'MSK_DPAR_INTPNT_TOL_DFEAS') + % MSK_DPAR_INTPNT_TOL_DFEAS + % Dual feasibility tolerance used by the interior-point optimizer for linear problems. + % Default + % 1.0e-8 + % Accepted + % [0.0; 1.0] + % Example + % param.MSK_DPAR_INTPNT_TOL_DFEAS = 1.0e-8 + % Groups + % Interior-point method, Termination criteria + param.MSK_DPAR_INTPNT_TOL_DFEAS=param.optTol; +end + +if ~isfield(param, 'MSK_DPAR_INTPNT_QO_TOL_DFEAS') + % Dual feasibility tolerance used by the interior-point optimizer for quadratic problems. + % Default + % 1.0e-8 + % Accepted + % [0.0; 1.0] + % Example + % param.MSK_DPAR_INTPNT_QO_TOL_DFEAS = 1.0e-8 + % See also + % MSK_DPAR_INTPNT_QO_TOL_NEAR_REL + % Groups + % Interior-point method, Termination criteria + param.MSK_DPAR_INTPNT_QO_TOL_DFEAS=param.optTol; +end + +if ~isfield(param, 'MSK_DPAR_INTPNT_CO_TOL_DFEAS') + % Dual feasibility tolerance used by the interior-point optimizer for linear problems. + % Default + % 1.0e-8 + % Accepted + % [0.0; 1.0] + % Example + % param.MSK_DPAR_INTPNT_TOL_DFEAS = 1.0e-8 + % Groups + % Interior-point method, Termination criteria + param.MSK_DPAR_INTPNT_CO_TOL_DFEAS=param.optTol; +end + +if isfield(param,'lifted') && param.lifted==1 + % Controls the maximum amount of fill-in that can be created by one pivot in the elimination phase of the presolve. + % A negative value means the parameter value is selected automatically. + % Default-1 + % Accepted [-inf; +inf] + % Example param.MSK_IPAR_PRESOLVE_ELIMINATOR_MAX_FILL = -1 + if ~isfield(param,'MSK_IPAR_PRESOLVE_ELIMINATOR_MAX_NUM_TRIES') + param.MSK_IPAR_PRESOLVE_ELIMINATOR_MAX_NUM_TRIES = 0; + end +end + +%turn on multiscale if infeasibilities after unscaling +if isfield(param,'multiscale') && param.multiscale==1 && param.lifted==0 + % Controls whether whether a new experimental linear dependency checker is employed. + % Default + % "OFF" + % Accepted + % "ON", "OFF" + % Example + % param.MSK_IPAR_PRESOLVE_LINDEP_NEW = 'MSK_OFF' + if ~isfield(param,'MSK_IPAR_PRESOLVE_LINDEP_NEW') + param.MSK_IPAR_PRESOLVE_LINDEP_NEW = 'MSK_OFF'; + end + + % MSK_IPAR_INTPNT_SCALING + % Controls how the problem is scaled before the interior-point optimizer is used. + % Default + % "FREE" + % Accepted + % "FREE", "NONE" + % param..MSK_IPAR_INTPNT_SCALING = 'MSK_SCALING_FREE'; + if ~isfield(param,'MSK_IPAR_INTPNT_SCALING') + param.MSK_IPAR_INTPNT_SCALING='MSK_SCALING_NONE'; + end + % MSK_IPAR_SIM_SCALING + % Controls how much effort is used in scaling the problem before a simplex optimizer is used. + % Default + % "FREE" + % Accepted + % "FREE", "NONE" + % Example + % param.MSK_IPAR_SIM_SCALING = 'MSK_SCALING_FREE' + if ~isfield(param,'MSK_IPAR_SIM_SCALING') + param.MSK_IPAR_SIM_SCALING='MSK_SCALING_NONE'; + end +end + +if isfield(param,'debug') && param.debug==1 + % https://docs.mosek.com/latest/rmosek/debugging-infeas.html + % Controls whether an infeasibility report is automatically produced after the optimization if the problem is primal or dual infeasible. + param.MSK_IPAR_INFEAS_REPORT_AUTO='MSK_ON'; +end + +if isfield(param,'strict') + % MSK_IPAR_BI_IGNORE_MAX_ITER + % If the parameter MSK_IPAR_INTPNT_BASIS has the value MSK_BI_NO_ERROR and + % the interior-point optimizer has terminated due to maximum number of + % iterations, then basis identification is performed if this parameter has + % the value MSK_ON. + % Possible Values: + % MSK_ON Switch the option on. + % MSK_OFF Switch the option off. + % Default value: + % MSK_OFF + if ~isfield(param,'MSK_IPAR_BI_IGNORE_MAX_ITER') + param.MSK_IPAR_BI_IGNORE_MAX_ITER='MSK_OFF'; + end + + %%%%%%%%%%% Solution Approach + % MSK_IPAR_INTPNT_SOLVE_FORM + % Controls whether the primal or the dual problem is solved. + % Possible Values: + % MSK_SOLVE_PRIMAL + % The optimizer should solve the primal problem. + % MSK_SOLVE_DUAL + % The optimizer should solve the dual problem. + % MSK_SOLVE_FREE + % The optimizer is free to solve either the primal or the dual problem. + % Default value:MSK_SOLVE_FREE + if ~isfield(param,'MSK_IPAR_INTPNT_SOLVE_FORM') + param.MSK_IPAR_INTPNT_SOLVE_FORM='MSK_SOLVE_FREE'; + %param.MSK_IPAR_INTPNT_SOLVE_FORM='MSK_SOLVE_PRIMAL'; + end + + %%%%%%% Infeasibility + % MSK_DPAR_INTPNT_TOL_INFEAS + % Controls when the optimizer declares the model primal or dual infeasible. + % A small number means the optimizer gets more conservative about declaring the model infeasible. + % Possible Values:Any number between 0.0 and 1.0. + % Default value: 1.0e-8 + if ~isfield(param,'MSK_DPAR_INTPNT_TOL_INFEAS') + % param.MSK_DPAR_INTPNT_TOL_INFEAS=1e-10; + param.MSK_DPAR_INTPNT_TOL_INFEAS=1e-8; + end +end + +% param.lpmethod='BARRIER'; +% param.qpmethod='BARRIER'; +% method = param.lpmethod; +% case 'ibm_cplex' +% param.lpmethod='BARRIER'; +% param.qpmethod='BARRIER'; +% method = param.lpmethod; +% case 'mosek' +% method = 'FREE'; +% method = 'INTPNT'; +% method = 'CONIC'; +% + + +%backward compatibility +if isfield(param,'method') + if isempty(param.method) + param = rmfield(param,'method'); + else + if ~isfield(param,[lower(param.problemType) 'method']) + param.([lower(param.problemType) 'method'])=param.method; + end + end +end + +switch param.problemType + case {'LP'} + if isfield(param,'lpmethod') + if contains(param.lpmethod,'MSK_OPTIMIZER_') + param.MSK_IPAR_OPTIMIZER=param.lpmethod; + else + param.MSK_IPAR_OPTIMIZER=['MSK_OPTIMIZER_' param.lpmethod]; + end + end + case {'QP'} + if isfield(param,'qpmethod') + if contains(param.qpmethod,'MSK_OPTIMIZER_') + param.MSK_IPAR_OPTIMIZER=param.qpmethod; + else + param.MSK_IPAR_OPTIMIZER=['MSK_OPTIMIZER_' param.qpmethod]; + end + end + case {'CLP'} + if isfield(param,'clpmethod') + if contains(param.qpmethod,'MSK_OPTIMIZER_') + param.MSK_IPAR_OPTIMIZER=param.clpmethod; + else + param.MSK_IPAR_OPTIMIZER=['MSK_OPTIMIZER_' param.clpmethod]; + end + end + case {'EP'} + if isfield(param,'epmethod') + if contains(param.epmethod,'MSK_OPTIMIZER_') + param.MSK_IPAR_OPTIMIZER=param.epmethod; + else + param.MSK_IPAR_OPTIMIZER=['MSK_OPTIMIZER_' param.epmethod]; + end + end + % MSK_IPAR_INTPNT_REGULARIZATION_USE + % Description:Controls whether regularization is allowed. + % Possible Values: MSK_ON Switch the option on. + % MSK_OFF Switch the option off. + % Default value: MSK_ON + if ~isfield(param,'MSK_IPAR_INTPNT_REGULARIZATION_USE') + param.MSK_IPAR_INTPNT_REGULARIZATION_USE='MSK_ON'; + end + + + + % MSK_IPAR_INTPNT_MAX_ITERATIONS + % Controls the maximum number of iterations allowed in the interior-point optimizer. + % Possible Values:Any number between 0 and +inf. + % Default value: 400 + if ~isfield(param,'MSK_IPAR_INTPNT_MAX_ITERATIONS') + param.MSK_IPAR_INTPNT_MAX_ITERATIONS=400; + end + + %%%%%%%%%%%%% NONLINEAR SOLVER INTEGER PARAM %%%%%%%%%%%%%%%%%%%%%%%%%%%% + + case {'VK'} + + % MSK_IPAR_SIM_SCALING_METHOD + % Controls how the problem is scaled before a simplex optimizer is used. + % Default + % "POW2" + % Accepted + % "POW2", "FREE" + % Example + % param.MSK_IPAR_SIM_SCALING_METHOD = 'MSK_SCALING_METHOD_POW2' + % param.MSK_IPAR_SIM_SCALING_METHOD='MSK_SCALING_METHOD_FREE'; +end + +% Remove outer function specific parameters to avoid crashing solver interfaces +mosekParam = mosekParamStrip(param); +if 0 + disp(mosekParam) +end diff --git a/src/base/solvers/param/getCobraSolverParams.m b/src/base/solvers/param/getCobraSolverParams.m index 96f94931bf..2ea0d2f56b 100644 --- a/src/base/solvers/param/getCobraSolverParams.m +++ b/src/base/solvers/param/getCobraSolverParams.m @@ -1,6 +1,6 @@ function varargout = getCobraSolverParams(problemType, paramNames, paramStructure) % This function gets the specified paramStructure in `paramNames` from -% paramStructure, the global cobra paramters variable or default values set within +% paramStructure, the global cobra parameters variable or default values set within % this script. % % It will use values with the following priority @@ -43,21 +43,47 @@ % - Richard Que (12/01/2009) % - Ronan (16/07/2013) default MPS paramStructure are no longer global variables -if nargin < 2 - error('getCobraSolverParams: No paramStructure specified') -end - -if nargin < 3 - paramStructure = []; +if ~exist('paramNames','var') || isempty(paramNames) + %get the names of all the parameters + defaultParamNames = getCobraSolverParamsOptionsForType(problemType); +else + if ~iscell(paramNames) + paramNames = {paramNames}; + end end -if ~isempty(paramStructure) - paramStructureOut = paramStructure; +if exist('paramStructure','var') + overRideDefaults = 1; + %overide defaults unless hack is asking for defaults + if isstruct(paramStructure) + paramStructureOut = paramStructure; + else + if ischar(paramStructure) && strcmp(paramStructure,'default') + overRideDefaults=1; + end + end +else + %the last argument will be the structure with defaults + overRideDefaults=0; + paramStructure=[]; end + % Persistence will make names specific to one type of solver. % Default Values % For descriptions of the different settings please have a look at % getCobraSolverParamsOptionsForType -valDef.minNorm = 0; + +% %These default tolerances are based on the default values for the Gurobi LP +% %solver. Do not change them without first consulting with other developers. +% %https://www.gurobi.com/documentation/9.0/refman/parameters.html +% % (primal) feasibility tolerance +% changeCobraSolverParams('LP', 'feasTol', 1e-6); +% % (dual) optimality tolerance +% changeCobraSolverParams('LP', 'optTol', 1e-6); +% +% % (primal) feasibility tolerance +% changeCobraSolverParams('EP', 'feasTol', 1e-8); +% % (dual) optimality tolerance +% changeCobraSolverParams('EP', 'optTol', 1e-12); %These default tolerances are based on the default values for the Gurobi LP solver %https://www.gurobi.com/documentation/9.0/refman/paramStructure.html @@ -76,20 +102,12 @@ valDef.timeLimit = 1e36; valDef.iterationLimit = 1000; -%valDef.logFile = ['Cobra' problemType 'Solver.log']; valDef.logFile = []; %log file should be empty to avoid creating it by default valDef.saveInput = []; -valDef.PbName = [problemType 'problem']; +valDef.problemType = problemType; valDef.debug = 0; valDef.lifting = 0; - -valDef.method = -1; - -% CPLEX paramStructure -valDef.DATACHECK = 1; -valDef.DEPIND = 1; -valDef.checkNaN = 0; -valDef.warning = 0; +valDef.multiscale=0; % true if problem is multiscale % tolerances valDef.intTol = 1e-12; @@ -97,10 +115,6 @@ valDef.absMipGapTol = 1e-12; valDef.NUMERICALEMPHASIS = 1; -if ~iscell(paramNames) - paramNames = {paramNames}; -end - switch problemType case 'LP' global CBT_LP_PARAMS @@ -117,9 +131,13 @@ valDef.feasTol = 1e-6; % (primal) feasibility tolerance valDef.optTol = 1e-6; % (dual) optimality tolerance valDef.solver='mosek'; + valDef.epmethod='CONIC'; case 'CLP' + % This is never used elsewhere except for parameter setting loop + % for backward compatibility global CBT_CLP_PARAMS parametersGlobal = CBT_CLP_PARAMS; + valDef.feasTol = 1e-6; % (primal) feasibility tolerance valDef.optTol = 1e-6; % (dual) optimality tolerance valDef.solver='mosek'; @@ -134,27 +152,44 @@ return; end +if exist('defaultParamNames','var') + % set parameter structure values to default, unless they are overridden by paramStructure + for i=1:length(defaultParamNames) + if isfield(valDef,defaultParamNames{i}) + %structure of default parameter values + paramStructureOut.(defaultParamNames{i})=valDef.(defaultParamNames{i}); + end + end +end + varargout = cell(1, numel(paramNames)); paramNames = columnVector(paramNames); for i=1:length(paramNames) - % set values to default + % first set each value to default if isfield(valDef,paramNames{i}) + %list of default parameter values varargout{i} = valDef.(paramNames{i}); + %structure of default parameter values + paramStructureOut.(paramNames{i})= valDef.(paramNames{i}); end - if ~strcmp(paramStructure,'default') % skip of using default values - % set values to global values + + if overRideDefaults + % second set values to global values if isfield(parametersGlobal,paramNames{i}) - varargout{i} = parametersGlobal.(paramNames{i}); + %list of default parameter values + varargout{i} = parametersGlobal.(paramNames{i}); + %structure of default parameter values + paramStructureOut.(paramNames{i})=parametersGlobal.(paramNames{i}); end - % set values to specified values + % third set values to specified values in paramStructure if isfield(paramStructure,paramNames{i}) - varargout{i} = paramStructure.(paramNames{i}); - %remove this field so paramStructureOut can be passed directly to the solver, without extra fields - paramStructureOut = rmfield(paramStructureOut,paramNames{i}); + varargout{i} = paramStructure.(paramNames{i}); + %structure of default parameter values + paramStructureOut.(paramNames{i})=paramStructure.(paramNames{i}); end end end -if ~isempty(paramStructure) +if overRideDefaults varargout{length(paramNames)+1}=paramStructureOut; end diff --git a/src/base/solvers/param/getCobraSolverParamsOptionsForType.m b/src/base/solvers/param/getCobraSolverParamsOptionsForType.m index d228910f21..cee1c4cc0c 100644 --- a/src/base/solvers/param/getCobraSolverParamsOptionsForType.m +++ b/src/base/solvers/param/getCobraSolverParamsOptionsForType.m @@ -7,24 +7,26 @@ % % INPUT: % problemType : One of the problem types available in the COBRA -% Toolbox ('LP','QP','MILP','MIQP','NLP') +% Toolbox ('LP','QP','EP','CLP','MILP','MIQP','NLP') % % OUPTUT: % paramNames: Cell array of names of parameters that can be set % for each problem type, independent of the specific % solver being used. -if iscell(problemType ) +if iscell(problemType) paramNames = {}; - for j = 1:numel(problemType ) - paramNames = [paramNames, getCobraSolverParamsOptionsForType(problemType {j})]; + for j = 1:numel(problemType) + paramNames = [paramNames, getCobraSolverParamsOptionsForType(problemType{j})]; end paramNames = unique(paramNames); return end switch problemType case 'LP' - paramNames = {'verify',... % verify that it is a suitable LP problem + paramNames = {'multiscale'... % true if problem is multiscale + 'problemType'... % problem type + 'verify',... % verify that it is a suitable LP problem 'minNorm', ... % type of normalization used. 'printLevel', ... % print Level 'primalOnly', ... % only solve for primal @@ -37,7 +39,9 @@ 'lifting', ... % whether to lift a problem 'method'}; % solver method: -1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent, 5 = Network Solver(if supported by the solver) case 'QP' - paramNames = {'verify',... % verify that it is a suitable QP problem + paramNames = {'multiscale'... % true if problem is multiscale + 'problemType'... % problem type + 'verify',... % verify that it is a suitable QP problem 'method', ... % solver method: -1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent, 5 = Network Solver(if supported by the solver) 'printLevel', ... % print level 'saveInput', ... % save the input to a file (specified) @@ -48,7 +52,9 @@ 'solver'}; % the solver to use case 'EP' - paramNames = {'verify',... % verify that it is a suitable EP problem + paramNames = {'multiscale'... % true if problem is multiscale + 'problemType'... % problem type + 'verify',... % verify that it is a suitable EP problem 'method', ... % solver method: -1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent, 5 = Network Solver(if supported by the solver) 'printLevel', ... % print level 'debug', ... % run debgugging code @@ -56,7 +62,9 @@ 'optTol',... % optimality tolerance 'solver'}; % the solver to use case 'CLP' - paramNames = {'verify',... % verify that it is a suitable CLP problem + paramNames = {'multiscale'... % true if problem is multiscale + 'problemType'... % problem type + 'verify',... % verify that it is a suitable CLP problem 'printLevel', ... % print level 'debug', ... % run debgugging code 'feasTol',... % feasibility tolerance @@ -65,7 +73,8 @@ }; case 'MILP' - paramNames = {'intTol', ... % integer tolerance (accepted derivation from integer numbers) + paramNames = {'problemType'... % problem type + 'intTol', ... % integer tolerance (accepted derivation from integer numbers) 'relMipGapTol', ... % relative MIP Gap tolerance 'absMipGapTol', ... % absolute MIP Gap tolerance 'timeLimit', ... % maximum time before stopping computation (if supported by the solver) @@ -78,7 +87,8 @@ 'debug'}; % run debgugging code case 'MIQP' - paramNames = {'timeLimit', ... % maximum time before stopping computation (if supported by the solver) + paramNames = {'problemType'... % problem type + 'timeLimit', ... % maximum time before stopping computation (if supported by the solver) 'method', ... % solver method: -1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent, 5 = Network Solver(if supported by the solver) 'feasTol',... % feasibility tolerance 'optTol',... % optimality tolerance @@ -91,7 +101,8 @@ 'solver'}; % the solver to use case 'NLP' - paramNames = {'warning', ... % whether to display warnings + paramNames = {'problemType'... % problem type + 'warning', ... % whether to display warnings 'checkNaN', ... % check for NaN solutions 'PbName', ... % name of the problem 'iterationLimit', ... % maximum number of iterations before stopping computation (if supported by the solver) diff --git a/src/base/solvers/param/mergeCobraParams.m b/src/base/solvers/param/mergeCobraParams.m new file mode 100644 index 0000000000..6e3e991a57 --- /dev/null +++ b/src/base/solvers/param/mergeCobraParams.m @@ -0,0 +1,11 @@ +function param = mergeCobraParams(param1,param2) +% Merge the structures with preference to param1 fields over param2 fields + +param = param2; % Start with the first structure +fields = fieldnames(param1); % Get the field names of the second structure +for i = 1:numel(fields) + param.(fields{i}) = param1.(fields{i}); % struct2 values will overwrite struct1 values if overlap occurs +end + +end + diff --git a/src/base/solvers/param/parseSolverParameters.m b/src/base/solvers/param/parseSolverParameters.m index 8198ec37ca..2eb8eae37f 100644 --- a/src/base/solvers/param/parseSolverParameters.m +++ b/src/base/solvers/param/parseSolverParameters.m @@ -103,35 +103,34 @@ end end -%move following set of parameters from solverOnlyParams to param -if isfield(solverOnlyParams,'maxConc') - param.maxConc = solverOnlyParams.maxConc; - solverOnlyParams = rmfield(solverOnlyParams,'maxConc'); -end -if isfield(solverOnlyParams,'method') - param.method = solverOnlyParams.method; - solverOnlyParams = rmfield(solverOnlyParams,'method'); -end -if isfield(solverOnlyParams,'maxUnidirectionalFlux') - param.maxUnidirectionalFlux = solverOnlyParams.maxUnidirectionalFlux; - solverOnlyParams = rmfield(solverOnlyParams,'maxUnidirectionalFlux'); -end -if isfield(solverOnlyParams,'minUnidirectionalFlux') - param.minUnidirectionalFlux = solverOnlyParams.minUnidirectionalFlux; - solverOnlyParams = rmfield(solverOnlyParams,'minUnidirectionalFlux'); -end -if isfield(solverOnlyParams,'internalNetFluxBounds') - param.internalNetFluxBounds = solverOnlyParams.internalNetFluxBounds; - solverOnlyParams = rmfield(solverOnlyParams,'internalNetFluxBounds'); -end -if isfield(solverOnlyParams,'externalNetFluxBounds') - param.externalNetFluxBounds = solverOnlyParams.externalNetFluxBounds; - solverOnlyParams = rmfield(solverOnlyParams,'externalNetFluxBounds'); -end -if isfield(solverOnlyParams,'rounding') - param.rounding = solverOnlyParams.rounding; - solverOnlyParams = rmfield(solverOnlyParams,'rounding'); -end -if isfield(param,'printLevel') - solverOnlyParams.printLevel = param.printLevel -1; -end +if 0 + %move following set of parameters from solverOnlyParams to param + if isfield(solverOnlyParams,'maxConc') + param.maxConc = solverOnlyParams.maxConc; + solverOnlyParams = rmfield(solverOnlyParams,'maxConc'); + end + if isfield(solverOnlyParams,'method') + param.method = solverOnlyParams.method; + solverOnlyParams = rmfield(solverOnlyParams,'method'); + end + if isfield(solverOnlyParams,'maxUnidirectionalFlux') + param.maxUnidirectionalFlux = solverOnlyParams.maxUnidirectionalFlux; + solverOnlyParams = rmfield(solverOnlyParams,'maxUnidirectionalFlux'); + end + if isfield(solverOnlyParams,'minUnidirectionalFlux') + param.minUnidirectionalFlux = solverOnlyParams.minUnidirectionalFlux; + solverOnlyParams = rmfield(solverOnlyParams,'minUnidirectionalFlux'); + end + if isfield(solverOnlyParams,'internalNetFluxBounds') + param.internalNetFluxBounds = solverOnlyParams.internalNetFluxBounds; + solverOnlyParams = rmfield(solverOnlyParams,'internalNetFluxBounds'); + end + if isfield(solverOnlyParams,'externalNetFluxBounds') + param.externalNetFluxBounds = solverOnlyParams.externalNetFluxBounds; + solverOnlyParams = rmfield(solverOnlyParams,'externalNetFluxBounds'); + end + if isfield(solverOnlyParams,'rounding') + param.rounding = solverOnlyParams.rounding; + solverOnlyParams = rmfield(solverOnlyParams,'rounding'); + end +end \ No newline at end of file diff --git a/src/base/solvers/solveCobraLP.m b/src/base/solvers/solveCobraLP.m index 7da85a4155..e0eb9c9d5c 100644 --- a/src/base/solvers/solveCobraLP.m +++ b/src/base/solvers/solveCobraLP.m @@ -61,7 +61,7 @@ % * .rcost: Reduced costs, dual solution to :math:`lb <= v <= ub` % * .dual: dual solution to `A*v ('E' | 'G' | 'L') b` % * .solver: Solver used to solve LP problem -% * .method: Algorithm used by solver to solve LP problem +% * .lpmethod: Algorithm used by solver to solve LP problem % * .stat: Solver status in standardized form % % * 0 - Infeasible problem @@ -222,7 +222,7 @@ stat = 0; origStat = []; origStatText = []; -method = ''; +lpmethod = ''; t_start = clock; if isempty(solver) @@ -545,6 +545,12 @@ case 'glpk' %% GLPK + % msglev (default: 1) + % Level of messages output by solver routines: + % 0 - No output. + % 1 - Error messages only. + % 2 - Normal output. + % 3 - Full output (includes informational messages). param.msglev = problemTypeParams.printLevel; % level of verbosity param.tolbnd = problemTypeParams.feasTol; % tolerance param.toldj = problemTypeParams.optTol; % tolerance @@ -630,55 +636,9 @@ w = []; case 'mosek' % mosek - % use msklpopt with full control over all mosek parameters - % http://docs.mosek.com/7.0/toolbox/Parameters.html - % see also - % http://docs.mosek.com/7.0/toolbox/A_guided_tour.html#SEC:VIEWSETPARAM - % e.g. - % http://docs.mosek.com/7.0/toolbox/MSK_IPAR_OPTIMIZER.html - - %[rcode,res] = mosekopt('param echo(0)',[],solverParams); - - % Remove outer function specific parameters to avoid crashing solver interfaces - solverParams = mosekParamStrip(solverParams); - - param = solverParams; - % only set the print level if not already set via solverParams structure - if ~isfield(solverParams, 'MSK_IPAR_LOG') - switch problemTypeParams.printLevel - case 0 - echolev = 0; - case 1 - echolev = 3; - case 2 - solverParams.MSK_IPAR_LOG_INTPNT = 1; - solverParams.MSK_IPAR_LOG_SIM = 1; - echolev = 3; - otherwise - echolev = 0; - end - if echolev == 0 - solverParams.MSK_IPAR_LOG = 0; - cmd = ['minimize echo(' int2str(echolev) ')']; - else - cmd = 'minimize'; - end - end - - %https://docs.mosek.com/8.1/toolbox/solving-linear.html - if ~isfield(solverParams, 'MSK_DPAR_INTPNT_TOL_PFEAS') - solverParams.MSK_DPAR_INTPNT_TOL_PFEAS=problemTypeParams.feasTol; - end - if ~isfield(solverParams, 'MSK_DPAR_INTPNT_TOL_DFEAS.') - solverParams.MSK_DPAR_INTPNT_TOL_DFEAS=problemTypeParams.feasTol; - end - %If the feasibility tolerance is changed by the solverParams - %struct, this needs to be forwarded to the cobra Params for the - %final consistency test! - if isfield(solverParams,'MSK_DPAR_INTPNT_TOL_PFEAS') - problemTypeParams.feasTol = solverParams.MSK_DPAR_INTPNT_TOL_PFEAS; - end - + param = mergeCobraParams(solverParams,problemTypeParams); + [cmd,mosekParam] = setMosekParam(param); + % basis reuse - TODO % http://docs.mosek.com/7.0/toolbox/A_guided_tour.html#section-node-_A%20guided%20tour_Advanced%20start%20%28hot-start%29 @@ -708,13 +668,12 @@ % blx and bux. Note -inf is allowed in blc and blx. % Similarly, inf is allowed in buc and bux. - + blc = b; buc = b; buc(csense == 'G') = inf; blc(csense == 'L') = -inf; - prob.c = osense * c; prob.a = A; prob.blc = blc; @@ -729,8 +688,18 @@ prob.sol.bas.xx = basis.xx; end - [rcode,res] = mosekopt(cmd,prob,solverParams); + if param.debug + probBeforeMosekopt = prob; + save('probBeforeMosekopt','probBeforeMosekopt'); + end + + [rcode,res] = mosekopt(cmd,prob,mosekParam); + if isfield(param,'lpmethod') + lpmethod = param.lpmethod; + else + lpmethod = 'FREE'; + end if rcode~=0 % MSK_RES_TRM_STALL @@ -745,7 +714,7 @@ end %parse mosek result structure - [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res,solverParams,problemTypeParams.printLevel); + [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res); if stat ==1 || stat ==3 f = c'*x; @@ -762,10 +731,7 @@ f = NaN; s = NaN*ones(size(A,1),1); end - - if isfield(param,'MSK_IPAR_OPTIMIZER') - method=param.MSK_IPAR_OPTIMIZER; - end + case 'mosek_linprog' %% mosek % if mosek is installed, and the paths are added ahead of matlab's @@ -822,88 +788,8 @@ case 'gurobi' % Free academic licenses for the Gurobi solver can be obtained from % http://www.gurobi.com/html/academic.html - % resultgurobi = struct('x',[],'objval',[],'pi',[]); - - % The params struct contains Gurobi parameters. A full list may be - % found on the Parameter page of the reference manual: - % https://www.gurobi.com/documentation/current/refman/parameter_descriptions.html - % MATLAB Parameter Examples - % In the MATLAB interface, parameters are passed to Gurobi through a struct. - % To modify a parameter, you create a field in the struct with the appropriate name, - % and set it to the desired value. For example, to set the TimeLimit parameter to 100 you'd do: - % - % params.timelimit = 100; - % The case of the parameter name is ignored, as are underscores. Thus, you could also do: - % params.timeLimit = 100; - % ...or... - % params.TIME_LIMIT = 100; - % All desired parameter changes should be stored in a single struct, which is passed as the second parameter to the gurobi function. - param=solverParams; - - % https://www.gurobi.com/documentation/current/refman/method.html - % params.method gives the method used to solve continuous models - % -1=automatic, - % 0=primal simplex, - % 1=dual simplex, - % 2=barrier, - % 3=concurrent, - % 4=deterministic concurrent - % i.e. params.method = 1; % use dual simplex method - if isfield(param,'lpmethod') - %gurobiAlgorithms = {'AUTOMATIC','PRIMAL','DUAL','BARRIER','CONCURRENT','CONCURRENT_DETERMINISTIC'}; - % -1=automatic, - % 0=primal simplex, - % 1=dual simplex, - % 2=barrier, - % 3=concurrent, - % 4=deterministic concurrent - switch param.lpmethod - case 'AUTOMATIC' - param.method = -1; - case 'PRIMAL' - param.method = 0; - case 'DUAL' - param.method = 2; - case 'BARRIER' - param.method = 2; - case 'CONCURRENT' - param.method = 3; - case 'DETERMINISTIC_CONCURRENT' - param.method = 4; - otherwise - error('Unrecognised param.lpmethod for gurobi') - end - param = rmfield(param,'lpmethod'); - end - - if ~isfield(param,'OutputFlag') - switch problemTypeParams.printLevel - case 0 - param.OutputFlag = 0; - param.DisplayInterval = 1; - case 1 - param.OutputFlag = 0; - param.DisplayInterval = 1; - otherwise - % silent - param.OutputFlag = 0; - param.DisplayInterval = 1; - end - end - - if isfield(param,'FeasibilityTol') - % update tolerance according to actual setting - problemTypeParams.feasTol = param.FeasibilityTol; - else - param.FeasibilityTol = problemTypeParams.feasTol; - end - - if isfield(param,'OptimalityTol') - % update tolerance according to actual setting - problemTypeParams.optTol = param.OptimalityTol; - else - param.OptimalityTol = problemTypeParams.optTol; - end + param = mergeCobraParams(solverParams,problemTypeParams); + gurobiParam = setGurobiParam(param); gurobiLP.sense(1:length(b),1) = '='; gurobiLP.sense(csense == 'L') = '<'; @@ -930,24 +816,14 @@ %gurobi wants a dense double vector as an objective gurobiLP.obj = double(c)+0;%full - % basis reuse - Ronan if ~isempty(basis) gurobiLP.cbasis = full(basis.cbasis); gurobiLP.vbasis = full(basis.vbasis); end - % set the solver specific parameters - param = updateStructData(param,solverParams); -% LPproblem = rmfield(LPproblem,'c'); -% LPproblem = rmfield(LPproblem,'b'); -% LPproblem = rmfield(LPproblem,'lb'); -% LPproblem = rmfield(LPproblem,'ub'); -% LPproblem = rmfield(LPproblem,'osense'); -% LPproblem = rmfield(LPproblem,'csense'); - % call the solver - resultgurobi = gurobi(gurobiLP,param); + resultgurobi = gurobi(gurobiLP,gurobiParam); % see the solvers original status -Ronan origStat = resultgurobi.status; @@ -1204,7 +1080,7 @@ [solution,LPprob] = solveCobraLPCPLEX(LPproblem,problemTypeParams.printLevel,1,[],[],minNorm); solution.basis = LPprob.LPBasis; solution.solver = solver; - solution.method = method; % dummy + solution.lpmethod = lpmethod; % dummy % solution.slack = []; if exist([pwd filesep 'clone1.log'],'file') delete('clone1.log') @@ -1317,7 +1193,7 @@ %this is the dual to the simple ineequality constraints : reduced costs w = lambda.lower - lambda.upper; - method = output.method; + lpmethod = output.method; if 0 %debug disp(method) norm(osense * c - A' * y - w,inf) @@ -1364,11 +1240,14 @@ x = CplexLPproblem.Solution.x; if 1 f = c'*x; + w = CplexLPproblem.Solution.reducedcost; + y = CplexLPproblem.Solution.dual; else f = CplexLPproblem.Solution.objval; %do not use as it gives the opposite sign for maximise + w = osense*CplexLPproblem.Solution.reducedcost; + y = osense*CplexLPproblem.Solution.dual; end - w = osense*CplexLPproblem.Solution.reducedcost; - y = osense*CplexLPproblem.Solution.dual; + %res1 = A*solution.full + solution.slack - b; s = b - A * x; % output the slack variables @@ -1417,6 +1296,8 @@ if isfield(CplexLPproblem.Solution ,'dual') y = osense*CplexLPproblem.Solution.dual; end + elseif origStat==101 + warning('101') else stat = -1; end @@ -1645,11 +1526,11 @@ if ~strcmp(solver,'cplex_direct') && ~strcmp(solver,'mps') % assign solution - t = etime(clock, t_start); + time = etime(clock, t_start); if ~exist('basis','var'), basis=[]; end [solution.full, solution.obj, solution.rcost, solution.dual, solution.slack, ... - solution.solver, solution.method, solution.stat, solution.origStat, ... - solution.origStatText,solution.time,solution.basis] = deal(x,f,w,y,s,solver,method,stat,origStat,origStatText,t,basis); + solution.solver, solution.lpmethod, solution.stat, solution.origStat, ... + solution.origStatText,solution.time,solution.basis] = deal(x,f,w,y,s,solver,lpmethod,stat,origStat,origStatText,time,basis); elseif strcmp(solver,'mps') solution = []; end diff --git a/src/base/solvers/solveCobraQP.m b/src/base/solvers/solveCobraQP.m index 07c959aba5..2670b7408f 100644 --- a/src/base/solvers/solveCobraQP.m +++ b/src/base/solvers/solveCobraQP.m @@ -62,6 +62,7 @@ % * .slack: slack variable such that :math:`A*x + s = b` % * .obj: Objective value % * .solver: Solver used to solve QP problem +% * .qpmethod: Solver qpmethod used to solve QP problem % * .origStat: Original status returned by the specific solver % * .time: Solve time in seconds % * .stat: Solver status in standardized form (see below) @@ -323,8 +324,9 @@ stat = 3; % Solution exists, but either scaling problems or not proven to be optimal else %(origStat >= 10) stat = -1; % No optimal solution found (time or other limits reached, other infeasibility problems) - end - + end + origStat = cplexProblem.Solution.statusstring; + %Update Tolerance According to actual setting problemTypeParams.feasTol = cplexProblem.Param.simplex.tolerances.feasibility.Cur; problemTypeParams.optTol = cplexProblem.Param.simplex.tolerances.optimality.Cur; @@ -372,12 +374,10 @@ %% case 'mosek' - if problemTypeParams.printLevel>0 - cmd='minimize'; - else - cmd='minimize echo(0)'; - end + param = mergeCobraParams(solverParams,problemTypeParams); + [cmd,mosekParam] = setMosekParam(param); + %matching bounds and zero diagonal of F at the same time bool = lb == ub & diag(F)==0; if any(bool) @@ -389,38 +389,6 @@ %} warning(['There are ' num2str(nnz(bool)) ' variables that have equal lower and upper bounds, and zero on the diagonal of F.']) end - - param = struct(); - % Set the printLevel, can be overwritten. - if ~isfield(param, 'MSK_IPAR_LOG') - switch problemTypeParams.printLevel - case 0 - echolev = 0; - case 1 - echolev = 3; - case 2 - param.MSK_IPAR_LOG_INTPNT = 1; - param.MSK_IPAR_LOG_SIM = 1; - echolev = 3; - otherwise - echolev = 0; - end - if echolev == 0 - param.MSK_IPAR_LOG = 0; - cmd = ['minimize echo(' int2str(echolev) ')']; - else - cmd = 'minimize'; - end - end - %remove parameter fields that mosek does not recognise - param.MSK_DPAR_INTPNT_QO_TOL_DFEAS = problemTypeParams.optTol; - param.MSK_DPAR_INTPNT_QO_TOL_PFEAS = problemTypeParams.feasTol; - - %Update with solver Specific Parameter struct - param = updateStructData(param,solverParams); - %problemTypeParams.feasTol = param.MSK_DPAR_INTPNT_NL_TOL_PFEAS; - - param = mosekParamStrip(param); blc = b; buc = b; @@ -438,8 +406,15 @@ %https://docs.mosek.com/latest/toolbox/data-types.html#prob [prob.qosubi,prob.qosubj,prob.qoval]=find(F); - [rcode,res] = mosekopt(cmd,prob,param); - + + [rcode,res] = mosekopt(cmd,prob,mosekParam); + + if isfield(param,'qpmethod') + qpmethod = param.qpmethod; + else + qpmethod = 'FREE'; + end + % stat Solver status % 1 Optimal solution found % 2 Unbounded solution @@ -460,7 +435,7 @@ end %parse mosek result structure - [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res,solverParams,problemTypeParams.printLevel); + [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res); %debugging if problemTypeParams.printLevel>2 @@ -493,8 +468,7 @@ f = NaN; s = NaN*ones(size(A,1),1); end - - + case 'pdco' %----------------------------------------------------------------------- % pdco.m: Primal-Dual Barrier Method for Convex Objectives (16 Dec 2008) @@ -633,74 +607,8 @@ % Free academic licenses for the Gurobi solver can be obtained from % http://www.gurobi.com/html/academic.html - % The param struct contains Gurobi parameters. A full list may be - % found on the Parameter page of the reference manual: - % https://www.gurobi.com/documentation/current/refman/parameter_descriptions.html - % MATLAB Parameter Examples - % In the MATLAB interface, parameters are passed to Gurobi through a struct. - % To modify a parameter, you create a field in the struct with the appropriate name, - % and set it to the desired value. For example, to set the TimeLimit parameter to 100 you'd do: - % - % param.timelimit = 100; - % The case of the parameter name is ignored, as are underscores. Thus, you could also do: - % param.timeLimit = 100; - % ...or... - % param.TIME_LIMIT = 100; - % All desired parameter changes should be stored in a single struct, which is passed as the second parameter to the gurobi function. - param=solverParams; - - % param.method gives the method used to solve continuous models - % -1=automatic, - % 0=primal simplex, - % 1=dual simplex, - % 2=barrier, - % 3=concurrent, - % 4=deterministic concurrent - % i.e. param.method = 1; % use dual simplex method - if isfield(param,'lpmethod') - %gurobiAlgorithms = {'AUTOMATIC','PRIMAL','DUAL','BARRIER','CONCURRENT','CONCURRENT_DETERMINISTIC'}; - % -1=automatic, - % 0=primal simplex, - % 1=dual simplex, - % 2=barrier, - % 3=concurrent, - % 4=deterministic concurrent - switch param.lpmethod - case 'AUTOMATIC' - param.method = -1; - case 'PRIMAL' - param.method = 0; - case 'DUAL' - param.method = 2; - case 'BARRIER' - param.method = 2; - otherwise - %https://www.gurobi.com/documentation/current/refman/method.html - %Concurrent methods aren't available for QP and QCP. - error('Unrecognised param.lpmethod for gurobi') - end - param = rmfield(param,'lpmethod'); - end - - switch problemTypeParams.printLevel - case 0 - param.OutputFlag = 0; - param.DisplayInterval = 1; - case problemTypeParams.printLevel>1 - param.OutputFlag = 1; - param.DisplayInterval = 5; - otherwise - param.OutputFlag = 0; - param.DisplayInterval = 1; - end - - param.FeasibilityTol = problemTypeParams.feasTol; - param.OptimalityTol = problemTypeParams.optTol; - %Update param struct with Solver Specific parameters - param = updateStructData(param,solverParams); - - %Update feasTol in case it is changed by the solver Parameters - problemTypeParams.feasTol = param.FeasibilityTol; + param = mergeCobraParams(solverParams,problemTypeParams); + gurobiParam = setGurobiParam(param); gurobiQP.sense(1:length(b),1) = '='; gurobiQP.sense(csense == 'L') = '<'; @@ -730,7 +638,7 @@ end try - resultgurobi = gurobi(gurobiQP,param); + resultgurobi = gurobi(gurobiQP,gurobiParam); catch ME if contains(ME.message,'Gurobi error 10020: Objective Q not PSD (negative diagonal entry)') warning('%s\n','Gurobi cannot solve a QP problem if it is given a diagonal Q with some of those diagonals equal to zero') @@ -748,7 +656,7 @@ [x,f,y,w,s] = deal(resultgurobi.x,resultgurobi.objval,osense*resultgurobi.pi,osense*resultgurobi.rc,resultgurobi.slack); - if problemTypeParams.printLevel>2 %|| 1 + if param.printLevel>2 %|| 1 res1 = A*x + s - b; disp('Check A*x + s - b = 0 (feasiblity):'); disp(norm(res1,inf)) @@ -784,7 +692,7 @@ % if the status becomes 'OPTIMAL', it is unbounded, otherwise it is infeasible. gurobiQP.obj(:) = 0; gurobiQP.F(:,:) = 0; - resultgurobi = gurobi(gurobiQP,param); + resultgurobi = gurobi(gurobiQP,gurobiParam); if strcmp(resultgurobi.status,'OPTIMAL') stat = 2; else @@ -1047,6 +955,11 @@ solution.slack = s; solution.dual = y; solution.rcost = w; +if exist('qpmethod','var') + solution.qpmethod = qpmethod; +else + solution.qpmethod = ''; +end if any(contains(solver,'cplex')) [ExitText,~] = cplexStatus(solution.origStat); solution.origStatText = ExitText; diff --git a/test/verifiedTests/base/testSolvers/createToyModelForLifting.m b/test/verifiedTests/base/testLifting/createToyModelForLifting.m similarity index 100% rename from test/verifiedTests/base/testSolvers/createToyModelForLifting.m rename to test/verifiedTests/base/testLifting/createToyModelForLifting.m diff --git a/test/verifiedTests/base/testSolvers/testLiftModel.m b/test/verifiedTests/base/testLifting/testLiftModel.m similarity index 100% rename from test/verifiedTests/base/testSolvers/testLiftModel.m rename to test/verifiedTests/base/testLifting/testLiftModel.m diff --git a/test/verifiedTests/base/testSolvers/runLPvariousSolvers.m b/test/verifiedTests/base/testSolvers/runLPvariousSolvers.m index 09f6493a77..30d7385c61 100644 --- a/test/verifiedTests/base/testSolvers/runLPvariousSolvers.m +++ b/test/verifiedTests/base/testSolvers/runLPvariousSolvers.m @@ -298,10 +298,10 @@ else tmp_rcost = full(solution{i}.rcost(randrcost)); end - fprintf('%3d%15f%15f%15f%15f%20s\t%30s\n', i, solution{i}.time, solution{i}.obj, tmp_dual, tmp_rcost, solution{i}.solver, solution{i}.method) + fprintf('%3d%15f%15f%15f%15f%20s\t%30s\n', i, solution{i}.time, solution{i}.obj, tmp_dual, tmp_rcost, solution{i}.solver, solution{i}.lpmethod) all_obj(i) = solution{i}.obj; else - fprintf('%3d%15f%15f%15f%15f%20s\t%30s\n', i, solution{i}.time, solution{i}.obj, NaN, NaN, solution{i}.solver, solution{i}.method) + fprintf('%3d%15f%15f%15f%15f%20s\t%30s\n', i, solution{i}.time, solution{i}.obj, NaN, NaN, solution{i}.solver, solution{i}.lpmethod) all_obj(i) = 0.0; end end diff --git a/test/verifiedTests/base/testSolvers/testSolveCobraLP.m b/test/verifiedTests/base/testSolvers/testSolveCobraLP.m index 2acfa54d09..e7d16cb9ee 100644 --- a/test/verifiedTests/base/testSolvers/testSolveCobraLP.m +++ b/test/verifiedTests/base/testSolvers/testSolveCobraLP.m @@ -168,8 +168,8 @@ osenseStr = 'max'; minNorm = 'zero'; allowLoops = 1; -param.zeroNormApprox = 'all'; -solution = optimizeCbModel(model, osenseStr,minNorm, allowLoops, param); +optimizeCbModel_param.zeroNormApprox = 'all'; +solution = optimizeCbModel(model, osenseStr,minNorm, allowLoops, optimizeCbModel_param); assert(solution.f0==1) %% % change the directory From 39b65b94577575046f17c8739549a5370bd4a3ba Mon Sep 17 00:00:00 2001 From: Ronan Fleming Date: Tue, 5 Nov 2024 20:16:36 +0000 Subject: [PATCH 2/3] updates to solvers --- src/analysis/FBA/optimizeCbModel.m | 21 +++++----- src/base/solvers/buildOptProblemFromModel.m | 14 +++++-- .../cplex/setCplexParametersForProblem.m | 19 +++++++++ .../entropicFBA/entropicFluxBalanceAnalysis.m | 9 +++-- .../entropicFBA/processConcConstraints.m | 23 ++++------- .../entropicFBA/processFluxConstraints.m | 7 ++-- src/base/solvers/entropicFBA/solveCobraEP.m | 36 +++++++++++------ .../solvers/getSetSolver/changeCobraSolver.m | 4 +- src/base/solvers/gurobi/setGurobiParam.m | 14 +++++++ src/base/solvers/mosek/parseMskResult.m | 21 +++++----- src/base/solvers/mosek/setMosekParam.m | 39 ++++++++++++------- src/base/solvers/param/getCobraSolverParams.m | 2 +- .../getCobraSolverParamsOptionsForType.m | 5 +-- .../solvers/param/parseSolverParameters.m | 32 --------------- src/base/solvers/solveCobraLP.m | 27 +++++++++---- src/base/solvers/solveCobraQP.m | 9 ++++- 16 files changed, 163 insertions(+), 119 deletions(-) diff --git a/src/analysis/FBA/optimizeCbModel.m b/src/analysis/FBA/optimizeCbModel.m index 2e5be448d2..72a9b8bdcf 100644 --- a/src/analysis/FBA/optimizeCbModel.m +++ b/src/analysis/FBA/optimizeCbModel.m @@ -418,10 +418,10 @@ % If this is a quadratically regularised LP, go straight to QP % TODO This is a hack of the param.minNorm to direct solution to QRLP or QRQP if isfield(param,'solveWBMmethod') - if any(strcmp(param.solveWBMmethod,{'QRLP','QRQP'})) + if any(strcmp(param.solveWBMmethod,{'QP','QRLP','QRQP'}))%TODO 'QRLP','QRQP' need coded in + model.c(:)=0; doLinearOptimisationFirst = 0; - param.minNormWBM = param.minNorm; - param.minNorm = param.solveWBMmethod; + minNorm = param.minNorm; else param.solveWBMmethod = []; end @@ -730,9 +730,8 @@ solution.q = -solution.q; % lb + p <= x <= ub + q elseif strcmp(minNorm, 'QRQP') - buildOptProblemFromModel_param = param; - buildOptProblemFromModel_param.minNorm = param.minNormWBM; - optProblem = buildOptProblemFromModel(model, 0, buildOptProblemFromModel_param); + minNorm = minNormWBM; + optProblem = buildOptProblemFromModel(model, 0, param); solution = solveCobraQP(optProblem); elseif length(minNorm)> 1 || minNorm > 0 @@ -781,7 +780,7 @@ %quadratic optimization will get rid of the loops unless you are maximizing a flux which is %part of a loop. By definition, exchange reactions are not part of these loops, more %properly called stoichiometrically balanced cycles. - solution = solveCobraQP(optProblem2); + solution = solveCobraQP(optProblem2,param); else %this is slow, but more useful than minimizing the Euclidean norm if one is trying to %maximize the flux through a reaction in a loop. e.g. in flux variablity analysis @@ -806,10 +805,10 @@ end end -%TODO fix this Hack in case param.minNorm is used again -if ~isempty(param.solveWBMmethod) - minNorm = param.minNormWBM; -end +% %TODO fix this Hack in case param.minNorm is used again +% if ~isempty(param.solveWBMmethod) +% param.minNorm = param.minNormWBM; +% end %dummy parts of the solution solution.f0 = NaN; diff --git a/src/base/solvers/buildOptProblemFromModel.m b/src/base/solvers/buildOptProblemFromModel.m index 5d47cece35..4ec8f19b31 100644 --- a/src/base/solvers/buildOptProblemFromModel.m +++ b/src/base/solvers/buildOptProblemFromModel.m @@ -72,6 +72,8 @@ % * `.osense`: Objective sense (`-1`: maximise (default); `1`: minimise) % * `.csense`: string with the constraint sense for each row in A ('E', equality, 'G' greater than, 'L' less than). % * `.F`: Positive semidefinite matrix for quadratic part of objective +% OPTIONAL OUTPUT: + [nMet,nRxn]=size(model.S); @@ -294,16 +296,22 @@ % var (cell) – a cell array where names.var{j} contains the name of the % -th variable. optProblem.names.name='optimizeCbModel'; - optProblem.names.obj=model.rxns{model.c~=0}; + if any(model.c~=0) + optProblem.names.obj=model.rxns{model.c~=0}; + else + optProblem.names.obj='noLP'; + end if isfield(model,'ctrs') optProblem.names.con=[model.mets;model.ctrs]; else optProblem.names.con=model.mets; end if isfield(model,'evars') - optProblem.names.con=[model.mets;model.ctrs]; + optProblem.names.var=[model.rxns;model.evars]; else - optProblem.names.con=model.mets; + optProblem.names.var=model.rxns; end end +else + optProblem.names=[]; end \ No newline at end of file diff --git a/src/base/solvers/cplex/setCplexParametersForProblem.m b/src/base/solvers/cplex/setCplexParametersForProblem.m index f53b7f0ca4..5b1d35b2fa 100644 --- a/src/base/solvers/cplex/setCplexParametersForProblem.m +++ b/src/base/solvers/cplex/setCplexParametersForProblem.m @@ -265,6 +265,25 @@ cplexProblem.Param.lpmethod.Cur=4;%BARRIER provided best benchmark performance on Harvetta end +if isfield(solverParams,'multiscale') && solverParams.multiscale==1 && 0 + % Decides how to scale the problem matrix. + % Value Meaning + % -1 No scaling + % 0 Equilibration scaling; default + % 1 More aggressive scaling + % https://www.ibm.com/docs/en/icos/12.10.0?topic=parameters-scale-parameter + cplexProblem.Param.read.scale.Cur = -1; + + % Emphasizes precision in numerically unstable or difficult problems. + % This parameter lets you specify to CPLEX that it should emphasize precision in + % numerically difficult or unstable problems, with consequent performance trade-offs in time and memory. + % Value Meaning + % 0 Do not emphasize numerical precision; default + % 1 Exercise extreme caution in computation + % https://www.ibm.com/docs/en/icos/12.10.0?topic=parameters-numerical-precision-emphasis + cplexProblem.Param.emphasis.numerical.Cur = 1; +end + if isfield(solverParams,'scaind') % Decides how to scale the problem matrix. % Value Meaning diff --git a/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m b/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m index c867a17a3b..8894c2e8bf 100644 --- a/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m +++ b/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m @@ -1,5 +1,8 @@ function [solution, modelOut] = entropicFluxBalanceAnalysis(model, param) -%% TBC +% Entropy maximisation of fluxes (or fluxes and concentrations) subject to +% mass balance, optionally coupling constraints, optionally quadratic +% penalisation of deviation from given fluxes. +% % minimize g.*vf'*(log(vf) -1) + (cf + ci)'*vf % vf,vr,w,x,x0 + g.*vr'*(log(vr) -1) + (cr - ci)'*vr % + f.*x' *(log(x) -1) + u0'*x @@ -219,13 +222,13 @@ %% processing for fluxes [vl,vu,vel,veu,vfl,vfu,vrl,vru,ci,ce,cf,cr,g] = processFluxConstraints(model,param); -if param.debug +if param.debug && 0 %TODO - remove this modelProcessed = model; modelProcessed.lb(model.SConsistentRxnBool)=vl; modelProcessed.lb(~model.SConsistentRxnBool)=vel; modelProcessed.ub(model.SConsistentRxnBool)=vu; modelProcessed.ub(~model.SConsistentRxnBool)=veu; - solutionLP = optimizeCbModel(modelProcessed); + solutionLP = optimizeCbModel(modelProcessed,param); switch solutionLP.stat case 0 diff --git a/src/base/solvers/entropicFBA/processConcConstraints.m b/src/base/solvers/entropicFBA/processConcConstraints.m index bef526d96e..ae486af86c 100644 --- a/src/base/solvers/entropicFBA/processConcConstraints.m +++ b/src/base/solvers/entropicFBA/processConcConstraints.m @@ -67,12 +67,15 @@ %% processing for concentrations if ~isfield(param,'maxConc') - param.maxConc=1e4; + param.maxConc=inf; end if ~isfield(param,'minConc') - param.minConc=1e-4; + param.minConc=0; end - +% %assume units are in mMol +% if ~isfield(param,'concUnit') +% param.concUnit = 10-3; +% end if ~isfield(param,'externalNetFluxBounds') if isfield(model,'dcl') || isfield(model,'dcu') @@ -256,10 +259,7 @@ end end -%assume concentrations are in uMol -if ~isfield(model,'concUnit') - concUnit = 10-6; -end + % Define constants if isfield(model,'gasConstant') && isfield(model,'T') @@ -306,12 +306,3 @@ error('f must all be finite') end end - -% %lower and upper bounds on logarithmic concentration -% pl = -log(param.maxConc*ones(m2,1)); -% if 1 -% pu = log(param.maxConc*ones(m2,1)); -% else -% %All potentials negative -% pu = zeros(m2,1); -% end \ No newline at end of file diff --git a/src/base/solvers/entropicFBA/processFluxConstraints.m b/src/base/solvers/entropicFBA/processFluxConstraints.m index c0944f5764..3886c5d222 100644 --- a/src/base/solvers/entropicFBA/processFluxConstraints.m +++ b/src/base/solvers/entropicFBA/processFluxConstraints.m @@ -56,7 +56,6 @@ % Author(s): Ronan Fleming %% processing for fluxes - if ~isfield(param,'maxUnidirectionalFlux') %try to set the maximum unidirectional flux based on the magnitude of the largest bound but dont have it greater than 1e5 param.maxUnidirectionalFlux=min(1e5,max(abs(model.ub))); @@ -103,8 +102,8 @@ error('internalBounds replaced by other parameter options') end -if isfield(param,'debug') && param.debug - solution_optimizeCbModel = optimizeCbModel(model); +if isfield(param,'debug') && param.debug && 0 + solution_optimizeCbModel = optimizeCbModel(model,param); switch solution_optimizeCbModel.stat case 0 disp(solution_optimizeCbModel.origStat) @@ -225,6 +224,7 @@ vfl = model.vfl; else vfl = max(param.minUnidirectionalFlux,vl); + %vfl = ones(n,1)*param.minUnidirectionalFlux; end if any(vfl<0) @@ -272,6 +272,7 @@ vrl = model.vrl; else vrl = max(param.minUnidirectionalFlux,-vu); + %vrl = ones(n,1)*param.minUnidirectionalFlux; end if any(vrl<0) error('lower bound on reverse flux cannot be less than zero') diff --git a/src/base/solvers/entropicFBA/solveCobraEP.m b/src/base/solvers/entropicFBA/solveCobraEP.m index 1a091ce0b1..9176859443 100644 --- a/src/base/solvers/entropicFBA/solveCobraEP.m +++ b/src/base/solvers/entropicFBA/solveCobraEP.m @@ -97,12 +97,16 @@ % *.objLinear osense*c'*x; % *.objEntropy d.*x'*(log(x) -1); % *.objQuadratic (1/2)*x'*Q*x; -% * .full: Primal sol vector -% * .slack: bl = A*x + s = bu -% * .rcost: Reduced costs, dual sol to :math:`lb <= x <= ub` -% * .dual: dual sol to constraints :math: `A*x ('E' | 'G' | 'L') b` -% -% * .solver: Solver used to solve EP problem +% *.v: n+k ×1 double +% *.vf: n × 1 double +% *.vr: n × 1 double +% *.vt: 1'*vt + 1'*vr +% *.y_N: m x 1 double dual sol to constraints :math: `A*x ('E' | 'G' | 'L') b` +% *.z_dx: 0 +% *.z_vf: n × 1 double dual sol to :math:`lb <= vr <= ub` +% *.z_vr: n × 1 double dual sol to :math:`lb <= vf <= ub` +% *.z_vi: n × 1 double dual sol to :math:`lb <= v <= ub` +% *.z_v: n + k × 1 double dual sol to :math:`lb <= w <= ub` % * .stat: Solver status in standardized form % * 0 - Infeasible problem % * 1 - Optimal sol @@ -112,7 +116,9 @@ % * .origStat: Original status returned by the specific solver % * .origStatText: Original status text returned by the specific solver % * .time: Solve time in seconds -% +% * .solver: Solver used to solve EP problem +% * .epmethod: solver method used e.g. 'CONIC' + % OPTIONAL OUTPUT (from conic optimisation with mosek): % sol.auxPrimal: auxiliary primal variable % sol.auxRcost: dual to auxiliary primal variable @@ -126,9 +132,8 @@ % % EXAMPLE: % -% NOTE: This code is a draft version released for the ELIXIR Fluxomic course and is not yet published and not to be redistributed without express permission of the author. % -% Author(s): Ronan M.T. Fleming, 2021 +% Author(s): Ronan M.T. Fleming, 2024 [problemTypeParams, solverParams] = parseSolverParameters('EP', varargin{:}); @@ -166,7 +171,7 @@ end %% if in debug mode, test to see if the LP part of the problem is feasible -if param.debug +if param.debug && ~isfield(param,'useTestVKSolution') && ~isfield(param,'VKproblem') %avoid LP if being called in driver_optimiseVKmodel switch param.solver case 'pdco' solutionLP2 = solveCobraLP(EPproblem); @@ -657,7 +662,7 @@ % Specify conic part of the problem % https://docs.mosek.com/9.2/toolbox/data-types.html#cones if param.printLevel>1 || param.debug - [~, res] = mosekopt('symbcon',mosekParam); + [~, res] = mosekopt('symbcon echo(0)'); else [~, res] = mosekopt('symbcon echo(0)'); end @@ -797,7 +802,14 @@ end - + if isfield(param,'saveProb') && param.saveProb + formattedTime = datestr(now, 'yyyymmddHHMMSS'); + EP.cmd=cmd; + EP.prob=prob; + EP.param=mosekParam; + save([formattedTime '_EP_probBeforeMosekopt'],"EP"); + end + %call mosek exponential cone solver tic; if 0 diff --git a/src/base/solvers/getSetSolver/changeCobraSolver.m b/src/base/solvers/getSetSolver/changeCobraSolver.m index 19bbc7c041..ceff065934 100644 --- a/src/base/solvers/getSetSolver/changeCobraSolver.m +++ b/src/base/solvers/getSetSolver/changeCobraSolver.m @@ -532,9 +532,9 @@ eval(['CBT_', problemType, '_SOLVER = solverName;']); % validate with a simple problem. if strcmp(solverName,'mosek') && strcmp(problemType,'CLP') || strcmp(problemType,'all') - problem = struct('A',[0 1],'b',0,'c',[1;1],'osense',-1,'lb',[0;0],'ub',[0;0],'csense','E','vartype',['C';'I'],'x0',[0;0]); + problem = struct('A',[0 1],'b',0,'c',[1;1],'osense',-1,'lb',[0;0],'ub',[0;0],'csense','E','vartype',['C';'I'],'x0',[0;0],'names',[]); else - problem = struct('A',[0 1],'b',0,'c',[1;1],'osense',-1,'F',speye(2),'lb',[0;0],'ub',[0;0],'csense','E','vartype',['C';'I'],'x0',[0;0]); + problem = struct('A',[0 1],'b',0,'c',[1;1],'osense',-1,'F',speye(2),'lb',[0;0],'ub',[0;0],'csense','E','vartype',['C';'I'],'x0',[0;0],'names',[]); end try %This is the code that actually tests if a solver is working diff --git a/src/base/solvers/gurobi/setGurobiParam.m b/src/base/solvers/gurobi/setGurobiParam.m index f5f47f0e94..edceb9cf66 100644 --- a/src/base/solvers/gurobi/setGurobiParam.m +++ b/src/base/solvers/gurobi/setGurobiParam.m @@ -13,6 +13,20 @@ param.ScaleFlag=0; end +if isfield(param,'lifted') && param.lifted==1 + param.Aggregate=1; + % Presolve + % Controls the presolve level + % Type: int + % Default value: -1 + % Minimum value: -1 + % Maximum value: 2 + % Controls the presolve level. A value of -1 corresponds to an automatic setting. + % Other options are off (0), conservative (1), or aggressive (2). More aggressive application + % of presolve takes more time, but can sometimes lead to a significantly tighter model. + param.Presolve=0; +end + %backward compatibility if isfield(param,'method') if isempty(param.method) diff --git a/src/base/solvers/mosek/parseMskResult.m b/src/base/solvers/mosek/parseMskResult.m index 8282b0cd1e..05eaf0d088 100644 --- a/src/base/solvers/mosek/parseMskResult.m +++ b/src/base/solvers/mosek/parseMskResult.m @@ -1,4 +1,4 @@ -function [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res) +function [stat,origStat,x,y,yl,yu,z,zl,zu,s,basis,pobjval,dobjval] = parseMskResult(res) %parse the res structure returned from mosek % INPUTS: % res: mosek results structure returned by mosekopt @@ -49,7 +49,7 @@ z = []; zl = []; zu = []; -k = []; +s = []; basis = []; pobjval =[]; dobjval =[]; @@ -172,7 +172,7 @@ zu=res.sol.itr.sux; %dual to x <= bux if isfield(res.sol.itr,'doty') % Dual variables to affine conic constraints - k = res.sol.itr.doty; + s = res.sol.itr.doty; end pobjval = res.sol.itr.pobjval; dobjval = res.sol.itr.dobjval; @@ -198,14 +198,14 @@ zu=res.sol.bas.sux; %dual to x <= bux if isfield(res.sol.bas,'s') % Dual variables to affine conic constraints - k = res.sol.bas.s; + s = res.sol.bas.s; end %https://docs.mosek.com/10.0/toolbox/advanced-hotstart.html - bas.skc = res.sol.bas.skc; - bas.skx = res.sol.bas.skx; - bas.xc = res.sol.bas.xc; - bas.xx = res.sol.bas.xx; + basis.skc = res.sol.bas.skc; + basis.skx = res.sol.bas.skx; + basis.xc = res.sol.bas.xc; + basis.xx = res.sol.bas.xx; pobjval = res.sol.bas.pobjval; dobjval = res.sol.bas.dobjval; otherwise @@ -213,6 +213,7 @@ end otherwise accessSolution = 'dontAccess'; + origStat = -1; end if strcmp(accessSolution,'dontAccess') @@ -223,11 +224,11 @@ case {'DUAL_INFEASIBLE_CER','MSK_SOL_STA_DUAL_INFEAS_CER','MSK_SOL_STA_NEAR_DUAL_INFEAS_CER'} stat=2; % Unbounded solution origStat = [origStat ' & ' res.rcodestr]; - case {'UNKNOWN','PRIM_ILLPOSED_CER','DUAL_ILLPOSED_CER','PRIM_FEAS','DUAL_FEAS','PRIM_AND_DUAL_FEAS'} + case {'UNKNOWN','PRIM_ILLPOSED_CER','DUAL_ILLPOSED_CER','PRIM_FEAS','DUAL_FEAS','PRIM_AND_DUAL_FEAS','DUAL_FEASIBLE'} stat=-1; %some other problem origStat = [origStat ' & ' res.rcodestr]; otherwise - warning(['Unrecognised res.sol.bas.solsta: ' origStat]) + warning(['Unrecognised res.sol.bas.solsta or res.sol.itr.solsta: ' origStat]) stat=-1; %some other problem fprintf('%s\n',res.rcode) fprintf('%s\n',res.rmsg) diff --git a/src/base/solvers/mosek/setMosekParam.m b/src/base/solvers/mosek/setMosekParam.m index d867594482..8b786f0b11 100644 --- a/src/base/solvers/mosek/setMosekParam.m +++ b/src/base/solvers/mosek/setMosekParam.m @@ -274,20 +274,6 @@ end end -% param.lpmethod='BARRIER'; -% param.qpmethod='BARRIER'; -% method = param.lpmethod; -% case 'ibm_cplex' -% param.lpmethod='BARRIER'; -% param.qpmethod='BARRIER'; -% method = param.lpmethod; -% case 'mosek' -% method = 'FREE'; -% method = 'INTPNT'; -% method = 'CONIC'; -% - - %backward compatibility if isfield(param,'method') if isempty(param.method) @@ -324,6 +310,16 @@ param.MSK_IPAR_OPTIMIZER=['MSK_OPTIMIZER_' param.clpmethod]; end end + + % MSK_IPAR_INTPNT_REGULARIZATION_USE + % Description:Controls whether regularization is allowed. + % Possible Values: MSK_ON Switch the option on. + % MSK_OFF Switch the option off. + % Default value: MSK_ON + if ~isfield(param,'MSK_IPAR_INTPNT_REGULARIZATION_USE') + param.MSK_IPAR_INTPNT_REGULARIZATION_USE='MSK_ON'; + end + case {'EP'} if isfield(param,'epmethod') if contains(param.epmethod,'MSK_OPTIMIZER_') @@ -366,6 +362,21 @@ % param.MSK_IPAR_SIM_SCALING_METHOD='MSK_SCALING_METHOD_FREE'; end + +if ~isfield(param,'MSK_IPAR_LOG_FEAS_REPAIR') && isfield(param,'repairInfeasibility') + % MSK_IPAR_LOG_FEAS_REPAIR + % Controls the amount of output printed when performing feasibility repair. A value higher than one means extensive logging. + % Default + % 1 + % Accepted + % [0; +inf] + % Example + % MSK_putintparam(task, MSK_IPAR_LOG_FEAS_REPAIR, 1) + % Groups + % Output information, Logging + param.MSK_IPAR_LOG_FEAS_REPAIR = param.repairInfeasibility; +end + % Remove outer function specific parameters to avoid crashing solver interfaces mosekParam = mosekParamStrip(param); if 0 diff --git a/src/base/solvers/param/getCobraSolverParams.m b/src/base/solvers/param/getCobraSolverParams.m index 2ea0d2f56b..92daa82e30 100644 --- a/src/base/solvers/param/getCobraSolverParams.m +++ b/src/base/solvers/param/getCobraSolverParams.m @@ -131,7 +131,7 @@ valDef.feasTol = 1e-6; % (primal) feasibility tolerance valDef.optTol = 1e-6; % (dual) optimality tolerance valDef.solver='mosek'; - valDef.epmethod='CONIC'; + case 'CLP' % This is never used elsewhere except for parameter setting loop % for backward compatibility diff --git a/src/base/solvers/param/getCobraSolverParamsOptionsForType.m b/src/base/solvers/param/getCobraSolverParamsOptionsForType.m index cee1c4cc0c..570eac374c 100644 --- a/src/base/solvers/param/getCobraSolverParamsOptionsForType.m +++ b/src/base/solvers/param/getCobraSolverParamsOptionsForType.m @@ -36,13 +36,11 @@ 'solver', ... % solver to use (overriding set solver) 'debug', ... % run debgugging code 'logFile', ... % file (location) to write logs to - 'lifting', ... % whether to lift a problem - 'method'}; % solver method: -1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent, 5 = Network Solver(if supported by the solver) + }; % case 'QP' paramNames = {'multiscale'... % true if problem is multiscale 'problemType'... % problem type 'verify',... % verify that it is a suitable QP problem - 'method', ... % solver method: -1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent, 5 = Network Solver(if supported by the solver) 'printLevel', ... % print level 'saveInput', ... % save the input to a file (specified) 'debug', ... % run debgugging code @@ -55,7 +53,6 @@ paramNames = {'multiscale'... % true if problem is multiscale 'problemType'... % problem type 'verify',... % verify that it is a suitable EP problem - 'method', ... % solver method: -1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent, 5 = Network Solver(if supported by the solver) 'printLevel', ... % print level 'debug', ... % run debgugging code 'feasTol',... % feasibility tolerance diff --git a/src/base/solvers/param/parseSolverParameters.m b/src/base/solvers/param/parseSolverParameters.m index 2eb8eae37f..930ea6b602 100644 --- a/src/base/solvers/param/parseSolverParameters.m +++ b/src/base/solvers/param/parseSolverParameters.m @@ -101,36 +101,4 @@ % otherwise use the default parameter param.(defaultParams{i,1}) = defaultParams{i,2}; end -end - -if 0 - %move following set of parameters from solverOnlyParams to param - if isfield(solverOnlyParams,'maxConc') - param.maxConc = solverOnlyParams.maxConc; - solverOnlyParams = rmfield(solverOnlyParams,'maxConc'); - end - if isfield(solverOnlyParams,'method') - param.method = solverOnlyParams.method; - solverOnlyParams = rmfield(solverOnlyParams,'method'); - end - if isfield(solverOnlyParams,'maxUnidirectionalFlux') - param.maxUnidirectionalFlux = solverOnlyParams.maxUnidirectionalFlux; - solverOnlyParams = rmfield(solverOnlyParams,'maxUnidirectionalFlux'); - end - if isfield(solverOnlyParams,'minUnidirectionalFlux') - param.minUnidirectionalFlux = solverOnlyParams.minUnidirectionalFlux; - solverOnlyParams = rmfield(solverOnlyParams,'minUnidirectionalFlux'); - end - if isfield(solverOnlyParams,'internalNetFluxBounds') - param.internalNetFluxBounds = solverOnlyParams.internalNetFluxBounds; - solverOnlyParams = rmfield(solverOnlyParams,'internalNetFluxBounds'); - end - if isfield(solverOnlyParams,'externalNetFluxBounds') - param.externalNetFluxBounds = solverOnlyParams.externalNetFluxBounds; - solverOnlyParams = rmfield(solverOnlyParams,'externalNetFluxBounds'); - end - if isfield(solverOnlyParams,'rounding') - param.rounding = solverOnlyParams.rounding; - solverOnlyParams = rmfield(solverOnlyParams,'rounding'); - end end \ No newline at end of file diff --git a/src/base/solvers/solveCobraLP.m b/src/base/solvers/solveCobraLP.m index e0eb9c9d5c..35667aa406 100644 --- a/src/base/solvers/solveCobraLP.m +++ b/src/base/solvers/solveCobraLP.m @@ -144,7 +144,7 @@ end % support for lifting of ill-scaled models -if problemTypeParams.lifting == 1 +if isfield(solverParams,'lifting') && solverParams.lifting == 1 largeNb = 1e4; % suitable for double precision solvers [LPproblem] = reformulate(LPproblem, largeNb, printLevel); end @@ -180,6 +180,9 @@ if ~isfield(LPproblem, 'modelID') LPproblem.modelID = 'aModelID'; end +if ~isfield(LPproblem, 'names') + LPproblem.names = []; +end %too time consuming % if any(~isfinite(LPproblem.A),'all') @@ -197,7 +200,7 @@ % extract the problem from the structure -[A, b, c, lb, ub, csense, osense, modelID] = deal(sparse(LPproblem.A), LPproblem.b, LPproblem.c, LPproblem.lb, LPproblem.ub, LPproblem.csense, LPproblem.osense, LPproblem.modelID); +[A, b, c, lb, ub, csense, osense, modelID, names] = deal(sparse(LPproblem.A), LPproblem.b, LPproblem.c, LPproblem.lb, LPproblem.ub, LPproblem.csense, LPproblem.osense, LPproblem.modelID,LPproblem.names); if isfield(LPproblem,'basis') && ~isempty(LPproblem.basis) basis = LPproblem.basis; @@ -688,13 +691,23 @@ prob.sol.bas.xx = basis.xx; end - if param.debug - probBeforeMosekopt = prob; - save('probBeforeMosekopt','probBeforeMosekopt'); + if isfield(param,'saveProb') && param.saveProb + formattedTime = datestr(now, 'yyyymmddHHMMSS'); + LP.cmd=cmd; + LP.prob=prob; + LP.param=mosekParam; + save([formattedTime '_LP_probBeforeMosekopt'],"LP"); end - - [rcode,res] = mosekopt(cmd,prob,mosekParam); + if isfield(param,'debug') && param.debug==1 + prob.names = names; + end + + if isfield(param,'defaultSolverParam') && param.defaultSolverParam + [rcode,res] = mosekopt(cmd,prob); + else + [rcode,res] = mosekopt(cmd,prob,mosekParam); + end if isfield(param,'lpmethod') lpmethod = param.lpmethod; else diff --git a/src/base/solvers/solveCobraQP.m b/src/base/solvers/solveCobraQP.m index 2670b7408f..12de31c3a5 100644 --- a/src/base/solvers/solveCobraQP.m +++ b/src/base/solvers/solveCobraQP.m @@ -406,6 +406,13 @@ %https://docs.mosek.com/latest/toolbox/data-types.html#prob [prob.qosubi,prob.qosubj,prob.qoval]=find(F); + if isfield(param,'saveProb') && param.saveProb + formattedTime = datestr(now, 'yyyymmddHHMMSS'); + QP.cmd=cmd; + QP.prob=prob; + QP.param=mosekParam; + save([formattedTime '_QP_probBeforeMosekopt'],"QP"); + end [rcode,res] = mosekopt(cmd,prob,mosekParam); @@ -981,7 +988,7 @@ fprintf('%s\n',['[' solver '] reports ' solution.origStat ' but Primal optimality condition in solveCobraQP not satisfied, residual = ' num2str(tmp1) ', while feasTol = ' num2str(problemTypeParams.feasTol)]) else if problemTypeParams.printLevel > 0 - fprintf(['\n > [' solver '] Primal optimality condition in solveCobraQP satisfied.']); + fprintf(['\n > [' solver '] Primal optimality condition in solveCobraQP satisfied.\n']); end end end From ab3ca80f60bfc00f7061a65aaffc84c2fe066b43 Mon Sep 17 00:00:00 2001 From: Ronan Fleming Date: Thu, 7 Nov 2024 09:28:14 +0000 Subject: [PATCH 3/3] mosek maxCardinalityConservationVector working --- src/base/solvers/cardOpt/sparseLP/optimizeCardinality.m | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/base/solvers/cardOpt/sparseLP/optimizeCardinality.m b/src/base/solvers/cardOpt/sparseLP/optimizeCardinality.m index 45dad396e5..b31014d0bf 100644 --- a/src/base/solvers/cardOpt/sparseLP/optimizeCardinality.m +++ b/src/base/solvers/cardOpt/sparseLP/optimizeCardinality.m @@ -91,9 +91,6 @@ [feasTol] = getCobraSolverParams('LP', 'feasTol'); global CBT_LP_SOLVER -if strcmp('mosek',CBT_LP_SOLVER) - error('mosek not working with optCard -TODO debug with testOptimizeCbModel') -end %% Check inputs if ~exist('param','var') || isempty(param)