wip: fix tsp

nucs/examples/tsp/__main__.py

@@ -15,7 +15,8 @@
 from nucs.constants import LOG_LEVEL_INFO, LOG_LEVELS, OPT_MODES, OPT_PRUNE
 from nucs.examples.tsp.tsp_instances import TSP_INSTANCES
 from nucs.examples.tsp.tsp_problem import TSPProblem
-from nucs.heuristics.heuristics import DOM_HEURISTIC_MIN_COST, VAR_HEURISTIC_MAX_REGRET
+from nucs.examples.tsp.tsp_var_heuristic import tsp_var_heuristic
+from nucs.heuristics.heuristics import DOM_HEURISTIC_MIN_COST, register_var_heuristic
 from nucs.solvers.backtrack_solver import BacktrackSolver
 from nucs.solvers.consistency_algorithms import CONSISTENCY_ALG_BC, CONSISTENCY_ALG_SHAVING
 from nucs.solvers.multiprocessing_solver import MultiprocessingSolver
@@ -28,21 +29,25 @@
     parser.add_argument("--name", choices=["GR17", "GR21", "GR24"], default="GR17")
     parser.add_argument("--opt_mode", choices=OPT_MODES, default=OPT_PRUNE)
     parser.add_argument("--processors", type=int, default=1)
-    parser.add_argument("--shaving", type=bool, action=argparse.BooleanOptionalAction, default=True)
+    parser.add_argument("--shaving", type=bool, action=argparse.BooleanOptionalAction, default=False)
     args = parser.parse_args()
     tsp_instance = TSP_INSTANCES[args.name]
+    n = len(tsp_instance)
+    decision_domains = list(range(0, 2 * n))
     problem = TSPProblem(tsp_instance)
