Feat: solve TSP using branch and bound method

.gitignore

@@ -142,3 +142,6 @@ dmypy.json
 
 # Cython debug symbols
 cython_debug/
+
+# Pycharm IDE
+.idea/

README.rst

@@ -196,7 +196,8 @@ There are two types of solvers available:
         - ``exact.solve_tsp_dynamic_programming``: uses a Dynamic Programming
           approach. It tends to be faster than the previous one, but it may
           demand more memory.
-
+        - ``exact.solve_tsp_branch_and_bound``: uses a Branch and Bound
+          approach, it is more scalable than previous methods.
 :Heuristics: These methods have no guarantees of finding the best solution,
              but usually return a good enough candidate in a more reasonable
              time for larger problems.

python_tsp/exact/__init__.py

@@ -14,19 +14,36 @@
 
 Functions present in this module are listed below.
 
+    solve_tsp_branch_and_bound
     solve_tsp_brute_force
     solve_tsp_dynamic_programming
 
 
 Tips on deciding the solver
 ---------------------------
 
-In general, ``solve_tsp_dynamic_programming`` is faster, specially with an
-appropriate `maxsize` input (the default is fine). However, because of its
-recursion, it may take more memory, particularly if the number of nodes grows
-large. If that becomes an issue and you still need a provably optimal solution,
-use the ``solve_tsp_brute_force``.
+The choice among the three exact methods depends on the specific
+characteristics of the Traveling Salesperson Problem (TSP) you are
+dealing with:
+
+If the TSP has only a few cities and the goal is a quick solution without
+worrying about scalability, ``solve_tsp_brute_force`` may be a simple
+and viable choice, but only for educational purposes or small cases.
+If the TSP is relatively small (with a few cities) and precision is
+essential, ``solve_tsp_dynamic_programming`` may be preferable, as
+long as the required memory and execution time are not prohibitive.
+If the TSP has many cities and an exact solution is required,
+``solve_tsp_branch_and_bound`` is more scalable and, therefore, more
+suitable for such scenarios.
+
+In general, ``solve_tsp_brute_force`` is not recommended for TSPs of
+significant size due to its exponential complexity.
+``solve_tsp_dynamic_programming`` and ``solve_tsp_branch_and_bound``
+are more efficient approaches to finding the optimal solution, but the
+choice between them will depend on the problem size and available
+computational resources.
 """
 
+from .branch_and_bound import solve_tsp_branch_and_bound  # noqa: F401
 from .brute_force import solve_tsp_brute_force  # noqa: F401
 from .dynamic_programming import solve_tsp_dynamic_programming  # noqa: F401

python_tsp/exact/branch_and_bound/__init__.py

@@ -0,0 +1,3 @@
+from .node import Node  # noqa: F401
+from .priority_queue import PriorityQueue  # noqa: F401
+from .solver import solve_tsp_branch_and_bound  # noqa: F401

python_tsp/exact/branch_and_bound/node.py

@@ -0,0 +1,151 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import List, Tuple
+
+import numpy as np
+
+
+@dataclass
+class Node:
+    """
+    Represents a node in the search tree for the Traveling Salesperson Problem.
+
+    Attributes
+    ----------
+    level : int
+        The level of the node in the search tree.
+    index : int
+        The index of the current city in the path.
+    path : List[int]
+        The list of city indices visited so far.
+    cost : int
+        The total cost of the path up to this node.
+    cost_matrix : numpy.ndarray
+        The cost matrix representing the distances between cities.
+
+    Methods
+    -------
+    compute_reduced_matrix(matrix: numpy.ndarray) -> Tuple[numpy.ndarray, int]:
+        Compute the reduced matrix and the cost of reducing it.
+    from_cost_matrix(cost_matrix: numpy.ndarray) -> Node:
+        Create a Node object from a given cost matrix.
+    from_parent(parent: Node, index: int) -> Node:
+        Create a new Node object based on a parent node and a city index.
+    """
+
+    level: int
+    index: int
+    path: List[int]
+    cost: int
+    cost_matrix: np.ndarray
+
+    @staticmethod
+    def compute_reduced_matrix(matrix: np.ndarray) -> Tuple[np.ndarray, int]:
+        """
+        Compute the reduced matrix and the cost of reducing it.
+
+        Parameters
+        ----------
+        matrix : numpy.ndarray
+            The cost matrix to compute the reductions.
+
+        Returns
+        -------
+        Tuple[numpy.ndarray, int]
+            A tuple containing the reduced matrix and the total
+            cost of reductions.
+        """
+        inf = np.iinfo(matrix.dtype).max
+        mask = matrix != inf
+        reduced_matrix = np.copy(matrix)
+
+        min_rows = np.min(reduced_matrix, axis=1, keepdims=True)
+        min_rows[min_rows == inf] = 0
+        if np.any(min_rows != 0):
+            reduced_matrix = np.where(
+                mask, reduced_matrix - min_rows, reduced_matrix
+            )
+
+        min_cols = np.min(reduced_matrix, axis=0, keepdims=True)
+        min_cols[min_cols == inf] = 0
+        if np.any(min_cols != 0):
+            reduced_matrix = np.where(
+                mask, reduced_matrix - min_cols, reduced_matrix
+            )
+
+        return reduced_matrix, min_rows.sum() + min_cols.sum()
+
+    @classmethod
+    def from_cost_matrix(cls, cost_matrix: np.ndarray) -> Node:
+        """
+        Create a Node object from a given cost matrix.
+
+        Parameters
+        ----------
+        cost_matrix : numpy.ndarray
+            The cost matrix representing the distances between cities.
+
+        Returns
+        -------
+        Node
+            A new Node object initialized with the reduced cost matrix.
+        """
+        _cost_matrix, _cost = cls.compute_reduced_matrix(matrix=cost_matrix)
+        return cls(
+            level=0,
+            index=0,
+            path=[0],
+            cost=_cost,
+            cost_matrix=_cost_matrix,
+        )
+
+    @classmethod
+    def from_parent(cls, parent: Node, index: int) -> Node:
+        """
+        Create a new Node object based on a parent node and a city index.
+
+        Parameters
+        ----------
+        parent : Node
+            The parent node.
+        index : int
+            The index of the new city to be added to the path.
+
+        Returns
+        -------
+        Node
+            A new Node object with the updated path and cost.
+        """
+        matrix = np.copy(parent.cost_matrix)
+        inf = np.iinfo(matrix.dtype).max
+        matrix[parent.index, :] = inf
+        matrix[:, index] = inf
+        matrix[index][0] = inf
+        _cost_matrix, _cost = cls.compute_reduced_matrix(matrix=matrix)
+        return cls(
+            level=parent.level + 1,
+            index=index,
+            path=parent.path[:] + [index],
+            cost=(
+                parent.cost + _cost + parent.cost_matrix[parent.index][index]
+            ),
+            cost_matrix=_cost_matrix,
+        )
+
+    def __lt__(self: Node, other: Node):
+        """
+        Compare two Node objects based on their costs.
+
+        Parameters
+        ----------
+        other : Node
+            The other Node object to compare with.
+
+        Returns
+        -------
+        bool
+            True if this Node's cost is less than the other Node's
+            cost, False otherwise.
+        """
+        return self.cost < other.cost

python_tsp/exact/branch_and_bound/priority_queue.py

@@ -0,0 +1,67 @@
+from dataclasses import dataclass, field
+from heapq import heappop, heappush
+from typing import List
+
+from python_tsp.exact.branch_and_bound import Node
+
+
+@dataclass
+class PriorityQueue:
+    """
+    A priority queue implementation using a binary heap
+    for efficient element retrieval.
+
+    Attributes
+    ----------
+    _container : List[Node]
+        The list that holds the elements in the priority queue.
+
+    Methods
+    -------
+    empty() -> bool:
+        Check if the priority queue is empty.
+    push(item: Node) -> None:
+        Push an item into the priority queue.
+    pop() -> Node:
+        Pop the item with the highest priority from the priority queue.
+    """
+
+    _container: List[Node] = field(default_factory=list)
+
+    @property
+    def empty(self) -> bool:
+        """
+        Check if the priority queue is empty.
+
+        Returns
+        -------
+        bool
+            True if the priority queue is empty, False otherwise.
+        """
+        return not self._container
+
+    def push(self, item: Node) -> None:
+        """
+        Push an item into the priority queue.
+
+        Parameters
+        ----------
+        item : Node
+            The item to be pushed into the priority queue.
+
+        Returns
+        -------
+        None
+        """
+        heappush(self._container, item)
+
+    def pop(self) -> Node:
+        """
+        Pop the item with the highest priority from the priority queue.
+
+        Returns
+        -------
+        Node
+            The node with the highest priority.
+        """
+        return heappop(self._container)

python_tsp/exact/branch_and_bound/solver.py

@@ -0,0 +1,47 @@
+from typing import List, Tuple
+
+import numpy as np
+
+from python_tsp.exact.branch_and_bound import Node, PriorityQueue
+
+
+def solve_tsp_branch_and_bound(
+    distance_matrix: np.ndarray,
+) -> Tuple[List[int], float]:
+    """
+    Solve the Traveling Salesperson Problem (TSP) using
+    the Branch and Bound algorithm.
+
+    Parameters
+    ----------
+    distance_matrix : numpy.ndarray
+        The distance matrix representing the distances between cities.
+
+    Returns
+    -------
+    Tuple[List[int], float]
+        A tuple containing the optimal path (list of city indices) and its
+        total cost. If the TSP cannot be solved, an empty path and cost of 0
+        will be returned.
+    """
+    num_cities = len(distance_matrix)
+    cost_matrix = np.copy(distance_matrix)
+    inf = np.iinfo(cost_matrix.dtype).max
+    np.fill_diagonal(cost_matrix, inf)
+
+    root = Node.from_cost_matrix(cost_matrix=cost_matrix)
+    pq = PriorityQueue([root])
+
+    while not pq.empty:
+        min_node = pq.pop()
+
+        if min_node.level == num_cities - 1:
+            return min_node.path, min_node.cost
+
+        for index in range(num_cities):
+            is_live_node = min_node.cost_matrix[min_node.index][index] != inf
+            if is_live_node:
+                live_node = Node.from_parent(parent=min_node, index=index)
+                pq.push(live_node)
+
+    return [], 0