diff --git a/discrete_optimization/coloring/solvers/coloring_lp_solvers.py b/discrete_optimization/coloring/solvers/coloring_lp_solvers.py index e39c8cc28..793cb64f8 100644 --- a/discrete_optimization/coloring/solvers/coloring_lp_solvers.py +++ b/discrete_optimization/coloring/solvers/coloring_lp_solvers.py @@ -11,6 +11,7 @@ import mip import networkx as nx from mip import BINARY, INTEGER, xsum +from ortools.math_opt.python import mathopt from discrete_optimization.coloring.coloring_model import ( ColoringProblem, @@ -34,6 +35,7 @@ GurobiMilpSolver, MilpSolver, MilpSolverName, + OrtoolsMathOptMilpSolver, PymipMilpSolver, ) @@ -483,3 +485,149 @@ def init_model(self, **kwargs: Any) -> None: self.description_constraint["constraints_neighbors"] = { "descr": "no neighbors can have same color" } + + +class ColoringLPMathOpt(OrtoolsMathOptMilpSolver, _BaseColoringLP): + """Coloring LP solver based on pymip library. + + Note: + Gurobi and CBC are available as backend solvers. + + + Attributes: + problem (ColoringProblem): coloring problem instance to solve + params_objective_function (ParamsObjectiveFunction): objective function parameters + (however this is just used for the ResultStorage creation, not in the optimisation) + + """ + + hyperparameters = _BaseColoringLP.hyperparameters + + problem: ColoringProblem + + def init_model(self, **kwargs: Any) -> None: + kwargs = self.complete_with_default_hyperparameters(kwargs) + greedy_start = kwargs["greedy_start"] + use_cliques = kwargs["use_cliques"] + if greedy_start: + logger.info("Computing greedy solution") + greedy_solver = GreedyColoring( + self.problem, + params_objective_function=self.params_objective_function, + ) + sol = greedy_solver.solve( + strategy=NXGreedyColoringMethod.best + ).get_best_solution() + if sol is None: + raise RuntimeError( + "greedy_solver.solve(strategy=NXGreedyColoringMethod.best).get_best_solution() " + "should not be None." + ) + if not isinstance(sol, ColoringSolution): + raise RuntimeError( + "greedy_solver.solve(strategy=NXGreedyColoringMethod.best).get_best_solution() " + "should be a ColoringSolution." + ) + self.start_solution = sol + else: + logger.info("Get dummy solution") + self.start_solution = self.problem.get_dummy_solution() + nb_colors = self.start_solution.nb_color + nb_colors_subset = nb_colors + if self.problem.use_subset: + nb_colors_subset = self.problem.count_colors(self.start_solution.colors) + nb_colors = self.problem.count_colors_all_index(self.start_solution.colors) + + if nb_colors is None: + raise RuntimeError("self.start_solution.nb_color should not be None.") + color_model = mathopt.Model(name="color") + colors_var = {} + range_node = self.nodes_name + range_color = range(nb_colors) + range_color_subset = range(nb_colors_subset) + range_color_per_node = {} + for node in self.nodes_name: + rng = self.get_range_color( + node_name=node, + range_color_subset=range_color_subset, + range_color_all=range_color, + ) + for color in rng: + colors_var[node, color] = color_model.add_binary_variable( + name="x_" + str((node, color)) + ) + range_color_per_node[node] = set(rng) + one_color_constraints = {} + for n in range_node: + one_color_constraints[n] = color_model.add_linear_constraint( + mathopt.LinearSum(colors_var[n, c] for c in range_color_per_node[n]) + == 1 + ) + cliques = [] + g = self.graph.to_networkx() + if use_cliques: + for c in nx.algorithms.clique.find_cliques(g): + cliques += [c] + cliques = sorted(cliques, key=lambda x: len(x), reverse=True) + else: + cliques = [[e[0], e[1]] for e in g.edges()] + cliques_constraint: Dict[Union[int, Tuple[int, int]], Any] = {} + index_c = 0 + opt = color_model.add_integer_variable(lb=0, ub=nb_colors, name="nb_colors") + color_model.minimize(opt) + if use_cliques: + for c in cliques[:100]: + cliques_constraint[index_c] = color_model.add_linear_constraint( + mathopt.LinearSum( + (color_i + 1) * colors_var[node, color_i] + for node in c + for color_i in range_color_per_node[node] + ) + >= sum([i + 1 for i in range(len(c))]) + ) + cliques_constraint[(index_c, 1)] = color_model.add_linear_constraint( + mathopt.LinearSum( + colors_var[node, color_i] + for node in c + for color_i in range_color_per_node[node] + ) + <= opt + ) + index_c += 1 + edges = g.edges() + constraints_neighbors = {} + for e in edges: + for c in range_color_per_node[e[0]]: + if c in range_color_per_node[e[1]]: + constraints_neighbors[ + (e[0], e[1], c) + ] = color_model.add_linear_constraint( + colors_var[e[0], c] + colors_var[e[1], c] <= 1 + ) + for n in range_node: + color_model.add_linear_constraint( + mathopt.LinearSum( + (color_i + 1) * colors_var[n, color_i] + for color_i in range_color_per_node[n] + ) + <= opt + ) + self.model = color_model + self.variable_decision = {"colors_var": colors_var} + self.constraints_dict = { + "one_color_constraints": one_color_constraints, + "constraints_neighbors": constraints_neighbors, + } + self.description_variable_description = { + "colors_var": { + "shape": (self.number_of_nodes, nb_colors), + "type": bool, + "descr": "for each node and each color," " a binary indicator", + } + } + self.description_constraint["one_color_constraints"] = { + "descr": "one and only one color " "should be assignated to a node" + } + self.description_constraint["constraints_neighbors"] = { + "descr": "no neighbors can have same color" + } diff --git a/discrete_optimization/generic_tools/lp_tools.py b/discrete_optimization/generic_tools/lp_tools.py index 56f3a6cab..cba679b42 100644 --- a/discrete_optimization/generic_tools/lp_tools.py +++ b/discrete_optimization/generic_tools/lp_tools.py @@ -1,13 +1,15 @@ # Copyright (c) 2022 AIRBUS and its affiliates. # This source code is licensed under the MIT license found in the # LICENSE file in the root directory of this source tree. - +import copy +import datetime import logging from abc import abstractmethod from enum import Enum from typing import Any, Callable, List, Optional, Tuple, Union import mip +from ortools.math_opt.python import mathopt from discrete_optimization.generic_tools.callbacks.callback import ( Callback, @@ -20,6 +22,9 @@ ) from discrete_optimization.generic_tools.do_solver import SolverDO from discrete_optimization.generic_tools.exceptions import SolveEarlyStop +from discrete_optimization.generic_tools.hyperparameters.hyperparameter import ( + EnumHyperparameter, +) from discrete_optimization.generic_tools.result_storage.multiobj_utils import ( TupleFitness, ) @@ -282,6 +287,213 @@ def nb_solutions(self) -> int: return self.model.num_solutions +class OrtoolsMathOptMilpSolver(MilpSolver): + """Milp solver wrapping a solver from pymip library.""" + + hyperparameters = [ + EnumHyperparameter( + name="mathopt_solver_type", + enum=mathopt.SolverType, + ) + ] + + def get_var_value_for_ith_solution(self, var: Any, i: int) -> float: + pass + + def get_obj_value_for_ith_solution(self, i: int) -> float: + pass + + @property + def nb_solutions(self) -> int: + pass + + model: Optional[mathopt.Model] = None + termination: mathopt.Termination + early_stopping_exception: Optional[Exception] = None + + def solve( + self, + callbacks: Optional[List[Callback]] = None, + parameters_milp: Optional[ParametersMilp] = None, + mathopt_solver_type: mathopt.SolverType = mathopt.SolverType.CP_SAT, + time_limit: Optional[float] = 30.0, + mathopt_enable_output: bool = False, + mathopt_model_parameters: Optional[mathopt.ModelSolveParameters] = None, + mathopt_additional_solve_parameters: Optional[mathopt.SolveParameters] = None, + **kwargs: Any, + ) -> ResultStorage: + """Solve with OR-Tools MathOpt API + + Args: + callbacks: list of callbacks used to hook into the various stage of the solve + parameters_milp: parameters passed to the MILP solver + mathopt_solver_type: underlying solver type to use. + Passed as `solver_type` to `mathopt.solve()` + time_limit: the solve process stops after this time limit (in seconds). + If None, no time limit is applied. + mathopt_enable_output: turn on mathopt logging + mathopt_model_parameters: passed to `mathopt.solve()` as `model_params` + mathopt_additional_solve_parameters: passed to `mathopt.solve()` as `params`, + except that parameters defined by above `time_limit` and `parameters_milp` + will be overriden by them. + **kwargs: passed to init_model() if model not yet existing + + Returns: + + """ + self.early_stopping_exception = None + callbacks_list = CallbackList(callbacks=callbacks) + + # callback: solve start + callbacks_list.on_solve_start(solver=self) + + # wrap user callback in a mathopt callback + mathopt_cb = MathOptCallback(do_solver=self, callback=callbacks_list) + + # optimize + self.optimize_model( + parameters_milp=parameters_milp, + time_limit=time_limit, + mathopt_solver_type=mathopt_solver_type, + mathopt_cb=mathopt_cb, + mathopt_enable_output=mathopt_enable_output, + mathopt_model_parameters=mathopt_model_parameters, + mathopt_additional_solve_parameters=mathopt_additional_solve_parameters, + **kwargs, + ) + + # raise potential exception found during callback (useful for optuna pruning, and debugging) + if self.early_stopping_exception: + if isinstance(self.early_stopping_exception, SolveEarlyStop): + logger.info(self.early_stopping_exception) + else: + raise self.early_stopping_exception + + # get result storage + res = mathopt_cb.res + + # callback: solve end + callbacks_list.on_solve_end(res=res, solver=self) + + return res + + def optimize_model( + self, + parameters_milp: Optional[ParametersMilp] = None, + mathopt_solver_type: mathopt.SolverType = mathopt.SolverType.CP_SAT, + time_limit: Optional[float] = 30.0, + mathopt_cb: Optional[mathopt.SolveCallback] = None, + mathopt_enable_output: bool = False, + mathopt_model_parameters: Optional[mathopt.ModelSolveParameters] = None, + mathopt_additional_solve_parameters: Optional[mathopt.SolveParameters] = None, + **kwargs: Any, + ) -> mathopt.SolveResult: + """ + + Args: + parameters_milp: parameters for the milp solver + mathopt_solver_type: underlying solver type to use. + Passed as `solver_type` to `mathopt.solve()` + time_limit: the solve process stops after this time limit (in seconds). + If None, no time limit is applied. + mathopt_cb: a mathopt callback passed to `mathopt.solve()` called at each new solution found + mathopt_enable_output: turn on mathopt logging + mathopt_model_parameters: passed to `mathopt.solve()` as `model_params` + mathopt_additional_solve_parameters: passed to `mathopt.solve()` as `params`, + except that parameters defined by above `time_limit`, `parameters_milp`, and `mathopt_enable_output` + will be overriden by them. + **kwargs: passed to init_model() if model not yet existing + + Returns: + + """ + """Optimize the mip Model. + + The solutions are yet to be retrieved via `self.retrieve_solutions()`. + + Args: + time_limit + + """ + if self.model is None: + self.init_model(**kwargs) + if self.model is None: + raise RuntimeError( + "self.model must not be None after self.init_model()." + ) + + if parameters_milp is None: + parameters_milp = ParametersMilp.default() + + if mathopt_additional_solve_parameters is None: + params = mathopt.SolveParameters() + else: + params = copy.deepcopy(mathopt_additional_solve_parameters) + params.time_limit = datetime.timedelta(seconds=time_limit) + params.absolute_gap_tolerance = parameters_milp.mip_gap_abs + params.relative_gap_tolerance = parameters_milp.mip_gap + params.enable_output = mathopt_enable_output + if mathopt_solver_type != mathopt.SolverType.HIGHS: + # solution_pool_size not supported for HIGHS solver + params.solution_pool_size = parameters_milp.pool_solutions + + callback_reg = mathopt.CallbackRegistration(events={mathopt.Event.MIP_SOLUTION}) + mathopt_res = mathopt.solve( + self.model, + solver_type=mathopt_solver_type, + params=params, + model_params=mathopt_model_parameters, + callback_reg=callback_reg, + cb=mathopt_cb, + ) + self.termination = mathopt_res.termination + logger.info(f"Solver found {len(mathopt_res.solutions)} solutions") + if mathopt_res.termination.reason in [ + mathopt.TerminationReason.OPTIMAL, + mathopt.TerminationReason.FEASIBLE, + ]: + logger.info(f"Objective : {mathopt_res.objective_value()}") + return mathopt_res + + +def _mathopt_cb_get_obj_value_for_current_solution(): + raise RuntimeError("Cannot retrieve objective!") + + +class MathOptCallback: + def __init__(self, do_solver: OrtoolsMathOptMilpSolver, callback: Callback): + self.do_solver = do_solver + self.callback = callback + self.res = do_solver.create_result_storage() + self.nb_solutions = 0 + + def __call__(self, callback_data: mathopt.CallbackData) -> mathopt.CallbackResult: + cb_sol = callback_data.solution + try: + # retrieve and store new solution + sol = self.do_solver.retrieve_current_solution( + get_var_value_for_current_solution=lambda var: cb_sol[var], + get_obj_value_for_current_solution=_mathopt_cb_get_obj_value_for_current_solution, + ) + fit = self.do_solver.aggreg_from_sol(sol) + self.res.append((sol, fit)) + self.nb_solutions += 1 + # end of step callback: stopping? + stopping = self.callback.on_step_end( + step=self.nb_solutions, res=self.res, solver=self.do_solver + ) + except Exception as e: + # catch exceptions because gurobi ignore them and do not stop solving + self.do_solver.early_stopping_exception = e + stopping = True + else: + if stopping: + self.do_solver.early_stopping_exception = SolveEarlyStop( + f"{self.do_solver.__class__.__name__}.solve() stopped by user callback." + ) + return mathopt.CallbackResult(terminate=stopping) + + class GurobiMilpSolver(MilpSolver): """Milp solver wrapping a solver from gurobi library.""" diff --git a/tests/coloring/test_coloring.py b/tests/coloring/test_coloring.py index 9316a4c6f..0d04988cf 100644 --- a/tests/coloring/test_coloring.py +++ b/tests/coloring/test_coloring.py @@ -34,6 +34,7 @@ ColoringCPSatSolver, ModelingCPSat, ) +from discrete_optimization.coloring.solvers.coloring_lp_solvers import ColoringLPMathOpt from discrete_optimization.coloring.solvers.greedy_coloring import ( GreedyColoring, NXGreedyColoringMethod, @@ -564,5 +565,73 @@ def test_color_lp_gurobi_cb_exception(): ) +def test_color_lp_ortools_mathopt(): + file = [f for f in get_data_available() if "gc_70_1" in f][0] + color_problem = parse_file(file) + solver = ColoringLPMathOpt( + color_problem, + params_objective_function=get_default_objective_setup(color_problem), + ) + result_store = solver.solve(parameters_milp=ParametersMilp.default()) + solution = result_store.get_best_solution_fit()[0] + assert color_problem.satisfy(solution) + assert len(result_store) > 1 + + # first solution is not start_solution + assert result_store[0][0].colors != solver.start_solution.colors + + # # warm start => first solution is start_solution + # solver.set_warm_start(solver.start_solution) + # result_store = solver.solve(parameters_milp=ParametersMilp.default()) + # assert result_store[0][0].colors == solver.start_solution.colors + + +def test_color_lp_ortools_mathopt_cb_log(): + file = [f for f in get_data_available() if "gc_70_1" in f][0] + color_problem = parse_file(file) + solver = ColoringLPMathOpt( + color_problem, + params_objective_function=get_default_objective_setup(color_problem), + ) + tracker = NbIterationTracker() + callbacks = [tracker] + result_store = solver.solve( + parameters_milp=ParametersMilp.default(), callbacks=callbacks + ) + solution = result_store.get_best_solution_fit()[0] + assert len(result_store) > 1 + # check tracker called at each solution found + assert tracker.nb_iteration == len(result_store) + + +def test_color_lp_ortools_mathopt_stop(): + file = [f for f in get_data_available() if "gc_70_1" in f][0] + color_problem = parse_file(file) + solver = ColoringLPMathOpt( + color_problem, + params_objective_function=get_default_objective_setup(color_problem), + ) + stopper = NbIterationStopper(nb_iteration_max=1) + callbacks = [stopper] + result_store = solver.solve( + parameters_milp=ParametersMilp.default(), callbacks=callbacks + ) + # check stop after 1st iteration + assert len(result_store) == 1 + + +def test_color_lp_ortools_mathopt_cb_exception(): + file = [f for f in get_data_available() if "gc_70_1" in f][0] + color_problem = parse_file(file) + solver = ColoringLPMathOpt( + color_problem, + params_objective_function=get_default_objective_setup(color_problem), + ) + with pytest.raises(RuntimeError, match="Explicit crash"): + solver.solve( + parameters_milp=ParametersMilp.default(), callbacks=[MyCallbackNok()] + ) + + if __name__ == "__main__": test_solvers()