Introduce a milp solver base class using ortools/mathopt api

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,
 )
 
@@ -110,6 +112,31 @@ def __init__(
         self.sense_optim = self.params_objective_function.sense_function
         self.start_solution: Optional[ColoringSolution] = None
 
+    def convert_to_variable_values(
+        self, solution: ColoringSolution
+    ) -> dict[Any, float]:
+        """Convert a solution to a mapping between model variables and their values.
+
+        Will be used by set_warm_start().
+
+        """
+        # Init all variables to 0
+        hinted_variables = {
+            var: 0 for var in self.variable_decision["colors_var"].values()
+        }
+
+        # Set var(node, color) to 1 according to the solution
+        for i, color in enumerate(solution.colors):
+            node = self.index_to_nodes_name[i]
+            variable_decision_key = (node, color)
+            hinted_variables[
+                self.variable_decision["colors_var"][variable_decision_key]
+            ] = 1
+
+        hinted_variables[self.variable_decision["nb_colors"]] = solution.nb_color
+
+        return hinted_variables
+
     def retrieve_current_solution(
         self,
         get_var_value_for_current_solution: Callable[[Any], float],
@@ -478,3 +505,162 @@ 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
+    solution_hint: Optional[dict[mathopt.Variable, float]] = None
+
+    def convert_to_variable_values(
+        self, solution: ColoringSolution
+    ) -> dict[mathopt.Variable, float]:
+        """Convert a solution to a mapping between model variables and their values.
+
+        Will be used by set_warm_start() to provide a suitable SolutionHint.variable_values.
+        See https://or-tools.github.io/docs/pdoc/ortools/math_opt/python/model_parameters.html#SolutionHint
+        for more information.
+
+        """
+        return _BaseColoringLP.convert_to_variable_values(self, solution)
+
+    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, "nb_colors": opt}
+        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"
+        }

discrete_optimization/facility/solvers/facility_lp_solver.py

@@ -12,6 +12,7 @@
 import numpy as np
 import numpy.typing as npt
 from ortools.linear_solver import pywraplp
+from ortools.math_opt.python import mathopt
 
 from discrete_optimization.facility.facility_model import (
     FacilityProblem,
@@ -31,6 +32,7 @@
     GurobiMilpSolver,
     MilpSolver,
     MilpSolverName,
+    OrtoolsMathOptMilpSolver,
     ParametersMilp,
     PymipMilpSolver,
 )
@@ -151,7 +153,9 @@ def __init__(
             problem=problem, params_objective_function=params_objective_function
         )
         self.model = None
