diff --git a/desdeo/problem/cvxpy_evaluator.py b/desdeo/problem/cvxpy_evaluator.py index 090297ac..883d2f48 100644 --- a/desdeo/problem/cvxpy_evaluator.py +++ b/desdeo/problem/cvxpy_evaluator.py @@ -34,6 +34,7 @@ class CVXPYEvaluator: variables: dict[str, cp.Variable] parameters: dict[str, cp.Parameter] constraints: dict[str, cp.Constraint] + constraint_expressions: dict[str, cp.Expression] objective_functions: dict[str, cp.Expression] scalarizations: dict[str, cp.Expression] extra_functions: dict[str, cp.Expression] @@ -55,6 +56,7 @@ def __init__(self, problem: Problem): self.variables = {} self.parameters = {} self.constraints = {} + self.constraint_expressions = {} self.objective_expr = None self.problem_model = None @@ -98,65 +100,9 @@ def init_variables(self, problem: Problem) -> dict[str, cp.Variable]: Returns: dict[str, cp.Variable]: the variables for the problem. """ - variables = {} for var in problem.variables: - if isinstance(var, Variable): - # handle regular variables - lowerbound = var.lowerbound - upperbound = var.upperbound - - # Set bounds - bounds = None - if lowerbound is not None or upperbound is not None: - lb = lowerbound if lowerbound is not None else -np.inf - ub = upperbound if upperbound is not None else np.inf - bounds = [lb, ub] - - # figure out the variable type - match var.variable_type: - case VariableTypeEnum.integer: - # variable is integer - cv_var = cp.Variable(integer=True, bounds=bounds, name=var.symbol) - case VariableTypeEnum.real: - # variable is real - cv_var = cp.Variable(bounds=bounds, name=var.symbol) - case VariableTypeEnum.binary: - cv_var = cp.Variable(boolean=True, bounds=bounds, name=var.symbol) - case _: - msg = f"Could not figure out the type for variable {var}." - raise CVXPYEvaluatorError(msg) - - variables[var.symbol] = cv_var - - elif isinstance(var, TensorVariable): - # handle tensor variables, i.e., vectors etc.. - lowerbounds = var.get_lowerbound_values() if var.lowerbounds is not None else None - upperbounds = var.get_upperbound_values() if var.upperbounds is not None else None - - # Set bounds - bounds = None - if lowerbounds is not None or upperbounds is not None: - lb = np.array(lowerbounds) if lowerbounds is not None else -np.inf - ub = np.array(upperbounds) if upperbounds is not None else np.inf - bounds = [lb, ub] - - # figure out the variable type - match var.variable_type: - case VariableTypeEnum.integer: - # variable is integer - cv_var = cp.Variable(shape=tuple(var.shape), integer=True, bounds=bounds, name=var.symbol) - case VariableTypeEnum.real: - # variable is real - cv_var = cp.Variable(shape=tuple(var.shape), bounds=bounds, name=var.symbol) - case VariableTypeEnum.binary: - cv_var = cp.Variable(shape=tuple(var.shape), boolean=True, bounds=bounds, name=var.symbol) - case _: - msg = f"Could not figure out the type for variable {var}." - raise CVXPYEvaluatorError(msg) - - variables[var.symbol] = cv_var - - return variables + self.add_variable(var) + return self.variables def init_parameters( self, problem: Problem @@ -221,22 +167,9 @@ def init_objectives(self, problem: Problem) -> dict[str, cp.Expression]: Returns: dict[str, cp.Expression]: dict containing the objective functions. """ - objective_functions: dict[str, cp.Expression] = {} for obj in problem.objectives: - expr = self.parse(obj.func, callback=self.get_expression_by_name) - if isinstance(expr, (int, float)): - warnings.warn( - "One or more of the problem objectives seems to be a constant.", - CVXPYEvaluatorWarning, - stacklevel=2, - ) - - objective_functions[obj.symbol] = expr - - # the obj.symbol_min objectives are used when optimizing and building scalarizations etc... - objective_functions[f"{obj.symbol}_min"] = -expr if obj.maximize else expr - - return objective_functions + self.add_objective(obj) + return self.objective_functions def init_constraints(self, problem: Problem) -> dict[str, cp.Constraint]: """Add constraint expressions to a CVXPY problem. @@ -250,23 +183,9 @@ def init_constraints(self, problem: Problem) -> dict[str, cp.Constraint]: Returns: dict[str, cp.Constraint]: dict of constraints keyed by symbol. """ - constraints = {} for cons in problem.constraints: - expr = self.parse(cons.func, callback=self.get_expression_by_name) - - match con_type := cons.cons_type: - case ConstraintTypeEnum.LTE: - # constraints in DESDEO are defined such that they must be less than zero - constraint = expr <= 0 - case ConstraintTypeEnum.EQ: - constraint = expr == 0 - case _: - msg = f"Constraint type of {con_type} not supported. Must be one of {ConstraintTypeEnum}." - raise CVXPYEvaluatorError(msg) - - constraints[cons.symbol] = constraint - - return constraints + self.add_constraint(cons) + return self.constraints def init_scalarizations(self, problem: Problem) -> dict[str, cp.Expression]: """Add scalarization expressions to a CVXPY evaluator. @@ -300,15 +219,15 @@ def add_constraint(self, constraint: Constraint) -> cp.Constraint: cp.Constraint: The CVXPY constraint that was added. """ expr = self.parse(constraint.func, self.get_expression_by_name) + self.constraint_expressions[constraint.symbol] = expr - match con_type := constraint.cons_type: + match constraint.cons_type: case ConstraintTypeEnum.LTE: - # constraints in DESDEO are defined such that they must be less than zero cvxpy_constraint = expr <= 0 case ConstraintTypeEnum.EQ: cvxpy_constraint = expr == 0 case _: - msg = f"Constraint type of {con_type} not supported. Must be one of {ConstraintTypeEnum}." + msg = f"Constraint type of {constraint.cons_type} not supported. Must be one of {ConstraintTypeEnum}." raise CVXPYEvaluatorError(msg) self.constraints[constraint.symbol] = cvxpy_constraint @@ -369,7 +288,7 @@ def add_variable(self, var: Variable | TensorVariable) -> cp.Variable: if lowerbound is not None or upperbound is not None: lb = lowerbound if lowerbound is not None else -cp.inf ub = upperbound if upperbound is not None else cp.inf - bounds = (lb, ub) + bounds = [lb, ub] # figure out the variable type match var.variable_type: @@ -395,7 +314,7 @@ def add_variable(self, var: Variable | TensorVariable) -> cp.Variable: if lowerbounds is not None or upperbounds is not None: lb = np.array(lowerbounds) if lowerbounds is not None else -np.inf ub = np.array(upperbounds) if upperbounds is not None else np.inf - bounds = (lb, ub) + bounds = [lb, ub] # figure out the variable type match var.variable_type: @@ -457,7 +376,7 @@ def get_values(self) -> dict[str, float | int | bool | list[float] | list[int] | Raises: CVXPYEvaluatorError: if the problem has not been solved yet. """ - if self.problem_model is None or self.problem_model.status not in ["optimal", "optimal_inaccurate"]: + if self.problem_model is None or self.problem_model.status not in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE]: msg = "Problem has not been solved yet or did not achieve optimal status." raise CVXPYEvaluatorError(msg) @@ -487,8 +406,7 @@ def get_values(self) -> dict[str, float | int | bool | list[float] | list[int] | if self.problem.constraints is not None: for con in self.problem.constraints: - constraint_expr = self.parse(con.func, callback=self.get_expression_by_name) - result_dict[con.symbol] = constraint_expr.value + result_dict[con.symbol] = self.constraint_expressions[con.symbol].value return result_dict @@ -497,8 +415,6 @@ def set_optimization_target(self, target: str): Args: target (str): a str representing a symbol. Needs to match an objective function or scalarization - function already found in the evaluator. - maximize (bool): If true, the target function is maximized instead of minimized. Raises: CVXPYEvaluatorError: the given target was not found in the evaluator. @@ -509,7 +425,7 @@ def set_optimization_target(self, target: str): elif target in self.scalarizations: maximize = False else: - msg = f"The gurobipy model has no objective or scalarization named {target}." + msg = f"The CVXPY model has no objective or scalarization named {target}." raise CVXPYEvaluatorError(msg) obj_expr = self.get_expression_by_name(target) @@ -532,7 +448,8 @@ def solve(self, **kwargs): if self.problem_model.is_dcp(): self.problem_model.solve(**kwargs) elif self.problem_model.is_dgp(): - self.problem_model.solve(gp=True, **kwargs) + kwargs["gp"] = True + self.problem_model.solve(**kwargs) else: warnings.warn( "The problem does not appear to be DCP or DGP. CVXPY may not be able to solve it.", diff --git a/desdeo/problem/testproblems/simple_problem.py b/desdeo/problem/testproblems/simple_problem.py index 58c1370f..63be57ac 100644 --- a/desdeo/problem/testproblems/simple_problem.py +++ b/desdeo/problem/testproblems/simple_problem.py @@ -250,6 +250,8 @@ def simple_constrained_quadratic_tensor_test_problem(dqp=False) -> Problem: constants=[mmult, bvector], constraints=[cons], objectives=[obj], + is_twice_differentiable=True, + is_convex=True, ) diff --git a/desdeo/tools/cvxpy_solver_interfaces.py b/desdeo/tools/cvxpy_solver_interfaces.py index eac4072f..a8b71f2b 100644 --- a/desdeo/tools/cvxpy_solver_interfaces.py +++ b/desdeo/tools/cvxpy_solver_interfaces.py @@ -1,5 +1,6 @@ """Defines solver interfaces for cvxpy.""" +import cvxpy as cp from pydantic import BaseModel, Field from desdeo.problem import ( @@ -79,18 +80,18 @@ def parse_cvxpy_optimizer_results(problem: Problem, evaluator: CVXPYEvaluator) - ) lagrange_multipliers = None - success = evaluator.problem_model.status in {"optimal", "optimal_inaccurate"} - if evaluator.problem_model.status == "optimal": + success = evaluator.problem_model.status in {cp.OPTIMAL, cp.OPTIMAL_INACCURATE} + if evaluator.problem_model.status == cp.OPTIMAL: status = "Optimal solution found." - elif evaluator.problem_model.status == "optimal_inaccurate": + elif evaluator.problem_model.status == cp.OPTIMAL_INACCURATE: status = "Optimal solution found (inaccurate)." - elif evaluator.problem_model.status == "infeasible": + elif evaluator.problem_model.status == cp.INFEASIBLE: status = "Problem is infeasible." - elif evaluator.problem_model.status == "unbounded": + elif evaluator.problem_model.status == cp.UNBOUNDED: status = "Problem is unbounded." - elif evaluator.problem_model.status == "infeasible_inaccurate": + elif evaluator.problem_model.status == cp.INFEASIBLE_INACCURATE: status = "Problem is infeasible (inaccurate)." - elif evaluator.problem_model.status == "unbounded_inaccurate": + elif evaluator.problem_model.status == cp.UNBOUNDED_INACCURATE: status = "Problem is unbounded (inaccurate)." else: status = f"Optimization ended with status: {evaluator.problem_model.status}" @@ -144,3 +145,18 @@ def solve(self, target: str) -> SolverResults: self.evaluator.set_optimization_target(target) self.evaluator.solve(**self.solve_options) return parse_cvxpy_optimizer_results(self.problem, self.evaluator) + + +def check_cvxpy_suitability(problem: Problem) -> bool: + """Checks whether a problem is suitable for being solved with CVXPY.""" + try: + evaluator = CVXPYEvaluator(problem) + for obj in problem.objectives: + evaluator.set_optimization_target(obj.symbol) + if not (evaluator.problem_model.is_dcp() or evaluator.problem_model.is_dgp()): + return False + break + else: + return True + except Exception: + return False diff --git a/desdeo/tools/gurobipy_solver_interfaces.py b/desdeo/tools/gurobipy_solver_interfaces.py index 56a5702e..bb26199b 100644 --- a/desdeo/tools/gurobipy_solver_interfaces.py +++ b/desdeo/tools/gurobipy_solver_interfaces.py @@ -1,5 +1,8 @@ """Defines solver interfaces for gurobipy.""" +import io +import sys + import gurobipy as gp from desdeo.problem import ( @@ -276,3 +279,28 @@ def solve(self, target: str) -> SolverResults: self.evaluator.set_optimization_target(target) self.evaluator.model.optimize() return parse_gurobipy_optimizer_results(self.problem, self.evaluator) + + +def check_gurobi_license(): + """Check if Gurobi is using a full license (not trial). + + Returns: + True if using full academic/commercial license + False if using trial license or no license found + """ + captured_output = io.StringIO() + original_stdout = sys.stdout + + try: + sys.stdout = captured_output + with gp.Env(empty=True) as env: + env.setParam("OutputFlag", 1) + env.start() + sys.stdout = original_stdout + + output = captured_output.getvalue() + return "Restricted license - for non-production use only" not in output # noqa: TRY300 + + except Exception: + sys.stdout = original_stdout + return False diff --git a/desdeo/tools/utils.py b/desdeo/tools/utils.py index 2a6b327d..f7e6bed7 100644 --- a/desdeo/tools/utils.py +++ b/desdeo/tools/utils.py @@ -14,8 +14,9 @@ numpy_array_to_objective_dict, variable_dimension_enumerate, ) +from desdeo.tools.cvxpy_solver_interfaces import CVXPYSolver, CVXPYSolverOptions, check_cvxpy_suitability from desdeo.tools.generics import BaseSolver -from desdeo.tools.gurobipy_solver_interfaces import GurobipySolver, PersistentGurobipySolver +from desdeo.tools.gurobipy_solver_interfaces import GurobipySolver, PersistentGurobipySolver, check_gurobi_license from desdeo.tools.ng_solver_interfaces import NevergradGenericOptions, NevergradGenericSolver from desdeo.tools.proximal_solver import ProximalSolver from desdeo.tools.pyomo_solver_interfaces import ( @@ -69,6 +70,10 @@ "constructor": PersistentGurobipySolver, "options": None, }, + "cvxpy": { + "constructor": CVXPYSolver, + "options": CVXPYSolverOptions, + }, } @@ -122,8 +127,12 @@ def find_compatible_solvers(problem: Problem) -> list[BaseSolver]: ): solvers.append(available_solvers["pyomo_ipopt"]["constructor"]) # ipopt has to be installed + # check if the problem is convex or log-log convex + if check_cvxpy_suitability(problem): + solvers.append(available_solvers["cvxpy"]["constructor"]) + # check if the problem is linear - if problem.is_linear: + if problem.is_linear and check_gurobi_license(): solvers.append(available_solvers["gurobipy"]["constructor"]) if problem.is_linear and shutil.which("gurobi"): solvers.append(available_solvers["pyomo_gurobi"]["constructor"]) # gurobi has to be installed @@ -197,9 +206,13 @@ def guess_best_solver(problem: Problem) -> BaseSolver: return available_solvers["proximal"]["constructor"] # check if the problem is linear - if problem.is_linear: + if problem.is_linear and check_gurobi_license(): return available_solvers["gurobipy"]["constructor"] + # check if the problem is convex or log-log convex + if check_cvxpy_suitability(problem): + return available_solvers["cvxpy"]["constructor"] + # check if the problem is differentiable and if it is mixed integer if ( problem.is_twice_differentiable diff --git a/tests/test_tools_utils.py b/tests/test_tools_utils.py index fb9d9a67..637c9f32 100644 --- a/tests/test_tools_utils.py +++ b/tests/test_tools_utils.py @@ -5,7 +5,12 @@ import pytest from fixtures import dtlz2_5x_3f_data_based # noqa: F401 -from desdeo.problem.testproblems import dtlz2, re21, river_pollution_problem +from desdeo.problem.testproblems import ( + dtlz2, + re21, + river_pollution_problem, + simple_constrained_quadratic_tensor_test_problem, +) from desdeo.tools.utils import ( available_solvers, find_compatible_solvers, @@ -52,10 +57,28 @@ def test_find_compatible_solvers(): else: assert len(solvers) == 3 + problem = simple_constrained_quadratic_tensor_test_problem(dqp=True) + + solvers = find_compatible_solvers(problem) + + correct_solvers = [ + available_solvers["pyomo_ipopt"]["constructor"], + available_solvers["cvxpy"]["constructor"], + ] + + # check that the solvers found are the correct ones + if shutil.which("ipopt"): + assert len(solvers) == 2 + assert all(solver in correct_solvers for solver in solvers) and all( + solver in solvers for solver in correct_solvers + ) + else: + assert len(solvers) == 1 + @pytest.mark.utils def test_payoff_dtlz2(): """Tests the payoff-table method with the dtlz2 problem.""" problem = dtlz2(6, 4) - ideal, nadir = payoff_table_method(problem) + ideal, nadir = payoff_table_method(problem) # noqa: RUF059