+    costs = tsp_instance + tsp_instance  # symmetry of costs
+    tsp_var_heuristic_idx = register_var_heuristic(tsp_var_heuristic)
     solver = (
         MultiprocessingSolver(
             [
                 BacktrackSolver(
                     prob,
                     consistency_alg_idx=CONSISTENCY_ALG_SHAVING if args.shaving else CONSISTENCY_ALG_BC,
-                    decision_domains=list(range(len(tsp_instance))),
-                    var_heuristic_idx=VAR_HEURISTIC_MAX_REGRET,  # register_var_heuristic(mrp_var_heuristic),
-                    var_heuristic_params=tsp_instance,
-                    dom_heuristic_idx=DOM_HEURISTIC_MIN_COST,  # register_dom_heuristic(mcp_dom_heuristic),
-                    dom_heuristic_params=tsp_instance,
+                    decision_domains=decision_domains,
+                    var_heuristic_idx=tsp_var_heuristic_idx,
+                    var_heuristic_params=costs,
+                    dom_heuristic_idx=DOM_HEURISTIC_MIN_COST,
+                    dom_heuristic_params=costs,
                     log_level=args.log_level,
                 )
                 for prob in problem.split(args.processors, 0)
@@ -52,11 +57,11 @@
         else BacktrackSolver(
             problem,
             consistency_alg_idx=CONSISTENCY_ALG_SHAVING if args.shaving else CONSISTENCY_ALG_BC,
-            decision_domains=list(range(len(tsp_instance))),
-            var_heuristic_idx=VAR_HEURISTIC_MAX_REGRET,  # register_var_heuristic(mrp_var_heuristic),
-            var_heuristic_params=tsp_instance,
-            dom_heuristic_idx=DOM_HEURISTIC_MIN_COST,  # register_dom_heuristic(mcp_dom_heuristic),
-            dom_heuristic_params=tsp_instance,
+            decision_domains=decision_domains,
+            var_heuristic_idx=tsp_var_heuristic_idx,
+            var_heuristic_params=costs,
+            dom_heuristic_idx=DOM_HEURISTIC_MIN_COST,
+            dom_heuristic_params=costs,
             log_level=args.log_level,
         )
     )

nucs/examples/tsp/total_cost_propagator.py

@@ -26,7 +26,7 @@ def get_complexity_total_cost(n: int, parameters: NDArray) -> float:
     :param parameters: the parameters, unused here
     :return: a float
     """
-    return 2 * n
+    return n * n
 
 
 def get_triggers_total_cost(n: int, parameters: NDArray) -> NDArray:

nucs/examples/tsp/tsp_problem.py

@@ -33,15 +33,17 @@ def __init__(self, cost_rows: List[List[int]]) -> None:
         super().__init__(n)
         max_costs = [max(cost_row) for cost_row in cost_rows]
         min_costs = [min([cost for cost in cost_row if cost > 0]) for cost_row in cost_rows]
-        cost_vars = self.add_variables([(min_costs[i], max_costs[i]) for i in range(n)])  # the costs
-        total_cost_var = self.add_variable((sum(min_costs), sum(max_costs)))  # the total cost
+        next_costs = self.add_variables([(min_costs[i], max_costs[i]) for i in range(n)])
+        prev_costs = self.add_variables([(min_costs[i], max_costs[i]) for i in range(n)])
+        total_cost = self.add_variable((sum(min_costs), sum(max_costs)))  # the total cost
         for i in range(n):
-            self.add_propagator(([i, cost_vars + i], ALG_ELEMENT_IV, cost_rows[i]))
-        self.add_propagator((list(range(cost_vars, cost_vars + n + 1)), ALG_AFFINE_EQ, [1] * n + [-1, 0]))
-        self.add_propagator(
-            (
-                list(range(0, n)) + [total_cost_var],
-                register_propagator(get_triggers_total_cost, get_complexity_total_cost, compute_domains_total_cost),
-                [cost for cost_row in cost_rows for cost in cost_row],
-            )
+            self.add_propagator(([i, next_costs + i], ALG_ELEMENT_IV, cost_rows[i]))
+            self.add_propagator(([n + i, prev_costs + i], ALG_ELEMENT_IV, cost_rows[i]))
+        self.add_propagator((list(range(next_costs, next_costs + n)) + [total_cost], ALG_AFFINE_EQ, [1] * n + [-1, 0]))
+        self.add_propagator((list(range(prev_costs, prev_costs + n)) + [total_cost], ALG_AFFINE_EQ, [1] * n + [-1, 0]))
+        total_cost_prop_idx = register_propagator(
+            get_triggers_total_cost, get_complexity_total_cost, compute_domains_total_cost
         )
+        costs = [cost for cost_row in cost_rows for cost in cost_row]
+        self.add_propagator((list(range(n)) + [total_cost], total_cost_prop_idx, costs))
+        self.add_propagator((list(range(n, 2 * n)) + [total_cost], total_cost_prop_idx, costs))

nucs/examples/tsp/tsp_var_heuristic.py

@@ -0,0 +1,52 @@
+###############################################################################
+# __   _            _____    _____
+# | \ | |          / ____|  / ____|
+# |  \| |  _   _  | |      | (___
+# | . ` | | | | | | |       \___ \
+# | |\  | | |_| | | |____   ____) |
+# |_| \_|  \__,_|  \_____| |_____/
+#
+# Fast constraint solving in Python  - https://github.com/yangeorget/nucs
+#
+# Copyright 2024 - Yan Georget
+###############################################################################
+import sys
+
+from numba import njit  # type: ignore
+from numpy.typing import NDArray
+
+from nucs.constants import MAX, MIN
+from nucs.heuristics.max_regret_var_heuristic import regret
+
+
+@njit(cache=True)
+def tsp_var_heuristic(
+    decision_domains: NDArray, shr_domains_stack: NDArray, stacks_top: NDArray, params: NDArray
+) -> int:
+    """
+    :param decision_domains: the indices of a subset of the shared domains
+    :param shr_domains_stack: the stack of shared domains
+    :param stacks_top: the index of the top of the stacks as a Numpy array
+    :param params: a two-dimensional (first dimension correspond to variables, second to values) costs array
+    :return: the index of the shared domain
+    """
+    best_score = -sys.maxsize
+    best_idx = -1
+    top = stacks_top[0]
+    for dom_idx in decision_domains:
+        shr_domain = shr_domains_stack[top, dom_idx]
+        if 0 < shr_domain[MAX] - shr_domain[MIN]:
+            score = compute_score(shr_domain, dom_idx, params)
+            if best_score < score:
+                best_idx = dom_idx
+                best_score = score
+    return best_idx
+
+
+@njit(cache=True)
+def compute_score(shr_domain: NDArray, dom_idx: int, params: NDArray) -> int:
+    """
+    Minimize [min(5, size(X)), -regret(X)] for lexicographic order.
+    """
+    size = min(5, 1 + shr_domain[MAX] - shr_domain[MIN])
+    return -size * 1024 + regret(shr_domain, dom_idx, params)

nucs/heuristics/greatest_domain_var_heuristic.py

@@ -28,13 +28,13 @@ def greatest_domain_var_heuristic(
     :param params: a two-dimensional parameters array, unused here
     :return: the index of the shared domain
     """
-    max_size = 0
-    max_idx = -1
+    best_score = 0
+    best_idx = -1
     top = stacks_top[0]
     for dom_idx in decision_domains:
         shr_domain = shr_domains_stack[top, dom_idx]
-        size = shr_domain[MAX] - shr_domain[MIN]  # actually this is size - 1
-        if max_size < size:
-            max_idx = dom_idx
-            max_size = size
-    return max_idx
+        score = shr_domain[MAX] - shr_domain[MIN]  # this is size - 1
+        if best_score < score:
+            best_idx = dom_idx
+            best_score = score
+    return best_idx

nucs/heuristics/max_regret_var_heuristic.py

@@ -30,24 +30,29 @@ def max_regret_var_heuristic(
     :param params: a two-dimensional (first dimension correspond to variables, second to values) costs array
     :return: the index of the shared domain
     """
-    max_regret = 0
+    best_score = 0
     best_idx = -1
     top = stacks_top[0]
     for dom_idx in decision_domains:
         shr_domain = shr_domains_stack[top, dom_idx]
         if 0 < shr_domain[MAX] - shr_domain[MIN]:
-            best_cost = sys.maxsize
-            second_cost = sys.maxsize
-            for value in range(shr_domain[MIN], shr_domain[MAX] + 1):
-                cost = params[dom_idx][value]
-                if cost > 0:
-                    if cost < best_cost:
-                        second_cost = best_cost
-                        best_cost = cost
-                    elif cost < second_cost:
-                        second_cost = cost
-            regret = second_cost - best_cost
-            if max_regret < regret:
+            score = regret(shr_domain, dom_idx, params)
+            if best_score < score:
                 best_idx = dom_idx
-                max_regret = regret
+                best_score = score
     return best_idx
+
+
+@njit(cache=True)
+def regret(shr_domain: NDArray, dom_idx: int, params: NDArray) -> int:
+    best_cost = sys.maxsize
+    second_cost = sys.maxsize
+    for value in range(shr_domain[MIN], shr_domain[MAX] + 1):
+        cost = params[dom_idx][value]
+        if cost > 0:
+            if cost < best_cost:
+                second_cost = best_cost
+                best_cost = cost
+            elif cost < second_cost:
+                second_cost = cost
+    return second_cost - best_cost

nucs/heuristics/smallest_domain_var_heuristic.py

@@ -30,13 +30,13 @@ def smallest_domain_var_heuristic(
     :param params: a two-dimensional parameters array, unused here
     :return: the index of the shared domain
     """
-    min_size = sys.maxsize
-    min_idx = -1
+    best_score = -sys.maxsize
+    best_idx = -1
     top = stacks_top[0]
     for dom_idx in decision_domains:
         shr_domain = shr_domains_stack[top, dom_idx]
-        size = shr_domain[MAX] - shr_domain[MIN]  # actually this is size - 1
-        if 0 < size < min_size:
-            min_idx = dom_idx
-            min_size = size
-    return min_idx
+        score = shr_domain[MIN] - shr_domain[MAX]
+        if best_score < score < 0:
+            best_idx = dom_idx
+            best_score = score
+    return best_idx

nucs/problems/circuit_problem.py

@@ -10,11 +10,11 @@
 #
 # Copyright 2024 - Yan Georget
 ###############################################################################
-from nucs.problems.problem import Problem
-from nucs.propagators.propagators import ALG_ALLDIFFERENT, ALG_NO_SUB_CYCLE
+from nucs.problems.permutation_problem import PermutationProblem
+from nucs.propagators.propagators import ALG_NO_SUB_CYCLE
 
 
-class CircuitProblem(Problem):
+class CircuitProblem(PermutationProblem):
     """
     A model for circuit.
     """
@@ -25,9 +25,8 @@ def __init__(self, n: int):
         :param n: the number of vertices
         """
         self.n = n
-        shr_domains = [(1, n - 1)] + [(0, n - 1)] * (n - 2) + [(0, n - 2)]
-        super().__init__(shr_domains)
-        s_indices = list(range(n))
-        self.add_propagator((s_indices, ALG_ALLDIFFERENT, []))
-        self.add_propagator((s_indices, ALG_NO_SUB_CYCLE, []))
-        # self.add_propagator((s_indices, ALG_SCC, []))  # not worth the cost
+        super().__init__(n)  # TODO: do we want to specify the domains more precisely ?
+        self.add_propagator((list(range(n)), ALG_NO_SUB_CYCLE, []))
+        self.add_propagator((list(range(n, 2 * n)), ALG_NO_SUB_CYCLE, []))
+        # self.add_propagator((list(range(n)), ALG_SCC, []))  # not worth the cost
+        # self.add_propagator((list(range(n, 2*n)), ALG_SCC, []))  # not worth the cost

nucs/problems/permutation_problem.py

@@ -0,0 +1,34 @@
+###############################################################################
+# __   _            _____    _____
+# | \ | |          / ____|  / ____|
+# |  \| |  _   _  | |      | (___
+# | . ` | | | | | | |       \___ \
+# | |\  | | |_| | | |____   ____) |
+# |_| \_|  \__,_|  \_____| |_____/
+#
+# Fast constraint solving in Python  - https://github.com/yangeorget/nucs
+#
+# Copyright 2024 - Yan Georget
+###############################################################################
+from nucs.problems.problem import Problem
+from nucs.propagators.propagators import ALG_ALLDIFFERENT, ALG_PERMUTATION_AUX
+
+
+class PermutationProblem(Problem):
+    """
+    A model for permutation.
+    """
+
+    def __init__(self, n: int):
+        """
+        Inits the permutation problem.
+        :param n: the number variables/values
+        """
+        self.n = n
+        shr_domains = [(0, n - 1)] * 2 * n
+        super().__init__(shr_domains)
+        for i in range(n):
+            self.add_propagator((list(range(n)) + [n + i], ALG_PERMUTATION_AUX, [i]))
+            self.add_propagator((list(range(n, 2 * n)) + [i], ALG_PERMUTATION_AUX, [i]))
+        self.add_propagator((list(range(n)), ALG_ALLDIFFERENT, []))
+        self.add_propagator((list(range(n, 2 * n)), ALG_ALLDIFFERENT, []))

nucs/propagators/no_sub_cycle_propagator.py

@@ -56,23 +56,26 @@ def compute_domains_no_sub_cycle(domains: NDArray, parameters: NDArray) -> int:
     while loop:
         loop = False
         for i in range(n):
-            if domains[i, MIN] == domains[i, MAX] and paths[i, PATH_END] == i:
+            if domains[i, MIN] == domains[i, MAX]:
                 j = domains[i, MIN]
-                end = paths[i, PATH_END] = paths[j, PATH_END]
-                start = paths[j, PATH_START] = paths[i, PATH_START]
-                paths[start, PATH_END] = end
-                paths[end, PATH_START] = start
-                length = paths[i, PATH_LENGTH] + 1 + paths[j, PATH_LENGTH]
-                paths[i, PATH_LENGTH] = paths[j, PATH_LENGTH] = paths[start, PATH_LENGTH] = paths[end, PATH_LENGTH] = (
-                    length
-                )
-                if length < n - 1:
-                    if domains[end, MIN] == start:
-                        domains[end, MIN] = start + 1
-                    if domains[end, MAX] == start:
-                        domains[end, MAX] = start - 1
-                    if domains[end, MIN] > domains[end, MAX]:
-                        return PROP_INCONSISTENCY
-                    if end < i:
-                        loop = True
+                if i == j:
+                    return PROP_INCONSISTENCY
+                if paths[i, PATH_END] == i:
+                    end = paths[i, PATH_END] = paths[j, PATH_END]
+                    start = paths[j, PATH_START] = paths[i, PATH_START]
+                    paths[start, PATH_END] = end
+                    paths[end, PATH_START] = start
+                    length = paths[i, PATH_LENGTH] + 1 + paths[j, PATH_LENGTH]
+                    paths[i, PATH_LENGTH] = paths[j, PATH_LENGTH] = paths[start, PATH_LENGTH] = paths[
+                        end, PATH_LENGTH
+                    ] = length
+                    if length < n - 1:
+                        if domains[end, MIN] == start:
+                            domains[end, MIN] = start + 1
+                        if domains[end, MAX] == start:
+                            domains[end, MAX] = start - 1
+                        if domains[end, MIN] > domains[end, MAX]:
+                            return PROP_INCONSISTENCY
+                        if end < i:
+                            loop = True
     return PROP_CONSISTENCY