-        self.variable_decision: dict[str, dict[tuple[int, int], Union[int, Any]]] = {}
+        self.variable_decision: dict[
+            str, dict[Union[int, tuple[int, int]], Union[int, Any]]
+        ] = {}
         self.constraints_dict: dict[str, dict[int, Any]] = {}
         self.description_variable_description = {
             "x": {
@@ -164,6 +168,24 @@ def __init__(
         }
         self.description_constraint: dict[str, dict[str, str]] = {}
 
+    def convert_to_variable_values(
+        self, solution: FacilitySolution
+    ) -> dict[Any, float]:
+        """Convert a solution to a mapping between model variables and their values.
+
+        Will be used by set_warm_start().
+
+        """
+        # Init all variables to 0
+        hinted_variables = {var: 0 for var in self.variable_decision["x"].values()}
+        # Set var(facility, customer) to 1 according to the solution
+        for c, f in enumerate(solution.facility_for_customers):
+            variable_decision_key = (f, c)
+            hinted_variables[self.variable_decision["x"][variable_decision_key]] = 1
+            hinted_variables[self.variable_decision["y"][f]] = 1
+
+        return hinted_variables
+
     def retrieve_current_solution(
         self,
         get_var_value_for_current_solution: Callable[[Any], float],
@@ -299,6 +321,111 @@ def set_warm_start(self, solution: FacilitySolution) -> None:
             self.variable_decision["x"][variable_decision_key].Start = 1
 
 
+class LP_Facility_Solver_MathOpt(OrtoolsMathOptMilpSolver, _LPFacilitySolverBase):
+    """Milp solver using gurobi library
+
+    Attributes:
+        coloring_model (FacilityProblem): facility problem instance to solve
+        params_objective_function (ParamsObjectiveFunction): objective function parameters
+                        (however this is just used for the ResultStorage creation, not in the optimisation)
+    """
+
+    def init_model(self, **kwargs: Any) -> None:
+        """
+
+        Keyword Args:
+            use_matrix_indicator_heuristic (bool): use the prune search method to reduce number of variable.
+            n_shortest (int): parameter for the prune search method
+            n_cheapest (int): parameter for the prune search method
+
+        Returns: None
+
+        """
+        nb_facilities = self.problem.facility_count
+        nb_customers = self.problem.customer_count
+        kwargs = self.complete_with_default_hyperparameters(kwargs)
+        use_matrix_indicator_heuristic = kwargs["use_matrix_indicator_heuristic"]
+        if use_matrix_indicator_heuristic:
+            n_shortest = kwargs["n_shortest"]
+            n_cheapest = kwargs["n_cheapest"]
+            matrix_fc_indicator, matrix_length = prune_search_space(
+                self.problem, n_cheapest=n_cheapest, n_shortest=n_shortest
+            )
+        else:
+            matrix_fc_indicator, matrix_length = prune_search_space(
+                self.problem,
+                n_cheapest=nb_facilities,
+                n_shortest=nb_facilities,
+            )
+        s = mathopt.Model(name="facilities")
+        x: dict[tuple[int, int], Union[int, Any]] = {}
+        for f in range(nb_facilities):
+            for c in range(nb_customers):
+                if matrix_fc_indicator[f, c] == 0:
+                    x[f, c] = 0
+                elif matrix_fc_indicator[f, c] == 1:
+                    x[f, c] = 1
+                elif matrix_fc_indicator[f, c] == 2:
+                    x[f, c] = s.add_binary_variable(name="x_" + str((f, c)))
+        facilities = self.problem.facilities
+        customers = self.problem.customers
+        used = [s.add_binary_variable(name=f"y_{i}") for i in range(nb_facilities)]
+        constraints_customer: dict[int, mathopt.LinearConstraint] = {}
+        for c in range(nb_customers):
+            constraints_customer[c] = s.add_linear_constraint(
+                mathopt.LinearSum(x[f, c] for f in range(nb_facilities)) == 1
+            )
+            # one facility
+        constraint_capacity: dict[int, mathopt.LinearConstraint] = {}
+        for f in range(nb_facilities):
+            for c in range(nb_customers):
+                s.add_linear_constraint(used[f] >= x[f, c])
+            constraint_capacity[f] = s.add_linear_constraint(
+                mathopt.LinearSum(
+                    x[f, c] * customers[c].demand for c in range(nb_customers)
+                )
+                <= facilities[f].capacity
+            )
+        new_obj_f = 0.0
+        new_obj_f += mathopt.LinearSum(
+            facilities[f].setup_cost * used[f] for f in range(nb_facilities)
+        )
+        new_obj_f += mathopt.LinearSum(
+            matrix_length[f, c] * x[f, c]
+            for f in range(nb_facilities)
+            for c in range(nb_customers)
+        )
+        s.minimize(new_obj_f)
+        self.model = s
+        self.variable_decision = {"x": x, "y": used}
+        self.constraints_dict = {
+            "constraint_customer": constraints_customer,
+            "constraint_capacity": constraint_capacity,
+        }
+        self.description_variable_description = {
+            "x": {
+                "shape": (nb_facilities, nb_customers),
+                "type": bool,
+                "descr": "for each facility/customer indicate"
+                " if the pair is active, meaning "
+                "that the customer c is dealt with facility f",
+            }
+        }
+        logger.info("Initialized")
+
+    def convert_to_variable_values(
+        self, solution: FacilitySolution
+    ) -> dict[mathopt.Variable, float]:
+        """Convert a solution to a mapping between model variables and their values.
+
+        Will be used by set_warm_start() to provide a suitable SolutionHint.variable_values.
+        See https://or-tools.github.io/docs/pdoc/ortools/math_opt/python/model_parameters.html#SolutionHint
+        for more information.
+
+        """
+        return _LPFacilitySolverBase.convert_to_variable_values(self, solution)
+
+
 class LP_Facility_Solver_CBC(SolverFacility):
     """Milp formulation using cbc solver."""