Feat: solve TSP using branch and bound method

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,150 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from math import inf
+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
+        The level of the node in the search tree.
+    index
+        The index of the current city in the path.
+    path
+        The list of city indices visited so far.
+    cost
+        The total cost of the path up to this node.
+    cost_matrix
+        The cost matrix representing the distances between cities.
+
+    Methods
+    -------
+    compute_reduced_matrix
+        Compute the reduced matrix and the cost of reducing it.
+    from_cost_matrix
+        Create a Node object from a given cost matrix.
+    from_parent
+        Create a new Node object based on a parent node and a city index.
+    """
+
+    level: int
+    index: int
+    path: List[int]
+    cost: float
+    cost_matrix: np.ndarray
+
+    @staticmethod
+    def compute_reduced_matrix(matrix: np.ndarray) -> Tuple[np.ndarray, float]:
+        """
+        Compute the reduced matrix and the cost of reducing it.
+
+        Parameters
+        ----------
+        matrix
+            The cost matrix to compute the reductions.
+
+        Returns
+        -------
+        Tuple
+            A tuple containing the reduced matrix and the total
+            cost of reductions.
+        """
+        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, np.sum(min_rows) + np.sum(min_cols)
+
+    @classmethod
+    def from_cost_matrix(cls, cost_matrix: np.ndarray) -> Node:
+        """
+        Create a Node object from a given cost matrix.
+
+        Parameters
+        ----------
+        cost_matrix
+            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
+            The parent node.
+        index
+            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)
+        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
+            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
+        The list that holds the elements in the priority queue.
+
+    Methods
+    -------
+    empty
+        Check if the priority queue is empty.
+    push
+        Push an item into the priority queue.
+    pop
+        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
+            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,68 @@
+from math import inf
+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
+        The distance matrix representing the distances between cities.
+
+    Returns
+    -------
+    Tuple
+        A tuple containing the optimal path (list of city indices) and its
+        total cost. If the TSP cannot be solved, an empty path and a cost
+        of positive infinity will be returned.
+
+    Notes
+    -----
+    The `distance_matrix` should be a square matrix with non-negative
+    values. The element `distance_matrix[i][j]` represents the distance from
+    city `i` to city `j`. If two cities are not directly connected, the
+    distance should be set to a float value of positive infinity
+    (float('inf')).
+
+    The path is represented as a list of city indices, and the total cost is a
+    float value indicating the sum of distances in the optimal path.
+
+    If the TSP cannot be solved (e.g., due to disconnected cities), the
+    function will return an empty path ([]) and a cost of positive infinity
+    (float('inf')).
+
+    References
+    ----------
+    .. [1] Horowitz, E., Sahni, S., & Rajasekaran, S. (1997).
+           Computer Algorithms. Chapter 8 - Branch and Bound. Section 8.3.
+           W. H. Freeman and Company.
+    """
+    num_cities = len(distance_matrix)
+    cost_matrix = np.copy(distance_matrix).astype(float)
+    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 [], inf