Christofides algorithm
You are encouraged to solve this task according to the task description, using any language you may know.
The Christofides algorithm is an approximation algorithm for solving the metric traveling salesman problem (TSP), a classic optimization problem in graph theory and computer science. Given a complete undirected graph with non-negative edge weights satisfying the triangle inequality, the algorithm produces a Hamiltonian circuit (a tour visiting each vertex exactly once) with a total weight at most 3/2 times the weight of the optimal tour.
- Problem Description
The metric TSP is defined as follows: given a set of n vertices and a distance function d(u, v) for each pair of vertices u and v, find a Hamiltonian circuit with the minimum total distance. The distance function must satisfy:
- d(u, v) ≥ 0 (non-negativity),
- d(u, v) = d(v, u) (symmetry),
- d(u, w) ≤ d(u, v) + d(v, w) (triangle inequality).
The TSP is NP-hard, meaning no known polynomial-time algorithm can solve it exactly for all cases unless P=NP. The Christofides algorithm provides a polynomial-time approximation with a guaranteed approximation ratio of 3/2 for metric TSP instances.
- Performance Guarantee
The algorithm guarantees that the weight of the resulting tour is at most 3/2 times the weight of the optimal TSP tour. This bound is derived as follows:
- The MST weight is at most the weight of the optimal tour (since removing one edge from the optimal tour yields a spanning tree).
- The minimum-weight perfect matching for the odd-degree vertices has a weight at most half the optimal tour weight (since the optimal tour induces a matching on these vertices).
- The Eulerian circuit’s weight is the sum of the MST and matching weights, and shortcutting does not increase the weight due to the triangle inequality.
Thus, the total weight is at most MST + Matching ≤ OPT + OPT/2 = 3/2 OPT, where OPT is the weight of the optimal tour.
- Complexity
The algorithm runs in polynomial time:
- MST computation: O(m log n) using Kruskal’s or Prim’s algorithm, where m is the number of edges and n is the number of vertices.
- Odd-degree vertex identification: O(n).
- Minimum-weight perfect matching: O(n^3) using algorithms like the blossom algorithm.
- Eulerian circuit construction: O(m).
- Shortcut to Hamiltonian circuit: O(n).
The overall time complexity is dominated by the matching step, resulting in O(n^3) for a complete graph where m = O(n^2).
- Limitations
While the Christofides algorithm provides a 3/2-approximation for metric TSP, it does not apply to non-metric instances (where the triangle inequality does not hold). Additionally, the 3/2 bound is not always tight; in some cases, the algorithm may produce a tour closer to the optimal. However, no polynomial-time algorithm is known to achieve a better approximation ratio for metric TSP unless P=NP.
- Applications
The Christofides algorithm is used in logistics, circuit design, and network optimization, where finding near-optimal tours is critical. It serves as a practical heuristic for problems like vehicle routing and scheduling, balancing computational efficiency with solution quality.
- Algorithm Overview
The Christofides algorithm, proposed by Nicos Christofides in 1976, proceeds as follows:
- Construct a minimum spanning tree (MST) of the graph using an algorithm like Kruskal's algorithm or Prim's algorithm.
- Identify the vertices with odd degree in the MST. Since a spanning tree has n-1 edges, the number of odd-degree vertices is even.
- Compute a minimum-weight perfect matching for these odd-degree vertices, forming a subgraph where each vertex has degree 2.
- Combine the MST and the matching to form a multigraph that is connected and has even degrees for all vertices.
- Find an Eulerian circuit in this multigraph (a circuit visiting each edge exactly once).
- Convert the Eulerian circuit into a Hamiltonian circuit by "shortcutting" repeated vertices, leveraging the triangle inequality to ensure the total weight does not increase.
C++
This is a translation of Python code
#include <iostream>
#include <vector>
#include <cmath>
#include <numeric>
#include <algorithm>
#include <map>
#include <list>
#include <stack>
#include <limits>
#include <random>
#include <utility> // For std::pair
#include <iomanip> // For std::fixed and std::setprecision
// --- Helper Structs ---
struct Point {
double x, y;
int id; // Original index
};
struct Edge {
int u, v;
double weight;
// For sorting edges by weight
bool operator<(const Edge& other) const {
return weight < other.weight;
}
};
// --- Helper Function to print vectors/lists ---
template <typename T>
void print_container(const T& container, const std::string& name) {
std::cout << name << ": [";
bool first = true;
for (const auto& item : container) {
if (!first) {
std::cout << ", ";
}
std::cout << item;
first = false;
}
std::cout << "]" << std::endl;
}
void print_edges(const std::vector<Edge>& edges, const std::string& name) {
std::cout << name << ": [";
bool first = true;
for (const auto& edge : edges) {
if (!first) {
std::cout << ", ";
}
std::cout << "(" << edge.u << ", " << edge.v << ", " << std::fixed << std::setprecision(2) << edge.weight << ")";
first = false;
}
std::cout << "]" << std::endl;
}
void print_graph(const std::vector<std::vector<double>>& graph, const std::string& name) {
std::cout << name << ": {" << std::endl;
int n = graph.size();
for (int i = 0; i < n; ++i) {
std::cout << " " << i << ": {";
bool first = true;
for (int j = 0; j < n; ++j) {
if (i != j) {
if (!first) {
std::cout << ", ";
}
std::cout << j << ": " << std::fixed << std::setprecision(2) << graph[i][j];
first = false;
}
}
std::cout << "}" << (i == n-1 ? "" : ",") << std::endl;
}
std::cout << "}" << std::endl;
}
// --- Euclidean Distance ---
double get_length(const Point& p1, const Point& p2) {
double dx = p1.x - p2.x;
double dy = p1.y - p2.y;
return std::sqrt(dx * dx + dy * dy);
}
// --- Build Complete Graph (Adjacency Matrix) ---
std::vector<std::vector<double>> build_graph(const std::vector<Point>& data) {
int n = data.size();
std::vector<std::vector<double>> graph(n, std::vector<double>(n, 0.0));
for (int i = 0; i < n; ++i) {
for (int j = i + 1; j < n; ++j) { // Only calculate upper triangle + diagonal is 0
double dist = get_length(data[i], data[j]);
graph[i][j] = dist;
graph[j][i] = dist; // Symmetric graph
}
}
return graph;
}
// --- Union-Find Data Structure ---
class UnionFind {
std::vector<int> parent;
std::vector<int> rank; // Or size
public:
UnionFind(int n) {
parent.resize(n);
std::iota(parent.begin(), parent.end(), 0); // Fill with 0, 1, 2, ...
rank.assign(n, 0);
}
int find(int i) {
if (parent[i] == i)
return i;
// Path compression
return parent[i] = find(parent[i]);
}
void unite(int i, int j) {
int root_i = find(i);
int root_j = find(j);
if (root_i != root_j) {
// Union by rank
if (rank[root_i] < rank[root_j]) {
parent[root_i] = root_j;
} else if (rank[root_i] > rank[root_j]) {
parent[root_j] = root_i;
} else {
parent[root_j] = root_i;
rank[root_i]++;
}
}
}
};
// --- Minimum Spanning Tree (Kruskal's Algorithm) ---
std::vector<Edge> minimum_spanning_tree(const std::vector<std::vector<double>>& graph) {
int n = graph.size();
if (n == 0) return {};
std::vector<Edge> edges;
for (int i = 0; i < n; ++i) {
for (int j = i + 1; j < n; ++j) { // Avoid duplicates and self-loops
edges.push_back({i, j, graph[i][j]});
}
}
// Sort edges by weight
std::sort(edges.begin(), edges.end());
std::vector<Edge> mst;
UnionFind uf(n);
int edges_count = 0;
for (const auto& edge : edges) {
if (uf.find(edge.u) != uf.find(edge.v)) {
mst.push_back(edge);
uf.unite(edge.u, edge.v);
edges_count++;
if (edges_count == n - 1) { // Optimization: MST has n-1 edges
break;
}
}
}
return mst;
}
// --- Find Vertices with Odd Degree in MST ---
std::vector<int> find_odd_vertexes(const std::vector<Edge>& mst, int n) {
std::vector<int> degree(n, 0);
for (const auto& edge : mst) {
degree[edge.u]++;
degree[edge.v]++;
}
std::vector<int> odd_vertices;
for (int i = 0; i < n; ++i) {
if (degree[i] % 2 != 0) {
odd_vertices.push_back(i);
}
}
return odd_vertices;
}
// --- Minimum Weight Matching (Greedy Heuristic) ---
// Note: This is a simplified greedy approach, not a true min-weight perfect matching.
// It modifies the mst vector by adding matching edges.
void minimum_weight_matching(std::vector<Edge>& mst, const std::vector<std::vector<double>>& graph, std::vector<int>& odd_vertices) {
// Use a copy to allow modification while iterating
std::vector<int> current_odd = odd_vertices;
// Shuffle for randomness (like the Python version)
// This makes the output non-deterministic unless seeded
std::random_device rd;
std::mt19937 g(rd()); // Mersenne Twister random number generator
std::shuffle(current_odd.begin(), current_odd.end(), g);
// Keep track of vertices already matched in this phase
std::vector<bool> matched(graph.size(), false);
for (int i = 0; i < current_odd.size(); ++i) {
int v = current_odd[i];
if (matched[v]) continue; // Skip if already matched
double min_length = std::numeric_limits<double>::infinity();
int closest_u = -1;
// Find the closest unmatched odd vertex
for (int j = i + 1; j < current_odd.size(); ++j) {
int u = current_odd[j];
if (!matched[u]) { // Check if 'u' is available
if (graph[v][u] < min_length) {
min_length = graph[v][u];
closest_u = u;
}
}
}
if (closest_u != -1) {
// Add the matching edge to the MST list (now a multigraph)
mst.push_back({v, closest_u, min_length});
matched[v] = true;
matched[closest_u] = true; // Mark both as matched
}
// Note: In a perfect matching on an even number of vertices,
// every vertex should find a match. If closest_u remains -1,
// something is wrong (e.g., odd number of odd_vertices input?)
// Christofides guarantees an even number of odd-degree vertices.
}
}
// --- Find Eulerian Tour (Hierholzer's Algorithm) ---
std::vector<int> find_eulerian_tour(const std::vector<Edge>& matched_mst, int n) {
if (matched_mst.empty()) return {};
// Build adjacency list representation of the multigraph (MST + matching)
// Use std::list for efficient removal of edges during tour construction
std::vector<std::list<std::pair<int, const Edge*>>> adj(n); // Store neighbor and pointer to original edge
std::map<const Edge*, bool> edge_used; // Track used edges (since multigraph)
for (const auto& edge : matched_mst) {
adj[edge.u].push_back({edge.v, &edge});
adj[edge.v].push_back({edge.u, &edge});
edge_used[&edge] = false;
}
std::vector<int> tour;
std::stack<int> current_path;
// Start at any vertex with edges (e.g., the first vertex of the first edge)
int start_node = matched_mst[0].u;
current_path.push(start_node);
int current_node = start_node;
while (!current_path.empty()) {
current_node = current_path.top();
bool found_edge = false;
// Find an unused edge from the current node
for (auto it = adj[current_node].begin(); it != adj[current_node].end(); ) {
int neighbor = it->first;
const Edge* edge_ptr = it->second;
if (!edge_used[edge_ptr]) {
edge_used[edge_ptr] = true; // Mark edge as used
// Push neighbor onto stack and move to it
current_path.push(neighbor);
// Remove the edge from the neighbor's list as well
// This part is tricky with pointers; let's simplify
// by just relying on the edge_used map. A more robust
// Hierholzer implementation might store iterators.
// We don't strictly need to remove from adj list if using edge_used map
found_edge = true;
break; // Move to the neighbor
} else {
++it; // Check next edge
}
}
// If no unused edge was found from current_node, backtrack
if (!found_edge) {
tour.push_back(current_path.top());
current_path.pop();
}
}
// The tour is constructed in reverse order
std::reverse(tour.begin(), tour.end());
return tour;
}
// --- Main TSP Function (Christofides Approximation) ---
std::pair<double, std::vector<int>> tsp(const std::vector<Point>& data) {
int n = data.size();
if (n == 0) {
return {0.0, {}};
}
if (n == 1) {
return {0.0, {data[0].id}}; // Or just {0} if using 0-based indexing internally
}
// Build a graph
std::vector<std::vector<double>> G = build_graph(data);
print_graph(G, "Graph");
// Build a minimum spanning tree
std::vector<Edge> MSTree = minimum_spanning_tree(G);
print_edges(MSTree, "MSTree");
// Find odd degree vertices
std::vector<int> odd_vertexes = find_odd_vertexes(MSTree, n);
print_container(odd_vertexes, "Odd vertexes in MSTree");
// Add minimum weight matching edges (using greedy heuristic)
// Note: This modifies MSTree, turning it into a multigraph representation
minimum_weight_matching(MSTree, G, odd_vertexes);
print_edges(MSTree, "Minimum weight matching (MST + Matching Edges)"); // Now contains MST + Matching
// Find an Eulerian tour in the combined graph
std::vector<int> eulerian_tour = find_eulerian_tour(MSTree, n);
print_container(eulerian_tour, "Eulerian tour");
// --- Create Hamiltonian Circuit by Skipping Visited Nodes ---
if (eulerian_tour.empty()) {
std::cerr << "Error: Eulerian tour could not be found." << std::endl;
return {-1.0, {}}; // Indicate error
}
std::vector<int> path;
double length = 0.0;
std::vector<bool> visited(n, false);
int current = eulerian_tour[0];
path.push_back(current);
visited[current] = true;
for (size_t i = 1; i < eulerian_tour.size(); ++i) {
int v = eulerian_tour[i];
if (!visited[v]) {
path.push_back(v);
visited[v] = true;
length += G[current][v]; // Add distance from previous node in path
current = v; // Update current node in path
}
}
// Add the edge back to the start
length += G[current][path[0]];
path.push_back(path[0]); // Complete the cycle
// Convert path indices back to original IDs if needed (assuming Point.id holds them)
// If the input 'data' vector's order corresponds to IDs 0..N-1, no conversion is needed.
// If Point.id was used differently, you'd map back here.
print_container(path, "Result path");
std::cout << "Result length of the path: " << std::fixed << std::setprecision(2) << length << std::endl;
return {length, path};
}
int main() {
// Input data matching the Python example
// Store as Point structs, retaining original order implicitly as ID 0..N-1
std::vector<std::pair<double, double>> raw_data = {
{1380, 939}, {2848, 96}, {3510, 1671}, {457, 334}, {3888, 666}, {984, 965}, {2721, 1482}, {1286, 525},
{2716, 1432}, {738, 1325}, {1251, 1832}, {2728, 1698}, {3815, 169}, {3683, 1533}, {1247, 1945}, {123, 862},
{1234, 1946}, {252, 1240}, {611, 673}, {2576, 1676}, {928, 1700}, {53, 857}, {1807, 1711}, {274, 1420},
{2574, 946}, {178, 24}, {2678, 1825}, {1795, 962}, {3384, 1498}, {3520, 1079}, {1256, 61}, {1424, 1728},
{3913, 192}, {3085, 1528}, {2573, 1969}, {463, 1670}, {3875, 598}, {298, 1513}, {3479, 821}, {2542, 236},
{3955, 1743}, {1323, 280}, {3447, 1830}, {2936, 337}, {1621, 1830}, {3373, 1646}, {1393, 1368},
{3874, 1318}, {938, 955}, {3022, 474}, {2482, 1183}, {3854, 923}, {376, 825}, {2519, 135}, {2945, 1622},
{953, 268}, {2628, 1479}, {2097, 981}, {890, 1846}, {2139, 1806}, {2421, 1007}, {2290, 1810}, {1115, 1052},
{2588, 302}, {327, 265}, {241, 341}, {1917, 687}, {2991, 792}, {2573, 599}, {19, 674}, {3911, 1673},
{872, 1559}, {2863, 558}, {929, 1766}, {839, 620}, {3893, 102}, {2178, 1619}, {3822, 899}, {378, 1048},
{1178, 100}, {2599, 901}, {3416, 143}, {2961, 1605}, {611, 1384}, {3113, 885}, {2597, 1830}, {2586, 1286},
{161, 906}, {1429, 134}, {742, 1025}, {1625, 1651}, {1187, 706}, {1787, 1009}, {22, 987}, {3640, 43},
{3756, 882}, {776, 392}, {1724, 1642}, {198, 1810}, {3950, 1558}
};
std::vector<Point> data_points;
data_points.reserve(raw_data.size());
for(int i = 0; i < raw_data.size(); ++i) {
data_points.push_back({raw_data[i].first, raw_data[i].second, i});
}
tsp(data_points);
return 0;
}
- Output:
Graph: ...... Eulerian tour: [14, 16, 98, 35, 37, 23, 17, 78, 52, 87, 15, 21, 69, 21, 93, 52, 18, 74, 96, 3, 64, 65, 64, 25, 96, 55, 79, 30, 88, 41, 7, 91, 0, 62, 5, 48, 89, 9, 83, 9, 71, 20, 73, 58, 73, 10, 31, 90, 97, 22, 59, 76, 59, 61, 85, 26, 11, 54, 82, 33, 28, 45, 2, 13, 70, 40, 70, 99, 47, 2, 42, 34, 85, 11, 19, 56, 6, 8, 86, 50, 60, 57, 27, 66, 27, 92, 60, 24, 80, 68, 72, 49, 43, 1, 53, 39, 63, 68, 72, 67, 84, 38, 29, 38, 95, 77, 4, 36, 32, 75, 81, 94, 12, 32, 77, 51, 44, 90, 46, 31, 10, 14] Result path: [14, 16, 98, 35, 37, 23, 17, 78, 52, 87, 15, 21, 69, 93, 18, 74, 96, 3, 64, 65, 25, 55, 79, 30, 88, 41, 7, 91, 0, 62, 5, 48, 89, 9, 83, 71, 20, 73, 58, 10, 31, 90, 97, 22, 59, 76, 61, 85, 26, 11, 54, 82, 33, 28, 45, 2, 13, 70, 40, 99, 47, 42, 34, 19, 56, 6, 8, 86, 50, 60, 57, 27, 66, 92, 24, 80, 68, 72, 49, 43, 1, 53, 39, 63, 67, 84, 38, 29, 95, 77, 4, 36, 32, 75, 81, 94, 12, 51, 44, 46, 14] Result length of the path: 27920.47
Go
package main
import (
"fmt"
"math"
"math/rand"
"sort"
"time"
)
// --- Helper Structs ---
type Point struct {
x, y float64
id int // Original index
}
type Edge struct {
u, v int
weight float64
}
// --- Helper Function to print vectors/lists ---
func printContainer[T any](container []T, name string) {
fmt.Printf("%s: [", name)
for i, item := range container {
if i > 0 {
fmt.Print(", ")
}
fmt.Print(item)
}
fmt.Println("]")
}
func printEdges(edges []Edge, name string) {
fmt.Printf("%s: [", name)
for i, edge := range edges {
if i > 0 {
fmt.Print(", ")
}
fmt.Printf("(%d, %d, %.2f)", edge.u, edge.v, edge.weight)
}
fmt.Println("]")
}
func printGraph(graph [][]float64, name string) {
fmt.Printf("%s: {\n", name)
n := len(graph)
for i := 0; i < n; i++ {
fmt.Printf(" %d: {", i)
first := true
for j := 0; j < n; j++ {
if i != j {
if !first {
fmt.Print(", ")
}
fmt.Printf("%d: %.2f", j, graph[i][j])
first = false
}
}
comma := ""
if i != n-1 {
comma = ","
}
fmt.Printf("}%s\n", comma)
}
fmt.Println("}")
}
// --- Euclidean Distance ---
func getLength(p1, p2 Point) float64 {
dx := p1.x - p2.x
dy := p1.y - p2.y
return math.Sqrt(dx*dx + dy*dy)
}
// --- Build Complete Graph (Adjacency Matrix) ---
func buildGraph(data []Point) [][]float64 {
n := len(data)
graph := make([][]float64, n)
for i := range graph {
graph[i] = make([]float64, n)
}
for i := 0; i < n; i++ {
for j := i + 1; j < n; j++ { // Only calculate upper triangle + diagonal is 0
dist := getLength(data[i], data[j])
graph[i][j] = dist
graph[j][i] = dist // Symmetric graph
}
}
return graph
}
// --- Union-Find Data Structure ---
type UnionFind struct {
parent []int
rank []int
}
func NewUnionFind(n int) *UnionFind {
parent := make([]int, n)
rank := make([]int, n)
for i := range parent {
parent[i] = i
}
return &UnionFind{parent, rank}
}
func (uf *UnionFind) Find(i int) int {
if uf.parent[i] == i {
return i
}
// Path compression
uf.parent[i] = uf.Find(uf.parent[i])
return uf.parent[i]
}
func (uf *UnionFind) Unite(i, j int) {
rootI := uf.Find(i)
rootJ := uf.Find(j)
if rootI != rootJ {
// Union by rank
if uf.rank[rootI] < uf.rank[rootJ] {
uf.parent[rootI] = rootJ
} else if uf.rank[rootI] > uf.rank[rootJ] {
uf.parent[rootJ] = rootI
} else {
uf.parent[rootJ] = rootI
uf.rank[rootI]++
}
}
}
// --- Minimum Spanning Tree (Kruskal's Algorithm) ---
func minimumSpanningTree(graph [][]float64) []Edge {
n := len(graph)
if n == 0 {
return []Edge{}
}
edges := []Edge{}
for i := 0; i < n; i++ {
for j := i + 1; j < n; j++ { // Avoid duplicates and self-loops
edges = append(edges, Edge{i, j, graph[i][j]})
}
}
// Sort edges by weight
sort.Slice(edges, func(i, j int) bool {
return edges[i].weight < edges[j].weight
})
mst := []Edge{}
uf := NewUnionFind(n)
edgesCount := 0
for _, edge := range edges {
if uf.Find(edge.u) != uf.Find(edge.v) {
mst = append(mst, edge)
uf.Unite(edge.u, edge.v)
edgesCount++
if edgesCount == n-1 { // Optimization: MST has n-1 edges
break
}
}
}
return mst
}
// --- Find Vertices with Odd Degree in MST ---
func findOddVertexes(mst []Edge, n int) []int {
degree := make([]int, n)
for _, edge := range mst {
degree[edge.u]++
degree[edge.v]++
}
oddVertices := []int{}
for i := 0; i < n; i++ {
if degree[i]%2 != 0 {
oddVertices = append(oddVertices, i)
}
}
return oddVertices
}
// --- Minimum Weight Matching (Greedy Heuristic) ---
func minimumWeightMatching(mst []Edge, graph [][]float64, oddVertices []int) []Edge {
// Create a copy to allow modification while iterating
currentOdd := make([]int, len(oddVertices))
copy(currentOdd, oddVertices)
// Shuffle for randomness
rand.Seed(time.Now().UnixNano())
rand.Shuffle(len(currentOdd), func(i, j int) {
currentOdd[i], currentOdd[j] = currentOdd[j], currentOdd[i]
})
// Keep track of vertices already matched in this phase
matched := make([]bool, len(graph))
result := make([]Edge, len(mst))
copy(result, mst)
for i := 0; i < len(currentOdd); i++ {
v := currentOdd[i]
if matched[v] {
continue // Skip if already matched
}
minLength := math.Inf(1)
closestU := -1
// Find the closest unmatched odd vertex
for j := i + 1; j < len(currentOdd); j++ {
u := currentOdd[j]
if !matched[u] { // Check if 'u' is available
if graph[v][u] < minLength {
minLength = graph[v][u]
closestU = u
}
}
}
if closestU != -1 {
// Add the matching edge to the MST list (now a multigraph)
result = append(result, Edge{v, closestU, minLength})
matched[v] = true
matched[closestU] = true // Mark both as matched
}
}
return result
}
// --- Find Eulerian Tour (Hierholzer's Algorithm) ---
func findEulerianTour(matchedMst []Edge, n int) []int {
if len(matchedMst) == 0 {
return []int{}
}
// Build adjacency list representation of the multigraph (MST + matching)
type EdgeInfo struct {
neighbor int
edgePtr *Edge
}
adj := make([][]EdgeInfo, n)
edgeUsed := make(map[*Edge]bool)
for i := range matchedMst {
edge := &matchedMst[i]
adj[edge.u] = append(adj[edge.u], EdgeInfo{edge.v, edge})
adj[edge.v] = append(adj[edge.v], EdgeInfo{edge.u, edge})
edgeUsed[edge] = false
}
tour := []int{}
currentPath := []int{}
// Start at any vertex with edges (e.g., the first vertex of the first edge)
startNode := matchedMst[0].u
currentPath = append(currentPath, startNode)
for len(currentPath) > 0 {
currentNode := currentPath[len(currentPath)-1]
foundEdge := false
// Find an unused edge from the current node
for i := 0; i < len(adj[currentNode]); i++ {
neighbor := adj[currentNode][i].neighbor
edgePtr := adj[currentNode][i].edgePtr
if !edgeUsed[edgePtr] {
edgeUsed[edgePtr] = true // Mark edge as used
// Push neighbor onto stack and move to it
currentPath = append(currentPath, neighbor)
foundEdge = true
break // Move to the neighbor
}
}
// If no unused edge was found from currentNode, backtrack
if !foundEdge {
tour = append(tour, currentPath[len(currentPath)-1])
currentPath = currentPath[:len(currentPath)-1]
}
}
// Reverse the tour
for i, j := 0, len(tour)-1; i < j; i, j = i+1, j-1 {
tour[i], tour[j] = tour[j], tour[i]
}
return tour
}
// --- Main TSP Function (Christofides Approximation) ---
func tsp(data []Point) (float64, []int) {
//fmt.Printf("%s\n", "step in `tsp` function")
n := len(data)
if n == 0 {
return 0.0, []int{}
}
if n == 1 {
return 0.0, []int{data[0].id}
}
// Build a graph
G := buildGraph(data)
//printGraph(G, "Graph")
// Build a minimum spanning tree
MSTree := minimumSpanningTree(G)
printEdges(MSTree, "MSTree")
// Find odd degree vertices
oddVertexes := findOddVertexes(MSTree, n)
printContainer(oddVertexes, "Odd vertexes in MSTree")
// Add minimum weight matching edges (using greedy heuristic)
// Note: This returns a new slice containing MST + matching edges
MSTreeWithMatching := minimumWeightMatching(MSTree, G, oddVertexes)
printEdges(MSTreeWithMatching, "Minimum weight matching (MST + Matching Edges)")
// Find an Eulerian tour in the combined graph
eulerianTour := findEulerianTour(MSTreeWithMatching, n)
printContainer(eulerianTour, "Eulerian tour")
// --- Create Hamiltonian Circuit by Skipping Visited Nodes ---
if len(eulerianTour) == 0 {
fmt.Println("Error: Eulerian tour could not be found.")
return -1.0, []int{}
}
path := []int{}
length := 0.0
visited := make([]bool, n)
current := eulerianTour[0]
path = append(path, current)
visited[current] = true
for i := 1; i < len(eulerianTour); i++ {
v := eulerianTour[i]
if !visited[v] {
path = append(path, v)
visited[v] = true
length += G[current][v] // Add distance from previous node in path
current = v // Update current node in path
}
}
// Add the edge back to the start
length += G[current][path[0]]
path = append(path, path[0]) // Complete the cycle
printContainer(path, "Result path")
fmt.Printf("Result length of the path: %.2f\n", length)
return length, path
}
func main() {
// Input data matching the C++ example
rawData := [][]float64{
{1380, 939}, {2848, 96}, {3510, 1671}, {457, 334}, {3888, 666}, {984, 965}, {2721, 1482}, {1286, 525},
{2716, 1432}, {738, 1325}, {1251, 1832}, {2728, 1698}, {3815, 169}, {3683, 1533}, {1247, 1945}, {123, 862},
{1234, 1946}, {252, 1240}, {611, 673}, {2576, 1676}, {928, 1700}, {53, 857}, {1807, 1711}, {274, 1420},
{2574, 946}, {178, 24}, {2678, 1825}, {1795, 962}, {3384, 1498}, {3520, 1079}, {1256, 61}, {1424, 1728},
{3913, 192}, {3085, 1528}, {2573, 1969}, {463, 1670}, {3875, 598}, {298, 1513}, {3479, 821}, {2542, 236},
{3955, 1743}, {1323, 280}, {3447, 1830}, {2936, 337}, {1621, 1830}, {3373, 1646}, {1393, 1368},
{3874, 1318}, {938, 955}, {3022, 474}, {2482, 1183}, {3854, 923}, {376, 825}, {2519, 135}, {2945, 1622},
{953, 268}, {2628, 1479}, {2097, 981}, {890, 1846}, {2139, 1806}, {2421, 1007}, {2290, 1810}, {1115, 1052},
{2588, 302}, {327, 265}, {241, 341}, {1917, 687}, {2991, 792}, {2573, 599}, {19, 674}, {3911, 1673},
{872, 1559}, {2863, 558}, {929, 1766}, {839, 620}, {3893, 102}, {2178, 1619}, {3822, 899}, {378, 1048},
{1178, 100}, {2599, 901}, {3416, 143}, {2961, 1605}, {611, 1384}, {3113, 885}, {2597, 1830}, {2586, 1286},
{161, 906}, {1429, 134}, {742, 1025}, {1625, 1651}, {1187, 706}, {1787, 1009}, {22, 987}, {3640, 43},
{3756, 882}, {776, 392}, {1724, 1642}, {198, 1810}, {3950, 1558},
}
dataPoints := make([]Point, len(rawData))
for i, point := range rawData {
dataPoints[i] = Point{point[0], point[1], i}
}
tsp(dataPoints)
}
- Output:
MSTree: [(14, 16, 13.04), (54, 82, 23.35), (51, 77, 40.00), (5, 48, 47.07), (27, 92, 47.68), (6, 8, 50.25), (24, 80, 51.48), (15, 87, 58.14), (20, 73, 66.01), (77, 95, 68.15), (4, 36, 69.23), (15, 21, 70.18), (39, 63, 80.45), (26, 85, 81.15), (40, 70, 82.68), (30, 79, 87.21), (58, 73, 89.00), (32, 75, 92.20), (6, 56, 93.05), (23, 37, 96.05), (90, 97, 99.41), (12, 32, 100.66), (39, 53, 103.59), (22, 97, 107.94), (10, 14, 113.07), (64, 65, 114.77), (70, 99, 121.43), (21, 93, 133.65), (11, 26, 136.49), (2, 45, 139.26), (9, 83, 140.04), (34, 85, 141.06), (33, 82, 145.96), (50, 86, 146.37), (3, 64, 147.18), (28, 45, 148.41), (59, 61, 151.05), (20, 71, 151.71), (67, 84, 153.40), (11, 19, 153.58), (5, 62, 157.26), (43, 49, 161.76), (24, 60, 164.71), (2, 42, 171.03), (44, 90, 179.04), (49, 72, 179.82), (41, 88, 180.42), (17, 23, 181.34), (21, 69, 186.13), (50, 60, 186.27), (30, 88, 187.77), (59, 76, 191.02), (8, 86, 195.49), (10, 31, 201.85), (19, 56, 203.75), (7, 91, 206.31), (48, 89, 208.12), (31, 90, 215.24), (12, 94, 215.64), (55, 96, 216.11), (2, 13, 221.30), (52, 78, 223.01), (35, 37, 227.76), (17, 78, 229.65), (52, 87, 229.75), (11, 54, 229.92), (18, 74, 234.08), (74, 96, 236.54), (4, 77, 242.17), (81, 94, 245.31), (7, 41, 247.78), (47, 99, 251.75), (1, 43, 256.56), (29, 38, 261.24), (67, 72, 266.72), (13, 70, 267.55), (9, 71, 269.65), (18, 52, 279.87), (55, 79, 280.80), (25, 64, 283.34), (38, 95, 283.64), (0, 62, 288.09), (68, 72, 292.88), (63, 68, 297.38), (35, 98, 299.71), (9, 89, 300.03), (28, 33, 300.50), (27, 66, 300.85), (0, 91, 302.55), (27, 57, 302.60), (68, 80, 303.12), (61, 85, 307.65), (3, 96, 324.23), (57, 60, 325.04), (10, 73, 328.69), (22, 59, 345.32), (31, 46, 361.33), (38, 84, 371.55), (32, 36, 407.77)] Odd vertexes in MSTree: [1, 2, 9, 10, 11, 16, 21, 25, 27, 29, 31, 32, 34, 38, 40, 42, 44, 46, 47, 51, 52, 53, 58, 59, 60, 64, 65, 66, 68, 69, 70, 72, 73, 75, 76, 77, 81, 83, 85, 90, 92, 93, 96, 98] Minimum weight matching (MST + Matching Edges): [(14, 16, 13.04), (54, 82, 23.35), (51, 77, 40.00), (5, 48, 47.07), (27, 92, 47.68), (6, 8, 50.25), (24, 80, 51.48), (15, 87, 58.14), (20, 73, 66.01), (77, 95, 68.15), (4, 36, 69.23), (15, 21, 70.18), (39, 63, 80.45), (26, 85, 81.15), (40, 70, 82.68), (30, 79, 87.21), (58, 73, 89.00), (32, 75, 92.20), (6, 56, 93.05), (23, 37, 96.05), (90, 97, 99.41), (12, 32, 100.66), (39, 53, 103.59), (22, 97, 107.94), (10, 14, 113.07), (64, 65, 114.77), (70, 99, 121.43), (21, 93, 133.65), (11, 26, 136.49), (2, 45, 139.26), (9, 83, 140.04), (34, 85, 141.06), (33, 82, 145.96), (50, 86, 146.37), (3, 64, 147.18), (28, 45, 148.41), (59, 61, 151.05), (20, 71, 151.71), (67, 84, 153.40), (11, 19, 153.58), (5, 62, 157.26), (43, 49, 161.76), (24, 60, 164.71), (2, 42, 171.03), (44, 90, 179.04), (49, 72, 179.82), (41, 88, 180.42), (17, 23, 181.34), (21, 69, 186.13), (50, 60, 186.27), (30, 88, 187.77), (59, 76, 191.02), (8, 86, 195.49), (10, 31, 201.85), (19, 56, 203.75), (7, 91, 206.31), (48, 89, 208.12), (31, 90, 215.24), (12, 94, 215.64), (55, 96, 216.11), (2, 13, 221.30), (52, 78, 223.01), (35, 37, 227.76), (17, 78, 229.65), (52, 87, 229.75), (11, 54, 229.92), (18, 74, 234.08), (74, 96, 236.54), (4, 77, 242.17), (81, 94, 245.31), (7, 41, 247.78), (47, 99, 251.75), (1, 43, 256.56), (29, 38, 261.24), (67, 72, 266.72), (13, 70, 267.55), (9, 71, 269.65), (18, 52, 279.87), (55, 79, 280.80), (25, 64, 283.34), (38, 95, 283.64), (0, 62, 288.09), (68, 72, 292.88), (63, 68, 297.38), (35, 98, 299.71), (9, 89, 300.03), (28, 33, 300.50), (27, 66, 300.85), (0, 91, 302.55), (27, 57, 302.60), (68, 80, 303.12), (61, 85, 307.65), (3, 96, 324.23), (57, 60, 325.04), (10, 73, 328.69), (22, 59, 345.32), (31, 46, 361.33), (38, 84, 371.55), (32, 36, 407.77), (76, 59, 191.02), (73, 58, 89.00), (44, 90, 179.04), (85, 34, 141.06), (51, 77, 40.00), (1, 53, 331.30), (27, 92, 47.68), (31, 10, 201.85), (2, 42, 171.03), (46, 16, 599.47), (69, 21, 186.13), (98, 83, 593.33), (60, 68, 435.39), (93, 52, 389.31), (65, 64, 114.77), (9, 96, 933.77), (40, 70, 82.68), (81, 75, 478.76), (29, 38, 261.24), (11, 72, 1147.97), (47, 32, 1126.68), (66, 25, 1861.10)] Eulerian tour: [14, 16, 46, 31, 10, 73, 58, 73, 20, 71, 9, 83, 98, 35, 37, 23, 17, 78, 52, 87, 15, 21, 69, 21, 93, 52, 18, 74, 96, 55, 79, 30, 88, 41, 7, 91, 0, 62, 5, 48, 89, 9, 96, 3, 64, 65, 64, 25, 66, 27, 92, 27, 57, 60, 24, 80, 68, 72, 49, 43, 1, 53, 39, 63, 68, 60, 50, 86, 8, 6, 56, 19, 11, 54, 82, 33, 28, 45, 2, 42, 2, 13, 70, 40, 70, 99, 47, 32, 75, 81, 94, 12, 32, 36, 4, 77, 51, 77, 95, 38, 29, 38, 84, 67, 72, 11, 26, 85, 34, 85, 61, 59, 76, 59, 22, 97, 90, 44, 90, 31, 10, 14] Result path: [14, 16, 46, 31, 10, 73, 58, 20, 71, 9, 83, 98, 35, 37, 23, 17, 78, 52, 87, 15, 21, 69, 93, 18, 74, 96, 55, 79, 30, 88, 41, 7, 91, 0, 62, 5, 48, 89, 3, 64, 65, 25, 66, 27, 92, 57, 60, 24, 80, 68, 72, 49, 43, 1, 53, 39, 63, 50, 86, 8, 6, 56, 19, 11, 54, 82, 33, 28, 45, 2, 42, 13, 70, 40, 99, 47, 32, 75, 81, 94, 12, 36, 4, 77, 51, 95, 38, 29, 84, 67, 26, 85, 34, 61, 59, 76, 22, 97, 90, 44, 14] Result length of the path: 26046.90
Java
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashSet;
import java.util.List;
import java.util.Optional;
import java.util.Set;
import java.util.Stack;
import java.util.function.BiFunction;
import java.util.stream.Collectors;
import java.util.stream.IntStream;
public class ChristofidesAlgorithm {
public static void main(String[] args) {
List<Pair> data = List.of(
new Pair(1380, 939), new Pair(2848, 96), new Pair(3510, 1671), new Pair(457, 334), new Pair(3888, 666),
new Pair(984, 965), new Pair(2721, 1482), new Pair(1286, 525), new Pair(2716, 1432), new Pair(738, 1325),
new Pair(1251, 1832), new Pair(2728, 1698), new Pair(3815, 169), new Pair(3683, 1533), new Pair(1247,1945),
new Pair(123, 862), new Pair(1234, 1946), new Pair(252, 1240), new Pair(611, 673), new Pair(2576, 1676),
new Pair(928, 1700), new Pair(53, 857), new Pair(1807, 1711), new Pair(274, 1420), new Pair(2574, 946),
new Pair(178, 24), new Pair(2678, 1825), new Pair(1795, 962), new Pair(3384, 1498), new Pair(3520, 1079),
new Pair(1256, 61), new Pair(1424, 1728), new Pair(3913, 192), new Pair(3085, 1528), new Pair(2573, 1969),
new Pair(463, 1670), new Pair(3875, 598), new Pair(298, 1513), new Pair(3479, 821), new Pair(2542, 236),
new Pair(3955, 1743), new Pair(1323, 280), new Pair(3447, 1830), new Pair(2936, 337), new Pair(1621, 1830),
new Pair(3373, 1646), new Pair(1393, 1368), new Pair(3874, 1318), new Pair(938, 955), new Pair(3022, 474),
new Pair(2482, 1183), new Pair(3854, 923), new Pair(376, 825), new Pair(2519, 135), new Pair(2945, 1622),
new Pair(953, 268), new Pair(2628, 1479), new Pair(2097, 981), new Pair(890, 1846), new Pair(2139, 1806),
new Pair(2421, 1007), new Pair(2290, 1810), new Pair(1115, 1052), new Pair(2588, 302), new Pair(327, 265),
new Pair(241, 341), new Pair(1917, 687), new Pair(2991, 792), new Pair(2573, 599), new Pair(19, 674),
new Pair(3911, 1673), new Pair(872, 1559), new Pair(2863, 558), new Pair(929, 1766), new Pair(839, 620),
new Pair(3893, 102), new Pair(2178, 1619), new Pair(3822, 899), new Pair(378, 1048), new Pair(1178, 100),
new Pair(2599, 901), new Pair(3416, 143), new Pair(2961, 1605), new Pair(611, 1384), new Pair(3113, 885),
new Pair(2597, 1830), new Pair(2586, 1286), new Pair(161, 906), new Pair(1429, 134), new Pair(742, 1025),
new Pair(1625, 1651), new Pair(1187, 706), new Pair(1787, 1009), new Pair(22, 987), new Pair(3640, 43),
new Pair(3756, 882), new Pair(776, 392), new Pair(1724, 1642), new Pair(198, 1810), new Pair(3950, 1558)
);
List<Point> points = IntStream.range(0, data.size()).mapToObj( i -> new Point(data.get(i), i) ).toList();
christofidesPath(points);
}
// Display and return an approximation to the optimum travelling salesman path using the Christofides algorithm
private static Result christofidesPath(List<Point> points) {
if ( points.isEmpty() ) {
return new Result(Collections.emptyList(), 0.0);
}
if ( points.size() == 1 ) {
return new Result(List.of(points.getFirst().id), 0.0);
}
Graph graph = new Graph(points);
graph.display();
List<Edge> minimumSpanningTree = graph.minimumSpanningTree();
System.out.println("Minimum spanning tree: " + minimumSpanningTree + System.lineSeparator());
List<Integer> oddVertices = graph.oddVertices(minimumSpanningTree);
System.out.println("Odd vertices in minimum spanning tree: " + oddVertices + System.lineSeparator());
List<Edge> minimumWeightMatching = graph.minimumWeightMatching(minimumSpanningTree, oddVertices);
System.out.println("Minimum weight matching: " + minimumWeightMatching + System.lineSeparator());
List<Integer> eulerianCircuit = graph.eulerianCircuit(minimumWeightMatching);
System.out.println("Eulerian circuit: " + eulerianCircuit + System.lineSeparator());
if ( eulerianCircuit.isEmpty() ) {
System.err.println("Error: Christofides path could not be found.");
return new Result(Collections.emptyList(), -1.0);
}
Result result = graph.hamiltonianCircuit(eulerianCircuit);
System.out.println("Result path: " + result.path + System.lineSeparator());
System.out.println(String.format("%s%.2f", "Length of the result path: ", result.length));
return result;
}
private record Pair(double x, double y) {}
private record Point(Pair pair, int id) {}
private record Edge(int u, int v, double weight) {
public String toString() {
return String.format("%s%d%s%d%s%.2f%s", "(", u, ", ", v , ", ", weight, ")");
}
}
private static final class Graph {
public Graph(List<Point> points) {
pointCount = points.size();
distanceLists = IntStream.range(0, pointCount).boxed()
.map( i -> new ArrayList<Double>(Collections.nCopies(pointCount, 0.0)) )
.collect(Collectors.toList());
BiFunction<Point, Point, Double> distance = (one, two) ->
Math.hypot(one.pair.x - two.pair.x, one.pair.y - two.pair.y);
IntStream.range(0, pointCount).forEach( i -> {
IntStream.range(i + 1, pointCount).forEach( j -> {
final double dist = distance.apply(points.get(i), points.get(j));
distanceLists.get(i).set(j, dist);
distanceLists.get(j).set(i, dist);
} );
} );
}
// Return the minimum spanning tree using Kruskal's algorithm
public List<Edge> minimumSpanningTree() {
List<Edge> edges = new ArrayList<Edge>();
if ( pointCount == 0 ) {
return edges;
}
IntStream.range(0, pointCount).forEach( i -> {
IntStream.range(i + 1, pointCount).forEach( j -> { // Avoids duplicates and self-loops
edges.addLast( new Edge(i, j, distanceLists.get(i).get(j)) );
} );
} );
// Sort edges by weight
Collections.sort(edges, (e1, e2) -> Double.compare(e1.weight, e2.weight));
List<Edge> minimumSpanningTree = new ArrayList<Edge>();
UnionFind unionFind = new UnionFind(pointCount);
int edgeCount = 0;
for ( Edge edge : edges ) {
if ( unionFind.unite(edge.u, edge.v) ) {
minimumSpanningTree.addLast(edge);
edgeCount += 1;
if ( edgeCount == pointCount - 1 ) {
break; // An optimisation, since the minimum spanning tree has n - 1 edges
}
}
}
return minimumSpanningTree;
}
// Return a list of vertices with odd degree in the minimum spanning tree
public List<Integer> oddVertices(List<Edge> minimumSpanningTree) {
List<Integer> degrees = new ArrayList<Integer>(Collections.nCopies(pointCount, 0));
minimumSpanningTree.forEach( edge -> {
degrees.set(edge.u, degrees.get(edge.u) + 1);
degrees.set(edge.v, degrees.get(edge.v) + 1);
} );
return IntStream.range(0, pointCount).filter( i -> degrees.get(i) % 2 == 1 ).boxed().toList();
}
// Return a minimum weight matching using a greedy heuristic
public List<Edge> minimumWeightMatching(List<Edge> minimumSpanningTree, List<Integer> oddVertices) {
List<Edge> minimumWeightMatching = new ArrayList<Edge>();
if ( oddVertices.isEmpty() ) {
return minimumWeightMatching;
}
// All elements of 'minimumSpanningTree' are included
minimumWeightMatching.addAll(minimumSpanningTree);
// Create a copy to prevent mutation of the argument 'oddVertices'
List<Integer> currentOdd = new ArrayList<Integer>(oddVertices);
Collections.shuffle(currentOdd); // Shuffle for randomness
// Maintain a record of the visited indices in the shuffled 'currentOdd' list
Set<Integer> visited = new HashSet<Integer>();
for ( int i = 0; i < currentOdd.size(); i++ ) {
if ( visited.contains(i) ) { // Omit a vertex which has already been processed
continue;
}
final int v = currentOdd.get(i);
double minimumDistance = Integer.MAX_VALUE;
Optional<Integer> closestUIndex = Optional.empty();
// Find the closest unmatched odd vertex occurring after 'v' in the shuffled 'currentOdd' list
for ( int j = i + 1; j < currentOdd.size(); j++ ) {
if ( ! visited.contains(j) ) { // Check whether a vertex is available
final int u = currentOdd.get(j);
if ( distanceLists.get(v).get(u) < minimumDistance ) {
minimumDistance = distanceLists.get(v).get(u);
closestUIndex = Optional.of(j);
}
}
}
if ( closestUIndex.isPresent() ) {
final int j = closestUIndex.get();
final int u = currentOdd.get(j);
minimumWeightMatching.addLast( new Edge(v, u, minimumDistance) );
visited.add(i); // Mark both vertices as processed
visited.add(j);
} else {
// This should not happen in the Christofides algorithm
// as the number of odd vertices is always even.
// If it does, it might indicate an issue with the input data
// or a graph where matching is not possible.
throw new AssertionError("Could not match an odd vertex in minimum weight matching: " + v);
}
}
return minimumWeightMatching;
}
// Return a list of vertices forming an Eulerian circuit using Hierholzer's algorithm
public List<Integer> eulerianCircuit(List<Edge> minimumWeightMatching) {
List<Integer> circuit = new ArrayList<Integer>();
if ( minimumWeightMatching.isEmpty() ) {
return circuit;
}
record Twin(int halfEdge, int index) {}
// Create adjacency lists for the argument 'minimumWeightMatching'
List<List<Twin>> adjacencyLists = IntStream.range(0, minimumWeightMatching.size())
.boxed().map( i -> new ArrayList<Twin>() ).collect(Collectors.toList());
IntStream.range(0, minimumWeightMatching.size()).forEach( index -> {
Edge edge = minimumWeightMatching.get(index);
adjacencyLists.get(edge.u).addLast( new Twin(edge.v, index) );
adjacencyLists.get(edge.v).addLast( new Twin(edge.u, index) );
} );
Set<Integer> edgesUsed = new HashSet<Integer>();
Stack<Integer> stack = new Stack<Integer>();
// Start from any vertex having edges.
// A suitable vertex is guaranteed to exist if 'minimumSpanningTree' is not empty
int currentVertex = minimumWeightMatching.getFirst().u;
stack.push(currentVertex);
while ( ! stack.isEmpty() ) {
currentVertex = stack.peek();
boolean foundEdge = false;
// Attempt to find an unused edge from the current vertex
while ( ! adjacencyLists.get(currentVertex).isEmpty() ) {
Twin twin = adjacencyLists.get(currentVertex).removeLast();
if ( ! edgesUsed.contains(twin.index) ) {
edgesUsed.add(twin.index);
stack.push(twin.halfEdge);
foundEdge = true;
break; // Move to the next vertex which is 'twin.halfEdge'
}
}
// If no unused edge was found from 'currentVertex',
// either the adjacency list was empty or all its edges have been used
if ( ! foundEdge ) {
circuit.addLast(stack.pop()); // Backtrack
}
}
Collections.reverse(circuit); // The circuit has been constructed in reverse order
return circuit;
}
public Result hamiltonianCircuit(List<Integer> eulerianCircuit) {
// Create a Hamiltonian circuit by removing any repeated visits to the same vertex
List<Integer> christofidesPath = new ArrayList<Integer>();
double length = 0.0;
Set<Integer> visited = new HashSet<Integer>();
int current = eulerianCircuit.getFirst();
christofidesPath.addLast(current);
visited.add(current);
for ( int vertex : eulerianCircuit ) {
if ( ! visited.contains(vertex) ) {
christofidesPath.addLast(vertex);
visited.add(vertex);
length += distanceLists.get(current).get(vertex); // Add distance from previous vertex in path
current = vertex; // Update current vertex in path
}
}
// Add the edge returning to the start vertex
length += distanceLists.get(current).get(christofidesPath.getFirst());
christofidesPath.addLast(christofidesPath.getFirst()); // Complete the cycle
return new Result(christofidesPath, length);
}
public void display() {
System.out.println("Graph: {");
IntStream.range(0, pointCount).forEach( u -> {
System.out.print(String.format("%3d%s", u, ": {"));
IntStream.range(0, pointCount).forEach( v -> {
if ( u != v ) {
if ( ! ( u == 0 && v == 1 ) && v > 0 ) {
System.out.print(", ");
}
System.out.print(String.format("%d%s%.2f", v, ": ", distanceLists.get(u).get(v)));
}
} );
System.out.println("}" + ( u == pointCount - 1 ? "" : "," ));
} );
System.out.println("}" + System.lineSeparator());
}
private static final class UnionFind {
public UnionFind(int number) {
parents = IntStream.range(0, number).boxed().collect(Collectors.toList());
ranks = new ArrayList<Integer>(Collections.nCopies(number, 0));
}
public int find(int n) {
if ( parents.get(n) == n ) {
return n;
}
// Path compression
parents.set(n, find(parents.get(n)));
return parents.get(n);
}
public boolean unite(int i, int j) {
final int rootI = find(i);
final int rootJ = find(j);
if ( rootI != rootJ ) {
switch ( Integer.compare(rootI, rootJ) ) {
case -1 -> parents.set(rootI, rootJ);
case +1 -> parents.set(rootJ, rootI);
case 0 -> { parents.set(rootJ, rootI); ranks.set(rootI, ranks.get(rootI) + 1); }
}
return true;
}
return false;
}
private List<Integer> parents;
private List<Integer> ranks;
}
private List<List<Double>> distanceLists;
private final int pointCount;
}
private record Result(List<Integer> path, double length) {}
}
- Output:
Graph: {
// omitted for brevity
}
Minimum spanning tree: [(14, 16, 13.04), (54, 82, 23.35), (51, 77, 40.00), (5, 48, 47.07), (27, 92, 47.68), (6, 8, 50.25), (24, 80, 51.48), (15, 87, 58.14), (20, 73, 66.01), (77, 95, 68.15), (4, 36, 69.23), (15, 21, 70.18), (39, 63, 80.45), (26, 85, 81.15), (40, 70, 82.68), (30, 79, 87.21), (58, 73, 89.00), (32, 75, 92.20), (6, 56, 93.05), (23, 37, 96.05), (90, 97, 99.41), (12, 32, 100.66), (39, 53, 103.59), (22, 97, 107.94), (10, 14, 113.07), (64, 65, 114.77), (70, 99, 121.43), (21, 93, 133.65), (11, 26, 136.49), (2, 45, 139.26), (9, 83, 140.04), (34, 85, 141.06), (33, 82, 145.96), (50, 86, 146.37), (3, 64, 147.18), (28, 45, 148.41), (59, 61, 151.05), (20, 71, 151.71), (67, 84, 153.40), (11, 19, 153.58), (5, 62, 157.26), (43, 49, 161.76), (24, 60, 164.71), (2, 42, 171.03), (44, 90, 179.04), (49, 72, 179.82), (41, 88, 180.42), (17, 23, 181.34), (21, 69, 186.13), (50, 60, 186.27), (30, 88, 187.77), (59, 76, 191.02), (8, 86, 195.49), (10, 31, 201.85), (19, 56, 203.75), (7, 91, 206.31), (48, 89, 208.12), (31, 90, 215.24), (12, 94, 215.64), (55, 96, 216.11), (2, 13, 221.30), (52, 78, 223.01), (35, 37, 227.76), (17, 78, 229.65), (52, 87, 229.75), (11, 54, 229.92), (18, 74, 234.08), (74, 96, 236.54), (4, 77, 242.17), (81, 94, 245.31), (7, 41, 247.78), (47, 99, 251.75), (1, 43, 256.56), (29, 38, 261.24), (67, 72, 266.72), (13, 70, 267.55), (9, 71, 269.65), (18, 52, 279.87), (55, 79, 280.80), (25, 64, 283.34), (38, 95, 283.64), (0, 62, 288.09), (68, 72, 292.88), (63, 68, 297.38), (35, 98, 299.71), (9, 89, 300.03), (28, 33, 300.50), (27, 66, 300.85), (0, 91, 302.55), (27, 57, 302.60), (68, 80, 303.12), (61, 85, 307.65), (3, 96, 324.23), (57, 60, 325.04), (10, 73, 328.69), (22, 59, 345.32), (31, 46, 361.33), (38, 84, 371.55), (32, 36, 407.77)]
Odd vertices in minimum spanning tree: [1, 2, 9, 10, 11, 16, 21, 25, 27, 29, 31, 32, 34, 38, 40, 42, 44, 46, 47, 51, 52, 53, 58, 59, 60, 64, 65, 66, 68, 69, 70, 72, 73, 75, 76, 77, 81, 83, 85, 90, 92, 93, 96, 98]
Minimum weight matching: [(14, 16, 13.04), (54, 82, 23.35), (51, 77, 40.00), (5, 48, 47.07), (27, 92, 47.68), (6, 8, 50.25), (24, 80, 51.48), (15, 87, 58.14), (20, 73, 66.01), (77, 95, 68.15), (4, 36, 69.23), (15, 21, 70.18), (39, 63, 80.45), (26, 85, 81.15), (40, 70, 82.68), (30, 79, 87.21), (58, 73, 89.00), (32, 75, 92.20), (6, 56, 93.05), (23, 37, 96.05), (90, 97, 99.41), (12, 32, 100.66), (39, 53, 103.59), (22, 97, 107.94), (10, 14, 113.07), (64, 65, 114.77), (70, 99, 121.43), (21, 93, 133.65), (11, 26, 136.49), (2, 45, 139.26), (9, 83, 140.04), (34, 85, 141.06), (33, 82, 145.96), (50, 86, 146.37), (3, 64, 147.18), (28, 45, 148.41), (59, 61, 151.05), (20, 71, 151.71), (67, 84, 153.40), (11, 19, 153.58), (5, 62, 157.26), (43, 49, 161.76), (24, 60, 164.71), (2, 42, 171.03), (44, 90, 179.04), (49, 72, 179.82), (41, 88, 180.42), (17, 23, 181.34), (21, 69, 186.13), (50, 60, 186.27), (30, 88, 187.77), (59, 76, 191.02), (8, 86, 195.49), (10, 31, 201.85), (19, 56, 203.75), (7, 91, 206.31), (48, 89, 208.12), (31, 90, 215.24), (12, 94, 215.64), (55, 96, 216.11), (2, 13, 221.30), (52, 78, 223.01), (35, 37, 227.76), (17, 78, 229.65), (52, 87, 229.75), (11, 54, 229.92), (18, 74, 234.08), (74, 96, 236.54), (4, 77, 242.17), (81, 94, 245.31), (7, 41, 247.78), (47, 99, 251.75), (1, 43, 256.56), (29, 38, 261.24), (67, 72, 266.72), (13, 70, 267.55), (9, 71, 269.65), (18, 52, 279.87), (55, 79, 280.80), (25, 64, 283.34), (38, 95, 283.64), (0, 62, 288.09), (68, 72, 292.88), (63, 68, 297.38), (35, 98, 299.71), (9, 89, 300.03), (28, 33, 300.50), (27, 66, 300.85), (0, 91, 302.55), (27, 57, 302.60), (68, 80, 303.12), (61, 85, 307.65), (3, 96, 324.23), (57, 60, 325.04), (10, 73, 328.69), (22, 59, 345.32), (31, 46, 361.33), (38, 84, 371.55), (32, 36, 407.77), (38, 29, 261.24), (77, 51, 40.00), (66, 27, 300.85), (21, 93, 133.65), (92, 46, 533.03), (90, 44, 179.04), (60, 68, 435.39), (10, 16, 115.26), (75, 32, 92.20), (85, 34, 141.06), (69, 52, 387.62), (59, 76, 191.02), (1, 53, 331.30), (96, 64, 466.62), (81, 72, 691.40), (2, 42, 171.03), (58, 73, 89.00), (31, 9, 795.62), (11, 70, 1183.26), (83, 98, 593.33), (47, 40, 432.65), (25, 65, 323.20)]
Eulerian circuit: [14, 10, 73, 58, 73, 20, 71, 9, 89, 48, 5, 62, 0, 91, 7, 41, 88, 30, 79, 55, 96, 64, 25, 65, 64, 3, 96, 74, 18, 52, 69, 21, 93, 21, 15, 87, 52, 78, 17, 23, 37, 35, 98, 83, 9, 31, 46, 92, 27, 66, 27, 57, 60, 68, 63, 39, 53, 1, 43, 49, 72, 81, 94, 12, 32, 75, 32, 36, 4, 77, 51, 77, 95, 38, 29, 38, 84, 67, 72, 68, 80, 24, 60, 50, 86, 8, 6, 56, 19, 11, 70, 99, 47, 40, 70, 13, 2, 42, 2, 45, 28, 33, 82, 54, 11, 26, 85, 34, 85, 61, 59, 76, 59, 22, 97, 90, 44, 90, 31, 10, 16, 14]
Result path: [14, 10, 73, 58, 20, 71, 9, 89, 48, 5, 62, 0, 91, 7, 41, 88, 30, 79, 55, 96, 64, 25, 65, 3, 74, 18, 52, 69, 21, 93, 15, 87, 78, 17, 23, 37, 35, 98, 83, 31, 46, 92, 27, 66, 57, 60, 68, 63, 39, 53, 1, 43, 49, 72, 81, 94, 12, 32, 75, 36, 4, 77, 51, 95, 38, 29, 84, 67, 80, 24, 50, 86, 8, 6, 56, 19, 11, 70, 99, 47, 40, 13, 2, 42, 45, 28, 33, 82, 54, 26, 85, 34, 61, 59, 76, 22, 97, 90, 44, 16, 14]
Length of the result path: 24502.99
JavaScript
// --- Helper Structs (using Objects) ---
// Point: { x: number, y: number, id: number }
// Edge: { u: number, v: number, weight: number, id?: number } // id added for Eulerian tour
// --- Helper Functions ---
function printContainer(container, name) {
console.log(`${name}: [${container.join(', ')}]`);
}
function printEdges(edges, name) {
const edgeStrings = edges.map(e => `(${e.u}, ${e.v}, ${e.weight.toFixed(2)})`);
console.log(`${name}: [${edgeStrings.join(', ')}]`);
}
function printGraph(graph, name) {
console.log(`${name}: {`);
const n = graph.length;
for (let i = 0; i < n; i++) {
const entries = [];
for (let j = 0; j < n; j++) {
if (i !== j) {
entries.push(`${j}: ${graph[i][j].toFixed(2)}`);
}
}
console.log(` ${i}: {${entries.join(', ')}}${i === n - 1 ? '' : ','}`);
}
console.log(`}`);
}
// --- Euclidean Distance ---
function getLength(p1, p2) {
const dx = p1.x - p2.x;
const dy = p1.y - p2.y;
return Math.sqrt(dx * dx + dy * dy);
}
// --- Build Complete Graph (Adjacency Matrix) ---
function buildGraph(data) {
const n = data.length;
// Initialize n x n array with 0s
const graph = Array(n).fill(0).map(() => Array(n).fill(0.0));
for (let i = 0; i < n; i++) {
for (let j = i + 1; j < n; j++) { // Only calculate upper triangle
const dist = getLength(data[i], data[j]);
graph[i][j] = dist;
graph[j][i] = dist; // Symmetric graph
}
}
return graph;
}
// --- Union-Find Data Structure ---
class UnionFind {
constructor(n) {
// Initialize parent array: each node is its own parent
this.parent = Array.from({ length: n }, (_, i) => i);
// Initialize rank (or size) array for optimization
this.rank = Array(n).fill(0);
}
find(i) {
if (this.parent[i] === i) {
return i;
}
// Path compression: point directly to the root
this.parent[i] = this.find(this.parent[i]);
return this.parent[i];
}
unite(i, j) {
let rootI = this.find(i);
let rootJ = this.find(j);
if (rootI !== rootJ) {
// Union by rank: attach smaller rank tree under larger rank tree
if (this.rank[rootI] < this.rank[rootJ]) {
this.parent[rootI] = rootJ;
} else if (this.rank[rootI] > this.rank[rootJ]) {
this.parent[rootJ] = rootI;
} else {
// Ranks are equal, choose one as parent and increment its rank
this.parent[rootJ] = rootI;
this.rank[rootI]++;
}
return true; // Successfully united
}
return false; // Already in the same set
}
}
// --- Minimum Spanning Tree (Kruskal's Algorithm) ---
function minimumSpanningTree(graph) {
const n = graph.length;
if (n === 0) return [];
const edges = [];
for (let i = 0; i < n; i++) {
for (let j = i + 1; j < n; j++) { // Avoid duplicates and self-loops
edges.push({ u: i, v: j, weight: graph[i][j] });
}
}
// Sort edges by weight in ascending order
edges.sort((a, b) => a.weight - b.weight);
const mst = [];
const uf = new UnionFind(n);
let edgesCount = 0;
for (const edge of edges) {
if (uf.unite(edge.u, edge.v)) { // If uniting forms a new connection
mst.push(edge);
edgesCount++;
if (edgesCount === n - 1) { // Optimization: MST has n-1 edges
break;
}
}
}
return mst;
}
// --- Find Vertices with Odd Degree in MST ---
function findOddVertexes(mst, n) {
const degree = Array(n).fill(0);
for (const edge of mst) {
degree[edge.u]++;
degree[edge.v]++;
}
const oddVertices = [];
for (let i = 0; i < n; i++) {
if (degree[i] % 2 !== 0) {
oddVertices.push(i);
}
}
return oddVertices;
}
// --- Minimum Weight Matching (Greedy Heuristic) ---
// Fisher-Yates (Knuth) Shuffle algorithm
function shuffleArray(array) {
for (let i = array.length - 1; i > 0; i--) {
const j = Math.floor(Math.random() * (i + 1));
[array[i], array[j]] = [array[j], array[i]]; // ES6 swap
}
}
// Note: This modifies the mst array by adding matching edges.
function minimumWeightMatching(mst, graph, oddVertices) {
// Shuffle for randomness (like the Python version)
shuffleArray(oddVertices);
const matched = new Set(); // Keep track of vertices already matched in this phase
for (let i = 0; i < oddVertices.length; i++) {
const v = oddVertices[i];
if (matched.has(v)) continue; // Skip if already matched
let minLength = Infinity;
let closestU = -1;
// Find the closest *unmatched* odd vertex
for (let j = i + 1; j < oddVertices.length; j++) {
const u = oddVertices[j];
if (!matched.has(u)) { // Check if 'u' is available
if (graph[v][u] < minLength) {
minLength = graph[v][u];
closestU = u;
}
}
}
if (closestU !== -1) {
// Add the matching edge to the MST list (now a multigraph)
mst.push({ u: v, v: closestU, weight: minLength });
matched.add(v);
matched.add(closestU); // Mark both as matched
}
// Christofides guarantees an even number of odd-degree vertices,
// so every vertex *should* find a match in a perfect matching scenario.
// The greedy approach might leave some unmatched if not careful, but this loop structure should work.
}
// No return value needed as mst is modified directly by reference
}
// --- Find Eulerian Tour (Hierholzer's Algorithm) ---
function findEulerianTour(matchedMST, n) {
if (matchedMST.length === 0) return [];
// Assign unique IDs to edges for tracking in the multigraph
// (essential because multiple edges can exist between nodes)
matchedMST.forEach((edge, index) => edge.id = index);
// Build adjacency list: adj[u] = [{ neighbor: v, edgeId: id }, ...]
const adj = Array(n).fill(0).map(() => []);
const edgeUsed = new Set(); // Store used edge IDs
for (const edge of matchedMST) {
adj[edge.u].push({ neighbor: edge.v, edgeId: edge.id });
adj[edge.v].push({ neighbor: edge.u, edgeId: edge.id });
}
const tour = [];
const stack = [];
// Start at any vertex with edges (guaranteed to exist if matchedMST is not empty)
const startNode = matchedMST[0].u;
stack.push(startNode);
let currentNode = startNode;
while (stack.length > 0) {
currentNode = stack[stack.length - 1]; // Peek top of stack
let foundUnusedEdge = false;
// Check outgoing edges from currentNode
while (adj[currentNode].length > 0) {
const edgeInfo = adj[currentNode][adj[currentNode].length - 1]; // Look at last edge (efficient removal)
if (!edgeUsed.has(edgeInfo.edgeId)) {
// Found an unused edge
edgeUsed.add(edgeInfo.edgeId); // Mark edge as used
stack.push(edgeInfo.neighbor); // Push neighbor onto stack
adj[currentNode].pop(); // Remove edge from current node's list
currentNode = edgeInfo.neighbor; // Move to the neighbor
foundUnusedEdge = true;
break; // Break inner loop to process the new currentNode
} else {
// This edge was already used (by traversal from the other side), remove it
adj[currentNode].pop();
}
}
if (!foundUnusedEdge) {
// If no unused edges from currentNode, backtrack
tour.push(stack.pop());
}
}
// The tour is constructed in reverse order by Hierholzer's
return tour.reverse();
}
// --- Main TSP Function (Christofides Approximation) ---
function tsp(data) {
const n = data.length;
if (n === 0) {
return { length: 0, path: [] };
}
if (n === 1) {
// If data points have explicit IDs use them, otherwise use index 0
const startId = data[0].id !== undefined ? data[0].id : 0;
return { length: 0, path: [startId] };
}
// Assign IDs if they don't exist, assuming order corresponds to 0..n-1
data.forEach((p, i) => { if (p.id === undefined) p.id = i; });
console.log("Building graph...");
const G = buildGraph(data);
// printGraph(G, "Graph"); // Often too large to print meaningfully
console.log("Finding Minimum Spanning Tree...");
const MSTree = minimumSpanningTree(G);
printEdges(MSTree, "MSTree");
console.log("Finding odd degree vertices...");
const odd_vertexes = findOddVertexes(MSTree, n);
printContainer(odd_vertexes, "Odd vertexes in MSTree");
console.log("Finding Minimum Weight Matching (greedy)...");
// Create a new array containing MST edges to avoid modifying the original MSTree variable directly
// The matching edges will be added to this new array.
const multigraphEdges = [...MSTree];
// Pass a copy of odd_vertexes as it might be shuffled in place
minimumWeightMatching(multigraphEdges, G, [...odd_vertexes]);
printEdges(multigraphEdges, "Minimum weight matching (MST + Matching Edges)");
console.log("Finding Eulerian Tour...");
const eulerian_tour = findEulerianTour(multigraphEdges, n);
printContainer(eulerian_tour, "Eulerian tour");
// --- Create Hamiltonian Circuit by Skipping Visited Nodes ---
if (eulerian_tour.length === 0 && n > 0) {
console.error("Error: Eulerian tour could not be found.");
return { length: -1, path: [] }; // Indicate error
}
console.log("Creating Hamiltonian path (shortcutting)...");
const path = [];
let length = 0.0;
const visited = new Set(); // Use Set for efficient O(1) average time complexity checks
let currentPathNode = -1; // Track the last node added to the *Hamiltonian* path
for (const node of eulerian_tour) {
if (!visited.has(node)) {
path.push(node);
visited.add(node);
if (currentPathNode !== -1) { // Add edge length from the previous node *in the path*
length += G[currentPathNode][node];
}
currentPathNode = node; // Update the last node added to the path
}
}
// Add the edge back to the start to complete the cycle
if (path.length > 0) {
length += G[currentPathNode][path[0]]; // Edge from last node in path to first node
path.push(path[0]); // Add the starting node again to show the cycle
}
printContainer(path, "Result path");
console.log(`Result length of the path: ${length.toFixed(2)}`);
return { length: length, path: path };
}
// --- Input Data ---
const raw_data = [
[1380, 939], [2848, 96], [3510, 1671], [457, 334], [3888, 666], [984, 965], [2721, 1482], [1286, 525],
[2716, 1432], [738, 1325], [1251, 1832], [2728, 1698], [3815, 169], [3683, 1533], [1247, 1945], [123, 862],
[1234, 1946], [252, 1240], [611, 673], [2576, 1676], [928, 1700], [53, 857], [1807, 1711], [274, 1420],
[2574, 946], [178, 24], [2678, 1825], [1795, 962], [3384, 1498], [3520, 1079], [1256, 61], [1424, 1728],
[3913, 192], [3085, 1528], [2573, 1969], [463, 1670], [3875, 598], [298, 1513], [3479, 821], [2542, 236],
[3955, 1743], [1323, 280], [3447, 1830], [2936, 337], [1621, 1830], [3373, 1646], [1393, 1368],
[3874, 1318], [938, 955], [3022, 474], [2482, 1183], [3854, 923], [376, 825], [2519, 135], [2945, 1622],
[953, 268], [2628, 1479], [2097, 981], [890, 1846], [2139, 1806], [2421, 1007], [2290, 1810], [1115, 1052],
[2588, 302], [327, 265], [241, 341], [1917, 687], [2991, 792], [2573, 599], [19, 674], [3911, 1673],
[872, 1559], [2863, 558], [929, 1766], [839, 620], [3893, 102], [2178, 1619], [3822, 899], [378, 1048],
[1178, 100], [2599, 901], [3416, 143], [2961, 1605], [611, 1384], [3113, 885], [2597, 1830], [2586, 1286],
[161, 906], [1429, 134], [742, 1025], [1625, 1651], [1187, 706], [1787, 1009], [22, 987], [3640, 43],
[3756, 882], [776, 392], [1724, 1642], [198, 1810], [3950, 1558]
];
// Convert raw data to point objects, using index as ID
const data_points = raw_data.map((coords, index) => ({
x: coords[0],
y: coords[1],
id: index // Use 0-based index as the vertex ID
}));
// --- Run TSP ---
tsp(data_points);
- Output:
Building graph... Finding Minimum Spanning Tree... MSTree: [(14, 16, 13.04), (54, 82, 23.35), (51, 77, 40.00), (5, 48, 47.07), (27, 92, 47.68), (6, 8, 50.25), (24, 80, 51.48), (15, 87, 58.14), (20, 73, 66.01), (77, 95, 68.15), (4, 36, 69.23), (15, 21, 70.18), (39, 63, 80.45), (26, 85, 81.15), (40, 70, 82.68), (30, 79, 87.21), (58, 73, 89.00), (32, 75, 92.20), (6, 56, 93.05), (23, 37, 96.05), (90, 97, 99.41), (12, 32, 100.66), (39, 53, 103.59), (22, 97, 107.94), (10, 14, 113.07), (64, 65, 114.77), (70, 99, 121.43), (21, 93, 133.65), (11, 26, 136.49), (2, 45, 139.26), (9, 83, 140.04), (34, 85, 141.06), (33, 82, 145.96), (50, 86, 146.37), (3, 64, 147.18), (28, 45, 148.41), (59, 61, 151.05), (20, 71, 151.71), (67, 84, 153.40), (11, 19, 153.58), (5, 62, 157.26), (43, 49, 161.76), (24, 60, 164.71), (2, 42, 171.03), (44, 90, 179.04), (49, 72, 179.82), (41, 88, 180.42), (17, 23, 181.34), (21, 69, 186.13), (50, 60, 186.27), (30, 88, 187.77), (59, 76, 191.02), (8, 86, 195.49), (10, 31, 201.85), (19, 56, 203.75), (7, 91, 206.31), (48, 89, 208.12), (31, 90, 215.24), (12, 94, 215.64), (55, 96, 216.11), (2, 13, 221.30), (52, 78, 223.01), (35, 37, 227.76), (17, 78, 229.65), (52, 87, 229.75), (11, 54, 229.92), (18, 74, 234.08), (74, 96, 236.54), (4, 77, 242.17), (81, 94, 245.31), (7, 41, 247.78), (47, 99, 251.75), (1, 43, 256.56), (29, 38, 261.24), (67, 72, 266.72), (13, 70, 267.55), (9, 71, 269.65), (18, 52, 279.87), (55, 79, 280.80), (25, 64, 283.34), (38, 95, 283.64), (0, 62, 288.09), (68, 72, 292.88), (63, 68, 297.38), (35, 98, 299.71), (9, 89, 300.03), (28, 33, 300.50), (27, 66, 300.85), (0, 91, 302.55), (27, 57, 302.60), (68, 80, 303.12), (61, 85, 307.65), (3, 96, 324.23), (57, 60, 325.04), (10, 73, 328.69), (22, 59, 345.32), (31, 46, 361.33), (38, 84, 371.55), (32, 36, 407.77)] Finding odd degree vertices... Odd vertexes in MSTree: [1, 2, 9, 10, 11, 16, 21, 25, 27, 29, 31, 32, 34, 38, 40, 42, 44, 46, 47, 51, 52, 53, 58, 59, 60, 64, 65, 66, 68, 69, 70, 72, 73, 75, 76, 77, 81, 83, 85, 90, 92, 93, 96, 98] Finding Minimum Weight Matching (greedy)... Minimum weight matching (MST + Matching Edges): [(14, 16, 13.04), (54, 82, 23.35), (51, 77, 40.00), (5, 48, 47.07), (27, 92, 47.68), (6, 8, 50.25), (24, 80, 51.48), (15, 87, 58.14), (20, 73, 66.01), (77, 95, 68.15), (4, 36, 69.23), (15, 21, 70.18), (39, 63, 80.45), (26, 85, 81.15), (40, 70, 82.68), (30, 79, 87.21), (58, 73, 89.00), (32, 75, 92.20), (6, 56, 93.05), (23, 37, 96.05), (90, 97, 99.41), (12, 32, 100.66), (39, 53, 103.59), (22, 97, 107.94), (10, 14, 113.07), (64, 65, 114.77), (70, 99, 121.43), (21, 93, 133.65), (11, 26, 136.49), (2, 45, 139.26), (9, 83, 140.04), (34, 85, 141.06), (33, 82, 145.96), (50, 86, 146.37), (3, 64, 147.18), (28, 45, 148.41), (59, 61, 151.05), (20, 71, 151.71), (67, 84, 153.40), (11, 19, 153.58), (5, 62, 157.26), (43, 49, 161.76), (24, 60, 164.71), (2, 42, 171.03), (44, 90, 179.04), (49, 72, 179.82), (41, 88, 180.42), (17, 23, 181.34), (21, 69, 186.13), (50, 60, 186.27), (30, 88, 187.77), (59, 76, 191.02), (8, 86, 195.49), (10, 31, 201.85), (19, 56, 203.75), (7, 91, 206.31), (48, 89, 208.12), (31, 90, 215.24), (12, 94, 215.64), (55, 96, 216.11), (2, 13, 221.30), (52, 78, 223.01), (35, 37, 227.76), (17, 78, 229.65), (52, 87, 229.75), (11, 54, 229.92), (18, 74, 234.08), (74, 96, 236.54), (4, 77, 242.17), (81, 94, 245.31), (7, 41, 247.78), (47, 99, 251.75), (1, 43, 256.56), (29, 38, 261.24), (67, 72, 266.72), (13, 70, 267.55), (9, 71, 269.65), (18, 52, 279.87), (55, 79, 280.80), (25, 64, 283.34), (38, 95, 283.64), (0, 62, 288.09), (68, 72, 292.88), (63, 68, 297.38), (35, 98, 299.71), (9, 89, 300.03), (28, 33, 300.50), (27, 66, 300.85), (0, 91, 302.55), (27, 57, 302.60), (68, 80, 303.12), (61, 85, 307.65), (3, 96, 324.23), (57, 60, 325.04), (10, 73, 328.69), (22, 59, 345.32), (31, 46, 361.33), (38, 84, 371.55), (32, 36, 407.77), (72, 68, 292.88), (11, 85, 185.97), (2, 42, 171.03), (70, 40, 82.68), (96, 64, 466.62), (53, 1, 331.30), (21, 93, 133.65), (75, 32, 92.20), (83, 9, 140.04), (77, 51, 40.00), (16, 10, 115.26), (73, 58, 89.00), (25, 65, 323.20), (46, 31, 361.33), (29, 38, 261.24), (92, 27, 47.68), (47, 81, 1261.11), (59, 76, 191.02), (69, 52, 387.62), (66, 60, 597.01), (44, 90, 179.04), (98, 34, 2380.32)] Finding Eulerian Tour... Eulerian tour: [14, 10, 73, 58, 73, 20, 71, 9, 83, 9, 89, 48, 5, 62, 0, 91, 7, 41, 88, 30, 79, 55, 96, 64, 25, 65, 64, 3, 96, 74, 18, 52, 69, 21, 93, 21, 15, 87, 52, 78, 17, 23, 37, 35, 98, 34, 85, 11, 54, 82, 33, 28, 45, 2, 42, 2, 13, 70, 40, 70, 99, 47, 81, 94, 12, 32, 75, 32, 36, 4, 77, 51, 77, 95, 38, 29, 38, 84, 67, 72, 68, 63, 39, 53, 1, 43, 49, 72, 68, 80, 24, 60, 66, 27, 92, 27, 57, 60, 50, 86, 8, 6, 56, 19, 11, 26, 85, 61, 59, 76, 59, 22, 97, 90, 44, 90, 31, 46, 31, 10, 16, 14] Creating Hamiltonian path (shortcutting)... Result path: [14, 10, 73, 58, 20, 71, 9, 83, 89, 48, 5, 62, 0, 91, 7, 41, 88, 30, 79, 55, 96, 64, 25, 65, 3, 74, 18, 52, 69, 21, 93, 15, 87, 78, 17, 23, 37, 35, 98, 34, 85, 11, 54, 82, 33, 28, 45, 2, 42, 13, 70, 40, 99, 47, 81, 94, 12, 32, 75, 36, 4, 77, 51, 95, 38, 29, 84, 67, 72, 68, 63, 39, 53, 1, 43, 49, 80, 24, 60, 66, 27, 92, 57, 50, 86, 8, 6, 56, 19, 26, 61, 59, 76, 22, 97, 90, 44, 31, 46, 16, 14] Result length of the path: 25448.96
Rust
use std::collections::HashSet;
use std::f64;
use rand::seq::SliceRandom;
use rand::thread_rng;
// --- Structs ---
#[derive(Debug, Clone, Copy, PartialEq)]
struct Point {
id: usize, // Original index or assigned ID
x: f64,
y: f64,
}
#[derive(Debug, Clone, Copy)]
struct Edge {
u: usize,
v: usize,
weight: f64,
}
// Implement PartialEq, Eq, PartialOrd, Ord for sorting Edges by weight
impl PartialEq for Edge {
fn eq(&self, other: &Self) -> bool {
self.weight == other.weight
}
}
impl Eq for Edge {}
impl PartialOrd for Edge {
fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
self.weight.partial_cmp(&other.weight)
}
}
impl Ord for Edge {
fn cmp(&self, other: &Self) -> std::cmp::Ordering {
self.partial_cmp(other).unwrap_or(std::cmp::Ordering::Equal) // Panic on NaN, ok for this problem
}
}
// --- Helper Functions ---
fn print_usize_vec(vec: &[usize], name: &str) {
let content: Vec<String> = vec.iter().map(|&x| x.to_string()).collect();
println!("{}: [{}]", name, content.join(", "));
}
fn print_edges(edges: &[Edge], name: &str) {
let edge_strings: Vec<String> = edges
.iter()
.map(|e| format!("({}, {}, {:.2})", e.u, e.v, e.weight))
.collect();
println!("{}: [{}]", name, edge_strings.join(", "));
}
fn print_graph(graph: &[Vec<f64>], name: &str) {
println!("{}: {{", name);
let n = graph.len();
for i in 0..n {
let entries: Vec<String> = (0..n)
.filter(|&j| i != j)
.map(|j| format!("{}: {:.2}", j, graph[i][j]))
.collect();
println!(" {}: {{{}}}{}", i, entries.join(", "), if i == n - 1 { "" } else { "," });
}
println!("}}");
}
// --- Euclidean Distance ---
fn get_length(p1: Point, p2: Point) -> f64 {
let dx = p1.x - p2.x;
let dy = p1.y - p2.y;
(dx * dx + dy * dy).sqrt()
}
// --- Build Complete Graph (Adjacency Matrix) ---
fn build_graph(data: &[Point]) -> Vec<Vec<f64>> {
let n = data.len();
let mut graph = vec![vec![0.0; n]; n];
for i in 0..n {
for j in i + 1..n { // Only calculate upper triangle
let dist = get_length(data[i], data[j]);
graph[i][j] = dist;
graph[j][i] = dist; // Symmetric graph
}
}
graph
}
// --- Union-Find Data Structure ---
struct UnionFind {
parent: Vec<usize>,
rank: Vec<usize>,
}
impl UnionFind {
fn new(n: usize) -> Self {
UnionFind {
parent: (0..n).collect(),
rank: vec![0; n],
}
}
fn find(&mut self, i: usize) -> usize {
if self.parent[i] == i {
i
} else {
// Path compression
self.parent[i] = self.find(self.parent[i]);
self.parent[i]
}
}
fn unite(&mut self, i: usize, j: usize) -> bool {
let root_i = self.find(i);
let root_j = self.find(j);
if root_i != root_j {
// Union by rank
if self.rank[root_i] < self.rank[root_j] {
self.parent[root_i] = root_j;
} else if self.rank[root_i] > self.rank[root_j] {
self.parent[root_j] = root_i;
} else {
self.parent[root_j] = root_i;
self.rank[root_i] += 1;
}
true
} else {
false
}
}
}
// --- Minimum Spanning Tree (Kruskal's Algorithm) ---
fn minimum_spanning_tree(graph: &[Vec<f64>]) -> Vec<Edge> {
let n = graph.len();
if n == 0 {
return vec![];
}
let mut edges = Vec::new();
for i in 0..n {
for j in i + 1..n { // Avoid duplicates and self-loops
edges.push(Edge { u: i, v: j, weight: graph[i][j] });
}
}
// Sort edges by weight
edges.sort_unstable(); // Uses the Ord implementation
let mut mst = Vec::new();
let mut uf = UnionFind::new(n);
let mut edges_count = 0;
for &edge in &edges {
if uf.unite(edge.u, edge.v) {
mst.push(edge);
edges_count += 1;
if edges_count == n - 1 { // Optimization: MST has n-1 edges
break;
}
}
}
mst
}
// --- Find Vertices with Odd Degree in MST ---
fn find_odd_vertexes(mst: &[Edge], n: usize) -> Vec<usize> {
let mut degree = vec![0; n];
for edge in mst {
degree[edge.u] += 1;
degree[edge.v] += 1;
}
(0..n).filter(|&i| degree[i] % 2 != 0).collect()
}
// --- Minimum Weight Matching (Greedy Heuristic) ---
// Note: This modifies the multigraph_edges vector by adding matching edges.
fn minimum_weight_matching(
multigraph_edges: &mut Vec<Edge>,
graph: &[Vec<f64>],
odd_vertices: &[usize],
) {
let n = graph.len();
if odd_vertices.is_empty() {
return;
}
let mut current_odd = odd_vertices.to_vec(); // Clone to allow shuffling and tracking
let mut rng = thread_rng();
current_odd.shuffle(&mut rng); // Shuffle for randomness
let mut matched = vec![false; n]; // Keep track of vertices matched in this phase
// We use a marker for visited indices in the shuffled list
let mut processed = vec![false; current_odd.len()];
for i in 0..current_odd.len() {
if processed[i] { continue; } // Skip if already processed (matched)
let v = current_odd[i];
let mut min_length = f64::INFINITY;
let mut closest_u_idx: Option<usize> = None; // Store index in current_odd
// Find the closest *unmatched* odd vertex *later* in the shuffled list
for j in i + 1..current_odd.len() {
if !processed[j] { // Check if 'u' (at current_odd[j]) is available
let u = current_odd[j];
if graph[v][u] < min_length {
min_length = graph[v][u];
closest_u_idx = Some(j);
}
}
}
if let Some(j) = closest_u_idx {
let u = current_odd[j];
// Add the matching edge to the MST list (now a multigraph)
multigraph_edges.push(Edge { u: v, v: u, weight: min_length });
// Mark both as processed in the shuffled list context
processed[i] = true;
processed[j] = true;
// Optionally mark in the global 'matched' array if needed elsewhere
matched[v] = true;
matched[u] = true;
} else {
// This *shouldn't* happen in Christofides as the number of odd vertices is always even.
// If it does, it might indicate an issue earlier or a graph where matching isn't possible?
eprintln!("Warning: Could not find match for odd vertex {} in greedy matching.", v);
}
}
}
// --- Find Eulerian Tour (Hierholzer's Algorithm) ---
fn find_eulerian_tour(multigraph_edges: &[Edge], n: usize) -> Vec<usize> {
if multigraph_edges.is_empty() {
return vec![];
}
// Build adjacency list: adj[u] = Vec<(neighbor, edge_index)>
let mut adj: Vec<Vec<(usize, usize)>> = vec![vec![]; n];
for (edge_idx, edge) in multigraph_edges.iter().enumerate() {
adj[edge.u].push((edge.v, edge_idx));
adj[edge.v].push((edge.u, edge_idx));
}
let mut edge_used = vec![false; multigraph_edges.len()];
let mut tour = Vec::new();
let mut stack = Vec::new();
// Start at any vertex with edges
let start_node = multigraph_edges[0].u;
stack.push(start_node);
while let Some(¤t_node) = stack.last() {
let mut found_edge = false;
// Try to find an unused edge using pop for efficiency
while let Some((neighbor, edge_idx)) = adj[current_node].pop() {
if !edge_used[edge_idx] {
edge_used[edge_idx] = true;
stack.push(neighbor);
found_edge = true;
break; // Move to the neighbor
}
// If edge was used, loop continues popping from adj[current_node]
}
// If no unused edge was found from current_node (adj list is empty or all remaining edges used)
if !found_edge {
tour.push(stack.pop().unwrap()); // Backtrack
}
}
// The tour is constructed in reverse order
tour.reverse();
tour
}
// --- Main TSP Function (Christofides Approximation) ---
fn tsp(data: &[Point]) -> (f64, Vec<usize>) {
let n = data.len();
if n == 0 {
return (0.0, vec![]);
}
if n == 1 {
return (0.0, vec![data[0].id]); // Use the point's ID
}
println!("Building graph...");
let graph = build_graph(data);
// print_graph(&graph, "Graph"); // Can be very large
println!("Finding Minimum Spanning Tree...");
let mst = minimum_spanning_tree(&graph);
print_edges(&mst, "MSTree");
println!("Finding odd degree vertices...");
let odd_vertices = find_odd_vertexes(&mst, n);
print_usize_vec(&odd_vertices, "Odd vertexes in MSTree");
println!("Finding Minimum Weight Matching (greedy)...");
// Clone MST edges to create the initial multigraph
let mut multigraph_edges = mst.clone();
// Add matching edges to the multigraph_edges vector
minimum_weight_matching(&mut multigraph_edges, &graph, &odd_vertices);
print_edges(&multigraph_edges, "Minimum weight matching (MST + Matching Edges)");
println!("Finding Eulerian Tour...");
let eulerian_tour = find_eulerian_tour(&multigraph_edges, n);
print_usize_vec(&eulerian_tour, "Eulerian tour");
// --- Create Hamiltonian Circuit by Skipping Visited Nodes ---
println!("Creating Hamiltonian path (shortcutting)...");
if eulerian_tour.is_empty() && n > 0 {
eprintln!("Error: Eulerian tour could not be found.");
return (-1.0, vec![]); // Indicate error
}
let mut path = Vec::new();
let mut length = 0.0;
let mut visited = HashSet::new(); // Use HashSet for O(1) average lookup
let mut last_node_in_path: Option<usize> = None;
for &node in &eulerian_tour {
// visited.insert returns true if the value was not already present
if visited.insert(node) {
if let Some(last_node) = last_node_in_path {
// Add distance from the *previous node added to the path*
length += graph[last_node][node];
}
path.push(node); // Add node to the Hamiltonian path
last_node_in_path = Some(node); // Update the last node *in the path*
}
}
// Add the edge back to the start to complete the cycle
if let (Some(last), Some(&first)) = (last_node_in_path, path.first()) {
length += graph[last][first];
path.push(first); // Add the starting node ID again to show the cycle
}
print_usize_vec(&path, "Result path");
println!("Result length of the path: {:.2}", length);
(length, path)
}
// --- Main Function ---
fn main() {
// Input data matching the Python/JS example
let raw_data: Vec<(f64, f64)> = vec![
(1380.0, 939.0), (2848.0, 96.0), (3510.0, 1671.0), (457.0, 334.0), (3888.0, 666.0),
(984.0, 965.0), (2721.0, 1482.0), (1286.0, 525.0), (2716.0, 1432.0), (738.0, 1325.0),
(1251.0, 1832.0), (2728.0, 1698.0), (3815.0, 169.0), (3683.0, 1533.0), (1247.0, 1945.0),
(123.0, 862.0), (1234.0, 1946.0), (252.0, 1240.0), (611.0, 673.0), (2576.0, 1676.0),
(928.0, 1700.0), (53.0, 857.0), (1807.0, 1711.0), (274.0, 1420.0), (2574.0, 946.0),
(178.0, 24.0), (2678.0, 1825.0), (1795.0, 962.0), (3384.0, 1498.0), (3520.0, 1079.0),
(1256.0, 61.0), (1424.0, 1728.0), (3913.0, 192.0), (3085.0, 1528.0), (2573.0, 1969.0),
(463.0, 1670.0), (3875.0, 598.0), (298.0, 1513.0), (3479.0, 821.0), (2542.0, 236.0),
(3955.0, 1743.0), (1323.0, 280.0), (3447.0, 1830.0), (2936.0, 337.0), (1621.0, 1830.0),
(3373.0, 1646.0), (1393.0, 1368.0), (3874.0, 1318.0), (938.0, 955.0), (3022.0, 474.0),
(2482.0, 1183.0), (3854.0, 923.0), (376.0, 825.0), (2519.0, 135.0), (2945.0, 1622.0),
(953.0, 268.0), (2628.0, 1479.0), (2097.0, 981.0), (890.0, 1846.0), (2139.0, 1806.0),
(2421.0, 1007.0), (2290.0, 1810.0), (1115.0, 1052.0), (2588.0, 302.0), (327.0, 265.0),
(241.0, 341.0), (1917.0, 687.0), (2991.0, 792.0), (2573.0, 599.0), (19.0, 674.0),
(3911.0, 1673.0), (872.0, 1559.0), (2863.0, 558.0), (929.0, 1766.0), (839.0, 620.0),
(3893.0, 102.0), (2178.0, 1619.0), (3822.0, 899.0), (378.0, 1048.0), (1178.0, 100.0),
(2599.0, 901.0), (3416.0, 143.0), (2961.0, 1605.0), (611.0, 1384.0), (3113.0, 885.0),
(2597.0, 1830.0), (2586.0, 1286.0), (161.0, 906.0), (1429.0, 134.0), (742.0, 1025.0),
(1625.0, 1651.0), (1187.0, 706.0), (1787.0, 1009.0), (22.0, 987.0), (3640.0, 43.0),
(3756.0, 882.0), (776.0, 392.0), (1724.0, 1642.0), (198.0, 1810.0), (3950.0, 1558.0),
];
// Convert raw data to Point structs, using index as ID
let data_points: Vec<Point> = raw_data
.into_iter()
.enumerate() // Get index along with coordinates
.map(|(i, (x, y))| Point { id: i, x, y })
.collect();
// --- Run TSP ---
let (_length, _path) = tsp(&data_points);
// Result is already printed within tsp function
}
- Output:
Building graph... Finding Minimum Spanning Tree... MSTree: [(14, 16, 13.04), (54, 82, 23.35), (51, 77, 40.00), (5, 48, 47.07), (27, 92, 47.68), (6, 8, 50.25), (24, 80, 51.48), (15, 87, 58.14), (20, 73, 66.01), (77, 95, 68.15), (4, 36, 69.23), (15, 21, 70.18), (39, 63, 80.45), (26, 85, 81.15), (40, 70, 82.68), (30, 79, 87.21), (58, 73, 89.00), (32, 75, 92.20), (6, 56, 93.05), (23, 37, 96.05), (90, 97, 99.41), (12, 32, 100.66), (39, 53, 103.59), (22, 97, 107.94), (10, 14, 113.07), (64, 65, 114.77), (70, 99, 121.43), (21, 93, 133.65), (11, 26, 136.49), (2, 45, 139.26), (9, 83, 140.04), (34, 85, 141.06), (33, 82, 145.96), (50, 86, 146.37), (3, 64, 147.18), (28, 45, 148.41), (59, 61, 151.05), (20, 71, 151.71), (67, 84, 153.40), (11, 19, 153.58), (5, 62, 157.26), (43, 49, 161.76), (24, 60, 164.71), (2, 42, 171.03), (44, 90, 179.04), (49, 72, 179.82), (41, 88, 180.42), (17, 23, 181.34), (21, 69, 186.13), (50, 60, 186.27), (30, 88, 187.77), (59, 76, 191.02), (8, 86, 195.49), (10, 31, 201.85), (19, 56, 203.75), (7, 91, 206.31), (48, 89, 208.12), (31, 90, 215.24), (12, 94, 215.64), (55, 96, 216.11), (2, 13, 221.30), (52, 78, 223.01), (35, 37, 227.76), (17, 78, 229.65), (52, 87, 229.75), (11, 54, 229.92), (18, 74, 234.08), (74, 96, 236.54), (4, 77, 242.17), (81, 94, 245.31), (7, 41, 247.78), (47, 99, 251.75), (1, 43, 256.56), (29, 38, 261.24), (67, 72, 266.72), (13, 70, 267.55), (9, 71, 269.65), (18, 52, 279.87), (55, 79, 280.80), (25, 64, 283.34), (38, 95, 283.64), (0, 62, 288.09), (68, 72, 292.88), (63, 68, 297.38), (35, 98, 299.71), (9, 89, 300.03), (28, 33, 300.50), (27, 66, 300.85), (0, 91, 302.55), (27, 57, 302.60), (68, 80, 303.12), (61, 85, 307.65), (3, 96, 324.23), (57, 60, 325.04), (10, 73, 328.69), (22, 59, 345.32), (31, 46, 361.33), (38, 84, 371.55), (32, 36, 407.77)] Finding odd degree vertices... Odd vertexes in MSTree: [1, 2, 9, 10, 11, 16, 21, 25, 27, 29, 31, 32, 34, 38, 40, 42, 44, 46, 47, 51, 52, 53, 58, 59, 60, 64, 65, 66, 68, 69, 70, 72, 73, 75, 76, 77, 81, 83, 85, 90, 92, 93, 96, 98] Finding Minimum Weight Matching (greedy)... Minimum weight matching (MST + Matching Edges): [(14, 16, 13.04), (54, 82, 23.35), (51, 77, 40.00), (5, 48, 47.07), (27, 92, 47.68), (6, 8, 50.25), (24, 80, 51.48), (15, 87, 58.14), (20, 73, 66.01), (77, 95, 68.15), (4, 36, 69.23), (15, 21, 70.18), (39, 63, 80.45), (26, 85, 81.15), (40, 70, 82.68), (30, 79, 87.21), (58, 73, 89.00), (32, 75, 92.20), (6, 56, 93.05), (23, 37, 96.05), (90, 97, 99.41), (12, 32, 100.66), (39, 53, 103.59), (22, 97, 107.94), (10, 14, 113.07), (64, 65, 114.77), (70, 99, 121.43), (21, 93, 133.65), (11, 26, 136.49), (2, 45, 139.26), (9, 83, 140.04), (34, 85, 141.06), (33, 82, 145.96), (50, 86, 146.37), (3, 64, 147.18), (28, 45, 148.41), (59, 61, 151.05), (20, 71, 151.71), (67, 84, 153.40), (11, 19, 153.58), (5, 62, 157.26), (43, 49, 161.76), (24, 60, 164.71), (2, 42, 171.03), (44, 90, 179.04), (49, 72, 179.82), (41, 88, 180.42), (17, 23, 181.34), (21, 69, 186.13), (50, 60, 186.27), (30, 88, 187.77), (59, 76, 191.02), (8, 86, 195.49), (10, 31, 201.85), (19, 56, 203.75), (7, 91, 206.31), (48, 89, 208.12), (31, 90, 215.24), (12, 94, 215.64), (55, 96, 216.11), (2, 13, 221.30), (52, 78, 223.01), (35, 37, 227.76), (17, 78, 229.65), (52, 87, 229.75), (11, 54, 229.92), (18, 74, 234.08), (74, 96, 236.54), (4, 77, 242.17), (81, 94, 245.31), (7, 41, 247.78), (47, 99, 251.75), (1, 43, 256.56), (29, 38, 261.24), (67, 72, 266.72), (13, 70, 267.55), (9, 71, 269.65), (18, 52, 279.87), (55, 79, 280.80), (25, 64, 283.34), (38, 95, 283.64), (0, 62, 288.09), (68, 72, 292.88), (63, 68, 297.38), (35, 98, 299.71), (9, 89, 300.03), (28, 33, 300.50), (27, 66, 300.85), (0, 91, 302.55), (27, 57, 302.60), (68, 80, 303.12), (61, 85, 307.65), (3, 96, 324.23), (57, 60, 325.04), (10, 73, 328.69), (22, 59, 345.32), (31, 46, 361.33), (38, 84, 371.55), (32, 36, 407.77), (81, 75, 478.76), (72, 68, 292.88), (60, 66, 597.01), (31, 10, 201.85), (65, 64, 114.77), (34, 85, 141.06), (77, 51, 40.00), (40, 70, 82.68), (59, 76, 191.02), (47, 29, 427.13), (73, 58, 89.00), (53, 1, 331.30), (42, 2, 171.03), (11, 90, 1104.00), (98, 83, 593.33), (93, 21, 133.65), (96, 52, 589.48), (92, 27, 47.68), (9, 46, 656.41), (44, 16, 404.01), (25, 69, 669.16), (38, 32, 764.20)] Finding Eulerian Tour... Eulerian tour: [14, 10, 31, 46, 9, 89, 48, 5, 62, 0, 91, 7, 41, 88, 30, 79, 55, 96, 52, 18, 74, 96, 3, 64, 65, 64, 25, 69, 21, 93, 21, 15, 87, 52, 78, 17, 23, 37, 35, 98, 83, 9, 71, 20, 73, 58, 73, 10, 31, 90, 11, 54, 82, 33, 28, 45, 2, 42, 2, 13, 70, 40, 70, 99, 47, 29, 38, 32, 12, 94, 81, 75, 32, 36, 4, 77, 51, 77, 95, 38, 84, 67, 72, 68, 63, 39, 53, 1, 43, 49, 72, 68, 80, 24, 60, 66, 27, 92, 27, 57, 60, 50, 86, 8, 6, 56, 19, 11, 26, 85, 34, 85, 61, 59, 76, 59, 22, 97, 90, 44, 16, 14] Creating Hamiltonian path (shortcutting)... Result path: [14, 10, 31, 46, 9, 89, 48, 5, 62, 0, 91, 7, 41, 88, 30, 79, 55, 96, 52, 18, 74, 3, 64, 65, 25, 69, 21, 93, 15, 87, 78, 17, 23, 37, 35, 98, 83, 71, 20, 73, 58, 90, 11, 54, 82, 33, 28, 45, 2, 42, 13, 70, 40, 99, 47, 29, 38, 32, 12, 94, 81, 75, 36, 4, 77, 51, 95, 84, 67, 72, 68, 63, 39, 53, 1, 43, 49, 80, 24, 60, 66, 27, 92, 57, 50, 86, 8, 6, 56, 19, 26, 85, 34, 61, 59, 76, 22, 97, 44, 16, 14] Result length of the path: 25535.83
Wren
import "./dynamic" for Struct
import "./set" for Set
import "random" for Random
/* Helper structs. */
var Point = Struct.create("Point", ["x", "y", "id"])
var Edge = Struct.create("Edge", ["u", "v", "weight", "id"])
var Neighbor = Struct.create("Neighbor", ["neighbor", "edgeId"])
/* Random number generator. */
var rand = Random.new()
/* Helper functions. */
var toFixed2 = Fn.new { |f| (f * 100).round / 100 }
var printContainer = Fn.new { |container, name|
System.print("%(name): [%(container.join(", "))]")
}
var printEdges = Fn.new { |edges, name|
var edgeStrings = edges.map { |e| "%(e.u), %(e.v), %(toFixed2.call(e.weight))" }
System.print("%(name): [%(edgeStrings.join(", "))]")
}
var printGraph = Fn.new { |graph, name|
System.print("%(name): {")
var n = graph.count
for (i in 0...n) {
var entries = []
for (j in 0...n) {
if (i != j) entries.add("%(j): %(toFixed2.call(graph[i][j]))")
}
System.print(" %(i): {%(entries.join(", "))}%(i == n - 1 ? "" : ",")")
}
System.print("}")
}
/* Euclidean distance. */
var getLength = Fn.new { |p1, p2|
var dx = p1.x - p2.x
var dy = p1.y - p2.y
return (dx * dx + dy * dy).sqrt
}
/* Build Complete Graph (Adjacency Matrix). */
var buildGraph = Fn.new { |data|
var n = data.count
// Initialize n x n list with 0's.
var graph = List.filled(n, null)
for (i in 0...n) graph[i] = List.filled(n, 0)
for (i in 0...n) {
// Only calculate upper triangle.
var j = i + 1
while (j < n) {
var dist = getLength.call(data[i], data[j])
graph[i][j] = dist
graph[j][i] = dist // symmetric graph
j = j + 1
}
}
return graph
}
/* Union-Find Data Structure. */
class UnionFind {
construct new(n) {
// Initialize parent list: each node is its own parent.
_parent = List.filled(n, 0)
for (i in 0...n) _parent[i] = i
// Initialize rank (or size) list for optimization.
_rank = List.filled(n, 0)
}
find(i) {
if (_parent[i] == i) return i
// Path compression: point directly to the root.
_parent[i] = find(_parent[i])
return _parent[i]
}
unite(i, j) {
var rootI = find(i)
var rootJ = find(j)
if (rootI != rootJ) {
// Union by rank: attach smaller rank tree under larger rank tree.
if (_rank[rootI] < _rank[rootJ]) {
_parent[rootI] = rootJ
} else if (_rank[rootI] > _rank[rootJ]) {
_parent[rootJ] = rootI
} else {
// Ranks are equal, choose one as parent and increment its rank.
_parent[rootJ] = rootI
_rank[rootI] = _rank[rootI] + 1
}
return true // successfully united
}
return false // already in the same set
}
}
/* Minimum Spanning Tree (Kruskal's Algorithm). */
var minimumSpanningTree = Fn.new { |graph|
var n = graph.count
if (n == 0) return []
var edges = []
for (i in 0...n) {
// Avoid duplicates and self-loops.
var j = i + 1
while (j < n) {
edges.add(Edge.new(i, j, graph[i][j], null))
j = j + 1
}
}
// Sort edges by weight in ascending order.
edges.sort { |a, b| a.weight < b.weight }
var mst = []
var uf = UnionFind.new(n)
var edgesCount = 0
for (edge in edges) {
// If uniting forms a new connection:
if (uf.unite(edge.u, edge.v)) {
mst.add(edge)
edgesCount = edgesCount + 1
if (edgesCount == n - 1) break // optimization: MST has n-1 edges
}
}
return mst
}
/* Find Vertices with Odd Degree in MST. */
var findOddVertices = Fn.new { |mst, n|
var degree = List.filled(n, 0)
for (edge in mst) {
degree[edge.u] = degree[edge.u] + 1
degree[edge.v] = degree[edge.v] + 1
}
var oddVertices = []
for (i in 0...n) {
if (degree[i] % 2 != 0) oddVertices.add(i)
}
return oddVertices
}
/* Minimum Weight Matching (Greedy Heuristic)
Fisher-Yates (Knuth) Shuffle algorithm. */
var shuffleArray = Fn.new { |array|
rand.shuffle(array) // uses F-Y uder the hood
}
// Note: This modifies the mst array by adding matching edges.
var minimumWeightMatching = Fn.new { |mst, graph, oddVertices|
// Shuffle for randomness.
shuffleArray.call(oddVertices)
// Keep track of vertices already matched in this phase.
var matched = Set.new()
for (i in 0...oddVertices.count) {
var v = oddVertices[i]
if (matched.contains(v)) continue // skip if already matched
var minLength = Num.infinity
var closestU = -1
// Find the closest *unmatched* odd vertex.
var j = i + 1
while (j < oddVertices.count) {
var u = oddVertices[j]
if (!matched.contains(u)) { // check if 'u' is available
if (graph[v][u] < minLength) {
minLength = graph[v][u]
closestU = u
}
}
j = j + 1
}
if (closestU != -1) {
// Add the matching edge to the MST list (now a multigraph).
mst.add(Edge.new(v, closestU, minLength, null))
matched.add(v)
matched.add(closestU) // mark both as matched
}
// Christofides guarantees an even number of odd-degree vertices,
// so every vertex *should* find a match in a perfect matching scenario.
// The greedy approach might leave some unmatched if not careful,
// but this loop structure should work.
}
}
/* Find Eulerian Tour (Hierholzer's Algorithm). */
var findEulerianTour = Fn.new { |matchedMST, n|
if (matchedMST.isEmpty) return []
// Assign unique IDs to edges for tracking in the multigraph
// (essential because multiple edges can exist between nodes).
for (i in 0...matchedMST.count) matchedMST[i].id = i
// Build adjacency list.
var adj = List.filled(n, null)
for (i in 0...n) adj[i] = []
var edgeUsed = Set.new() // store used edge IDs
for (edge in matchedMST) {
adj[edge.u].add(Neighbor.new(edge.v, edge.id))
adj[edge.v].add(Neighbor.new(edge.u, edge.id))
}
var tour = []
var stack = []
// Start at any vertex with edges (guaranteed to exist if matchedMST is not empty).
var startNode = matchedMST[0].u
stack.add(startNode)
var currentNode = startNode
while (!stack.isEmpty) {
currentNode = stack[-1] // peek top of stack
var foundUnusedEdge = false
// Check outgoing edges from currentNode.
while (!adj[currentNode].isEmpty) {
// Look at last edge (efficient removal).
var edgeInfo = adj[currentNode][adj[currentNode].count - 1]
if (!edgeUsed.contains(edgeInfo.edgeId)) {
// Found an unused edge.
edgeUsed.add(edgeInfo.edgeId) // mark edge as used
stack.add(edgeInfo.neighbor) // push neighbor onto stack
adj[currentNode].removeAt(-1) // remove edge for current node's list
currentNode = edgeInfo.neighbor // move to the neighbor
foundUnusedEdge = true
break // break inner loop to process the new currentNode
} else {
// This edge was already used (by traversal from the other side), remove it.
adj[currentNode].removeAt(-1)
}
}
if (!foundUnusedEdge) {
// If no unused edges from currentNode, backtrack.
tour.add(stack.removeAt(-1))
}
}
// The tour is constructed in reverse order by Hierholzer's algorithm.
return tour[-1..0]
}
/* Main TSP Function (Christofides Approximation). */
var tsp = Fn.new { |data|
var n = data.count
if (n == 0) return [0, []]
if (n == 1) {
// If data points have explicit IDs use them, otherwise use index 0.
var startId = data[0].id ? data[0].id : 0
return [0, [startId]]
}
// Assign IDs if they don't exist, assuming order corresponds to 0..n-1.
for (i in 0...data.count) if (!data[i].id) data[i].id = i
System.print("Bulding graph...")
var G = buildGraph.call(data)
// printGraph.call(G, "Graph") // often too large to print meaningfully
System.print("\nFinding Minimum Spanning Tree...")
var MSTree = minimumSpanningTree.call(G)
printEdges.call(MSTree, "MSTree")
System.print("\nFinding odd degree vertices...")
var oddVertices = findOddVertices.call(MSTree, n)
printContainer.call(oddVertices, "Odd vertices in MSTree")
System.print("\nFinding Minimum Weight Matching (greedy)...")
// Create a new list containing MST edges to avoid modifying the original
// MSTree variable directly. The matching edges will be added to this new list.
var multigraphEdges = MSTree.toList
// Pass a copy of oddVertices as it might be shuffled in place.
minimumWeightMatching.call(multigraphEdges, G, oddVertices.toList)
printEdges.call(multigraphEdges, "Minimum weight matching (MST + Matching Edges)")
System.print("\nFinding Eulerian Tour...")
var eulerianTour = findEulerianTour.call(multigraphEdges, n)
printContainer.call(eulerianTour, "Eulerian tour")
/* Create Hamiltonian Circuit by Skipping Visited Nodes. */
if (eulerianTour.isEmpty && n > 0) {
System.print("Error: Eulerian tour could not be found.")
return [-1, []] // indicate error
}
System.print("\nCreating Hamiltonian path (shortcutting)...")
var path = []
var length = 0
var visited = Set.new() // use Set for efficient O(1) average time complexity checks
var currentPathNode = -1 // track the last node added to the *Hamiltonian* path
for (node in eulerianTour) {
if (!visited.contains(node)) {
path.add(node)
visited.add(node)
if (currentPathNode != -1) { // add edge length from the previous node *in the path*
length = length + G[currentPathNode][node]
}
currentPathNode = node // update the last node added to the path
}
}
// Add the edge back to the start to complete the cycle.
if (!path.isEmpty) {
// Edge from last node in path to first node.
length = length + G[currentPathNode][path[0]]
// Add the starting node again to show the cycle.
path.add(path[0])
}
printContainer.call(path, "Result path")
System.print("Result length of the path: %(toFixed2.call(length))")
}
// Input Data
var rawData = [
[1380, 939], [2848, 96], [3510, 1671], [457, 334], [3888, 666], [984, 965], [2721, 1482], [1286, 525],
[2716, 1432], [738, 1325], [1251, 1832], [2728, 1698], [3815, 169], [3683, 1533], [1247, 1945], [123, 862],
[1234, 1946], [252, 1240], [611, 673], [2576, 1676], [928, 1700], [53, 857], [1807, 1711], [274, 1420],
[2574, 946], [178, 24], [2678, 1825], [1795, 962], [3384, 1498], [3520, 1079], [1256, 61], [1424, 1728],
[3913, 192], [3085, 1528], [2573, 1969], [463, 1670], [3875, 598], [298, 1513], [3479, 821], [2542, 236],
[3955, 1743], [1323, 280], [3447, 1830], [2936, 337], [1621, 1830], [3373, 1646], [1393, 1368],
[3874, 1318], [938, 955], [3022, 474], [2482, 1183], [3854, 923], [376, 825], [2519, 135], [2945, 1622],
[953, 268], [2628, 1479], [2097, 981], [890, 1846], [2139, 1806], [2421, 1007], [2290, 1810], [1115, 1052],
[2588, 302], [327, 265], [241, 341], [1917, 687], [2991, 792], [2573, 599], [19, 674], [3911, 1673],
[872, 1559], [2863, 558], [929, 1766], [839, 620], [3893, 102], [2178, 1619], [3822, 899], [378, 1048],
[1178, 100], [2599, 901], [3416, 143], [2961, 1605], [611, 1384], [3113, 885], [2597, 1830], [2586, 1286],
[161, 906], [1429, 134], [742, 1025], [1625, 1651], [1187, 706], [1787, 1009], [22, 987], [3640, 43],
[3756, 882], [776, 392], [1724, 1642], [198, 1810], [3950, 1558]
]
// Convert raw data to Point objects, using index as ID.
var dataPoints = (0...rawData.count).map { |i|
var x = rawData[i][0]
var y = rawData[i][1]
var id = i // Use 0-based index as the vertex ID
return Point.new(x, y, id)
}.toList
// Runs TSP
tsp.call(dataPoints)
- Output:
Sample run.
Bulding graph... Finding Minimum Spanning Tree... MSTree: [14, 16, 13.04, 54, 82, 23.35, 51, 77, 40, 5, 48, 47.07, 27, 92, 47.68, 6, 8, 50.25, 24, 80, 51.48, 15, 87, 58.14, 20, 73, 66.01, 77, 95, 68.15, 4, 36, 69.23, 15, 21, 70.18, 39, 63, 80.45, 26, 85, 81.15, 40, 70, 82.68, 30, 79, 87.21, 58, 73, 89, 32, 75, 92.2, 6, 56, 93.05, 23, 37, 96.05, 90, 97, 99.41, 12, 32, 100.66, 39, 53, 103.59, 22, 97, 107.94, 10, 14, 113.07, 64, 65, 114.77, 70, 99, 121.43, 21, 93, 133.65, 11, 26, 136.49, 2, 45, 139.26, 9, 83, 140.04, 34, 85, 141.06, 33, 82, 145.96, 50, 86, 146.37, 3, 64, 147.18, 28, 45, 148.41, 59, 61, 151.05, 20, 71, 151.71, 67, 84, 153.4, 11, 19, 153.58, 5, 62, 157.26, 43, 49, 161.76, 24, 60, 164.71, 2, 42, 171.03, 44, 90, 179.04, 49, 72, 179.82, 41, 88, 180.42, 17, 23, 181.34, 21, 69, 186.13, 50, 60, 186.27, 30, 88, 187.77, 59, 76, 191.02, 8, 86, 195.49, 10, 31, 201.85, 19, 56, 203.75, 7, 91, 206.31, 48, 89, 208.12, 31, 90, 215.24, 12, 94, 215.64, 55, 96, 216.11, 2, 13, 221.3, 52, 78, 223.01, 35, 37, 227.76, 17, 78, 229.65, 52, 87, 229.75, 11, 54, 229.92, 18, 74, 234.08, 74, 96, 236.54, 4, 77, 242.17, 81, 94, 245.31, 7, 41, 247.78, 47, 99, 251.75, 1, 43, 256.56, 29, 38, 261.24, 67, 72, 266.72, 13, 70, 267.55, 9, 71, 269.65, 18, 52, 279.87, 55, 79, 280.8, 25, 64, 283.34, 38, 95, 283.64, 0, 62, 288.09, 68, 72, 292.88, 63, 68, 297.38, 35, 98, 299.71, 9, 89, 300.03, 28, 33, 300.5, 27, 66, 300.85, 0, 91, 302.55, 27, 57, 302.6, 68, 80, 303.12, 61, 85, 307.65, 3, 96, 324.23, 57, 60, 325.04, 10, 73, 328.69, 22, 59, 345.32, 31, 46, 361.33, 38, 84, 371.55, 32, 36, 407.77] Finding odd degree vertices... Odd vertices in MSTree: [1, 2, 9, 10, 11, 16, 21, 25, 27, 29, 31, 32, 34, 38, 40, 42, 44, 46, 47, 51, 52, 53, 58, 59, 60, 64, 65, 66, 68, 69, 70, 72, 73, 75, 76, 77, 81, 83, 85, 90, 92, 93, 96, 98] Finding Minimum Weight Matching (greedy)... Minimum weight matching (MST + Matching Edges): [14, 16, 13.04, 54, 82, 23.35, 51, 77, 40, 5, 48, 47.07, 27, 92, 47.68, 6, 8, 50.25, 24, 80, 51.48, 15, 87, 58.14, 20, 73, 66.01, 77, 95, 68.15, 4, 36, 69.23, 15, 21, 70.18, 39, 63, 80.45, 26, 85, 81.15, 40, 70, 82.68, 30, 79, 87.21, 58, 73, 89, 32, 75, 92.2, 6, 56, 93.05, 23, 37, 96.05, 90, 97, 99.41, 12, 32, 100.66, 39, 53, 103.59, 22, 97, 107.94, 10, 14, 113.07, 64, 65, 114.77, 70, 99, 121.43, 21, 93, 133.65, 11, 26, 136.49, 2, 45, 139.26, 9, 83, 140.04, 34, 85, 141.06, 33, 82, 145.96, 50, 86, 146.37, 3, 64, 147.18, 28, 45, 148.41, 59, 61, 151.05, 20, 71, 151.71, 67, 84, 153.4, 11, 19, 153.58, 5, 62, 157.26, 43, 49, 161.76, 24, 60, 164.71, 2, 42, 171.03, 44, 90, 179.04, 49, 72, 179.82, 41, 88, 180.42, 17, 23, 181.34, 21, 69, 186.13, 50, 60, 186.27, 30, 88, 187.77, 59, 76, 191.02, 8, 86, 195.49, 10, 31, 201.85, 19, 56, 203.75, 7, 91, 206.31, 48, 89, 208.12, 31, 90, 215.24, 12, 94, 215.64, 55, 96, 216.11, 2, 13, 221.3, 52, 78, 223.01, 35, 37, 227.76, 17, 78, 229.65, 52, 87, 229.75, 11, 54, 229.92, 18, 74, 234.08, 74, 96, 236.54, 4, 77, 242.17, 81, 94, 245.31, 7, 41, 247.78, 47, 99, 251.75, 1, 43, 256.56, 29, 38, 261.24, 67, 72, 266.72, 13, 70, 267.55, 9, 71, 269.65, 18, 52, 279.87, 55, 79, 280.8, 25, 64, 283.34, 38, 95, 283.64, 0, 62, 288.09, 68, 72, 292.88, 63, 68, 297.38, 35, 98, 299.71, 9, 89, 300.03, 28, 33, 300.5, 27, 66, 300.85, 0, 91, 302.55, 27, 57, 302.6, 68, 80, 303.12, 61, 85, 307.65, 3, 96, 324.23, 57, 60, 325.04, 10, 73, 328.69, 22, 59, 345.32, 31, 46, 361.33, 38, 84, 371.55, 32, 36, 407.77, 98, 83, 593.33, 58, 73, 89, 52, 21, 324.58, 42, 2, 171.03, 81, 75, 478.76, 93, 69, 313.01, 25, 64, 283.34, 70, 40, 82.68, 34, 85, 141.06, 65, 96, 537.43, 16, 10, 115.26, 68, 72, 292.88, 59, 76, 191.02, 92, 27, 47.68, 29, 38, 261.24, 9, 46, 656.41, 53, 1, 331.3, 51, 77, 40, 47, 32, 1126.68, 90, 44, 179.04, 11, 60, 756.13, 31, 66, 1151.84] Finding Eulerian Tour... Eulerian tour: [14, 10, 73, 58, 73, 20, 71, 9, 89, 48, 5, 62, 0, 91, 7, 41, 88, 30, 79, 55, 96, 65, 64, 25, 64, 3, 96, 74, 18, 52, 21, 69, 93, 21, 15, 87, 52, 78, 17, 23, 37, 35, 98, 83, 9, 46, 31, 66, 27, 92, 27, 57, 60, 11, 54, 82, 33, 28, 45, 2, 42, 2, 13, 70, 40, 70, 99, 47, 32, 12, 94, 81, 75, 32, 36, 4, 77, 51, 77, 95, 38, 29, 38, 84, 67, 72, 68, 63, 39, 53, 1, 43, 49, 72, 68, 80, 24, 60, 50, 86, 8, 6, 56, 19, 11, 26, 85, 34, 85, 61, 59, 76, 59, 22, 97, 90, 44, 90, 31, 10, 16, 14] Creating Hamiltonian path (shortcutting)... Result path: [14, 10, 73, 58, 20, 71, 9, 89, 48, 5, 62, 0, 91, 7, 41, 88, 30, 79, 55, 96, 65, 64, 25, 3, 74, 18, 52, 21, 69, 93, 15, 87, 78, 17, 23, 37, 35, 98, 83, 46, 31, 66, 27, 92, 57, 60, 11, 54, 82, 33, 28, 45, 2, 42, 13, 70, 40, 99, 47, 32, 12, 94, 81, 75, 36, 4, 77, 51, 95, 38, 29, 84, 67, 72, 68, 63, 39, 53, 1, 43, 49, 80, 24, 50, 86, 8, 6, 56, 19, 26, 85, 34, 61, 59, 76, 22, 97, 90, 44, 16, 14] Result length of the path: 25358.71