graph_advanced/
tsp.rs

1/// Solves the **Travelling Salesman Problem** exactly using the Held-Karp
2/// bitmask dynamic programming algorithm.
3///
4/// Given a complete weighted directed graph represented as a distance matrix,
5/// finds the minimum-cost Hamiltonian circuit starting and ending at node `0`.
6///
7/// # Algorithm
8///
9/// Let `dp[mask][v]` be the minimum cost to reach node `v` having visited
10/// exactly the nodes encoded in `mask` (bit `i` set ⟹ node `i` visited),
11/// starting from node `0`.
12///
13/// **Recurrence:**  
14/// `dp[mask | (1 << u)][u] = min(dp[mask][v] + dist[v][u])`  
15/// for all `v` in `mask` and `u` not in `mask`.
16///
17/// **Answer:**  
18/// `min over v of (dp[all_visited][v] + dist[v][0])`
19///
20/// # Arguments
21///
22/// - `dist` — `n × n` distance matrix. `dist[i][j]` is the cost of going
23///   from node `i` to node `j`. Use `f64::INFINITY` for missing edges.
24///
25/// # Returns
26///
27/// `Some((cost, path))` — the minimum tour cost and the node visit order
28/// (starts and ends at node `0`).  
29/// `None` — no Hamiltonian circuit exists (some required edge is `INFINITY`),
30/// or `dist` is empty.
31///
32/// # Complexity
33///
34/// O(2^n · n²) time, O(2^n · n) space. Practical for n ≤ 20.
35///
36/// # Panics
37///
38/// Panics if `dist` is not square (all rows must have the same length as the
39/// number of rows).
40///
41/// # Examples
42///
43/// ```
44/// use graph_advanced::tsp_held_karp;
45///
46/// // 4-node complete graph with known optimal tour cost 35.
47/// //        0   1   2   3
48/// let dist = vec![
49///     vec![0.0, 10.0, 15.0, 20.0], // from 0
50///     vec![10.0, 0.0, 35.0, 25.0], // from 1
51///     vec![15.0, 35.0, 0.0, 30.0], // from 2
52///     vec![20.0, 25.0, 30.0, 0.0], // from 3
53/// ];
54///
55/// let (cost, path) = tsp_held_karp(&dist).unwrap();
56/// assert_eq!(cost, 80.0); // 0→1(10)+1→3(25)+3→2(30)+2→0(15)
57/// assert_eq!(path.first(), Some(&0));
58/// assert_eq!(path.last(), Some(&0));
59/// assert_eq!(path.len(), 5); // 4 nodes + return to 0
60/// ```
61pub fn tsp_held_karp(dist: &[Vec<f64>]) -> Option<(f64, Vec<usize>)> {
62    let n = dist.len();
63    if n == 0 {
64        return None;
65    }
66    assert!(
67        dist.iter().all(|row| row.len() == n),
68        "distance matrix must be square"
69    );
70    if n == 1 {
71        return Some((0.0, vec![0, 0]));
72    }
73
74    let full_mask = (1usize << n) - 1;
75
76    // dp[mask][v] = min cost to reach v having visited exactly nodes in mask.
77    // parent[mask][v] = the node we came from to achieve dp[mask][v].
78    let mut dp = vec![vec![f64::INFINITY; n]; 1 << n];
79    let mut parent: Vec<Vec<usize>> = vec![vec![usize::MAX; n]; 1 << n];
80
81    // Start at node 0.
82    dp[1][0] = 0.0;
83
84    for mask in 1..=(full_mask) {
85        // Only process masks that include node 0 (bit 0 always set).
86        if mask & 1 == 0 {
87            continue;
88        }
89        for v in 0..n {
90            if mask & (1 << v) == 0 {
91                continue; // v not in mask
92            }
93            if dp[mask][v] == f64::INFINITY {
94                continue;
95            }
96            // Try extending to each unvisited node u.
97            for u in 0..n {
98                if mask & (1 << u) != 0 {
99                    continue; // u already visited
100                }
101                let cost = dp[mask][v] + dist[v][u];
102                let next_mask = mask | (1 << u);
103                if cost < dp[next_mask][u] {
104                    dp[next_mask][u] = cost;
105                    parent[next_mask][u] = v;
106                }
107            }
108        }
109    }
110
111    // Find the minimum cost to complete the circuit: return to node 0.
112    let mut best_cost = f64::INFINITY;
113    let mut last_node = usize::MAX;
114
115    for v in 1..n {
116        if dp[full_mask][v] == f64::INFINITY {
117            continue;
118        }
119        let total = dp[full_mask][v] + dist[v][0];
120        if total < best_cost {
121            best_cost = total;
122            last_node = v;
123        }
124    }
125
126    if best_cost == f64::INFINITY {
127        return None; // No Hamiltonian circuit exists.
128    }
129
130    // Reconstruct path by following parent pointers.
131    let mut path = Vec::with_capacity(n + 1);
132    let mut mask = full_mask;
133    let mut curr = last_node;
134
135    while curr != usize::MAX && mask != 0 {
136        path.push(curr);
137        let prev = parent[mask][curr];
138        mask ^= 1 << curr;
139        curr = prev;
140    }
141
142    path.reverse();
143    path.push(0); // Return to start.
144
145    Some((best_cost, path))
146}
graph_advanced/tsp.rs

graph_advanced/
tsp.rs
