All Downloads are FREE. Search and download functionalities are using the official Maven repository.

smile.graph.Graph Maven / Gradle / Ivy

There is a newer version: 4.2.0
Show newest version
/*
 * Copyright (c) 2010-2021 Haifeng Li. All rights reserved.
 *
 * Smile is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Smile is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Smile.  If not, see .
 */

package smile.graph;

import java.util.*;
import java.util.stream.DoubleStream;
import smile.math.MathEx;
import smile.math.matrix.IMatrix;
import smile.sort.Sort;
import smile.util.ArrayElementConsumer;
import smile.util.ArrayElementFunction;
import smile.util.PriorityQueue;
import smile.util.Strings;

/**
 * A graph is an abstract representation of a set of objects where some pairs
 * of the objects are connected by links. The interconnected objects are
 * represented by mathematical abstractions called vertices, and the links
 * that connect some pairs of vertices are called edges. The edges may be
 * directed (asymmetric) or undirected (symmetric). A graph is a weighted graph
 * if a number (weight) is assigned to each edge. Such weights might represent,
 * for example, costs, lengths or capacities, etc., depending on the problem.
 *
 * @author Haifeng Li
 */
public abstract class Graph {
    private static final org.slf4j.Logger logger = org.slf4j.LoggerFactory.getLogger(Graph.class);
    /**
     * Is the graph directed?
     */
    private final boolean digraph;

    /**
     * Graph edge.
     * @param u the vertex id. For directed graph,
     *          this is the tail of arc.
     * @param v the other vertex id. For directed graph,
     *          this is the head of arc.
     * @param weight the weight of edge. For unweighted graph,
     *               this is always 1.
     */
    public record Edge(int u, int v, double weight) implements Comparable {
        /**
         * Constructor of unweighted edge.
         * @param u the vertex id.
         * @param v the other vertex id.
         */
        public Edge(int u, int v) {
            this(u, v, 1.0);
        }

        @Override
        public int compareTo(Edge o) {
            return Double.compare(weight, o.weight);
        }
    }

    /**
     * Constructor.
     * @param digraph true if this is a directed graph.
     */
    public Graph(boolean digraph) {
        this.digraph = digraph;
    }

    /**
     * Return true if the graph is directed.
     * @return true if the graph is directed.
     */
    public boolean isDigraph() {
        return digraph;
    }

    /**
     * Returns the graphic representation in Graphviz dot format.
     * Try http://viz-js.com/
     * to visualize the returned string.
     * @return the graphic representation in Graphviz dot format.
     */
    public String dot() {
        return dot(null, null);
    }

    /**
     * Returns the graphic representation in Graphviz dot format.
     * Try http://viz-js.com/
     * to visualize the returned string.
     * @param name the graph name.
     * @param label the label of nodes.
     * @return the graphic representation in Graphviz dot format.
     */
    public String dot(String name, String[] label) {
        StringBuilder builder = new StringBuilder();
        builder.append(digraph ? "digraph " : "graph ");
        if (name != null) builder.append(name);
        builder.append(" {\n");
        builder.append("  node [shape=box, style=\"rounded\", color=\"black\", fontname=helvetica];\n");
        builder.append("  edge [fontname=helvetica];\n");

        if (label != null) {
            for (int i = 0; i < label.length; i++) {
                builder.append(String.format("  %d [label=\"%s\"];\n", i, label[i]));
            }
        }

        int n = getVertexCount();
        String edge = digraph ? "->" : "--";
        for (int i = 0; i < n; i++) {
            int u = i;
            forEachEdge(i, (v, w) -> {
                if (digraph || v >= u) {
                    builder.append(String.format("  %d %s %d [label=\"%s\"];\n", u, edge, v, Strings.format(w)));
                }
            });
        }

        builder.append("}");
        return builder.toString();
    }

    /**
     * Returns the (dense or sparse) matrix representation of the graph.
     * @return the matrix representation of the graph.
     */
    public abstract IMatrix toMatrix();

    /**
     * Returns a subgraph containing all given vertices.
     * @param vertices the vertices to be included in subgraph.
     * @return a subgraph containing all given vertices
     */
    public abstract Graph subgraph(int[] vertices);
    
    /**
     * Returns the number of vertices.
     * @return the number of vertices.
     */
    public abstract int getVertexCount();

    /**
     * Returns true if and only if this graph contains an edge going
     * from the source vertex to the target vertex. In undirected graphs the
     * same result is obtained when source and target are inverted.
     *
     * @param source the id of source vertex of the edge.
     * @param target the id of target vertex of the edge.
     * @return true if this graph contains the specified edge.
     */
    public abstract boolean hasEdge(int source, int target);

    /**
     * Returns the weight assigned to a given edge. Unweighted graphs always
     * return 1.0. If the edge doesn't exist, it returns zero. For minimization
     * problems such as TSP, use getDistance as it will return infinity instead.
     *
     * @param source the id of source vertex of the edge.
     * @param target the id of target vertex of the edge.
     * @return the edge weight
     */
    public abstract double getWeight(int source, int target);

    /**
     * Returns the distance between two vertices.
     * @param source the id of source vertex of the edge.
     * @param target the id of target vertex of the edge.
     * @return the distance between two vertices.
     */
    public double getDistance(int source, int target) {
        double weight = getWeight(source, target);
        return weight == 0 ? Double.POSITIVE_INFINITY : weight;
    }

    /**
     * Sets the weight assigned to a given edge.
     *
     * @param source the id of source vertex of the edge.
     * @param target the id of target vertex of the edge.
     * @param weight the edge weight
     * @return this graph.
     */
    public abstract Graph setWeight(int source, int target, double weight);

    /**
     * Returns the edges from the specified vertex. If no edges are
     * touching the specified vertex returns an empty set.
     *
     * @param vertex the id of vertex for which a set of touching edges is to be
     * returned.
     * @return the edges touching the specified vertex.
     */
    public abstract List getEdges(int vertex);

    /**
     * Performs an action for each edge of a vertex.
     * @param vertex the vertex id.
     * @param action a non-interfering action to perform on the edges.
     */
    public abstract void forEachEdge(int vertex, ArrayElementConsumer action);

    /**
     * Returns a stream consisting of the results of applying the given
     * function to the edge weights of a vertex.
     * @param vertex the vertex id.
     * @param mapper a non-interfering, stateless function to map each
     *               edge weight of a vertex.
     * @return the stream of the new values of edge weights.
     */
    public abstract DoubleStream mapEdges(int vertex, ArrayElementFunction mapper);

    /**
     * Updates the edge weights of a vertex.
     * @param vertex the vertex id.
     * @param mapper a function to map each edge weight to new value.
     */
    public abstract void updateEdges(int vertex, ArrayElementFunction mapper);

    /**
     * Creates a new edge in this graph, going from the source vertex to the
     * target vertex, and returns the created edge.
     *
     * @param source the id of source vertex of the edge.
     * @param target the id of target vertex of the edge.
     */
    public Graph addEdge(int source, int target) {
        return addEdge(source, target, 1.0);
    }

    /**
     * Creates a new edge in this graph, going from the source vertex to the
     * target vertex, and returns the created edge.
     *
     * @param source the id of source vertex of the edge.
     * @param target the id of target vertex of the edge.
     * @param weight the weight of edge.
     */
    public Graph addEdge(int source, int target, double weight) {
        return setWeight(source, target, weight);
    }

    /**
     * Adds a set of edges to the graph.
     *
     * @param edges edges to be added to this graph.
     */
    public Graph addEdges(Collection edges) {
        for (Edge edge : edges) {
            setWeight(edge.u, edge.v, edge.weight);
        }
        return this;
    }

    /**
     * Removes a set of edges from the graph.
     *
     * @param edges edges to be removed from this graph.
     */
    public Graph removeEdges(Collection edges) {
        for (Edge edge : edges) {
            removeEdge(edge.u, edge.v);
        }
        return this;
    }

    /**
     * In a simple graph, removes and returns the edge going from the specified source
     * vertex to the specified target vertex.
     *
     * @param source the id of source vertex of the edge.
     * @param target the id of target vertex of the edge.
     */
    public Graph removeEdge(int source, int target) {
        return setWeight(source, target, 0.0);
    }

    /**
     * Returns the degree of the specified vertex. A degree of a vertex in an
     * undirected graph is the number of edges touching that vertex.
     *
     * @param vertex the id of vertex.
     * @return the degree of the specified vertex.
     */
    public int getDegree(int vertex) {
        return digraph ? getInDegree(vertex) + getOutDegree(vertex) : getOutDegree(vertex);
    }

    /**
     * Returns the in-degree of the specified vertex. An in-degree of a vertex in a
     * directed graph is the number of edges head to that vertex.
     *
     * @param vertex the id of vertex.
     * @return the degree of the specified vertex.
     */
    public abstract int getInDegree(int vertex);

    /**
     * Returns the out-degree of the specified vertex. An out-degree of a vertex in a
     * directed graph is the number of edges from that vertex.
     *
     * @param vertex the id of vertex.
     * @return the degree of the specified vertex.
     */
    public abstract int getOutDegree(int vertex);

    /**
     * Reverse topological sort digraph by depth-first search of graph.
     * @param v the start vertex.
     * @param visited the flag if vertex has been visited.
     * @param order the array to store the reverse topological order.
     * @param count the number of vertices have been visited before this search.
     *              It will be updated after this search.
     */
    private void dfsort(int v, boolean[] visited, int[] order, int[] count) {
        visited[v] = true;

        forEachEdge(v, (u, w) -> {
            if (!visited[u]) {
                dfsort(u, visited, order, count);
            }
        });

        order[count[0]++] = v;
    }

    /**
     * Reverse topological sort digraph by depth-first search of graph.
     *
     * @return the vertices in the reverse topological order.
     */
    public int[] dfsort() {
        if (!digraph) {
            throw new UnsupportedOperationException("Topological sort cannot be applied on undirected graph.");
        }

        int n = getVertexCount();
        boolean[] visited = new boolean[n];
        int[] order = new int[n];
        Arrays.fill(order, -1);

        int[] count = new int[1];
        for (int i = 0; i < n; i++) {
            if (!visited[i]) {
                dfsort(i, visited, order, count);
            }
        }

        return order;
    }

    /**
     * Depth-first search connected components of graph.
     * @param v the start vertex.
     * @param cc the array to store the connected component id of each vertex.
     * @param id the current component id.
     */
    private void dfcc(int v, int[] cc, int id) {
        cc[v] = id;
        forEachEdge(v, (u, w) -> {
            if (cc[u] == -1) {
                dfcc(u, cc, id);
            }
        });
    }

    /**
     * Returns the connected components by depth-first search.
     *
     * @return a two-dimensional array of which each row is the vertices
     *         in the same connected component.
     */
    public int[][] dfcc() {
        if (digraph) {
            throw new UnsupportedOperationException("Connected components algorithm cannot be applied on digraph");
        }

        int n = getVertexCount();
        int[] cc = new int[n];
        Arrays.fill(cc, -1);

        int id = 0;
        for (int i = 0; i < n; i++) {
            if (cc[i] == -1) {
                dfcc(i, cc, id++);
            }
        }

        return connectedComponents(id, cc);
    }

    /**
     * Return the connected components.
     * @param numComponents the number of connected components.
     * @param cc the component id of each vertex.
     * @return the connected components.
     */
    private int[][] connectedComponents(int numComponents, int[] cc) {
        int n = cc.length;
        int[] size = new int[numComponents];
        for (int c : cc) {
            size[c]++;
        }

        int[][] components = new int[numComponents][];
        for (int i = 0; i < numComponents; i++) {
            components[i] = new int[size[i]];
            for (int j = 0, k = 0; j < n; j++) {
                if (cc[j] == i) {
                    components[i][k++] = j;
                }
            }
            Arrays.sort(components[i]);
        }

        return components;
    }

    /**
     * Depth-first search of graph.
     * @param v the start vertex.
     * @param visited the flag if vertex has been visited.
     */
    private void dfs(VertexVisitor visitor, int v, boolean[] visited) {
        visitor.accept(v);
        visited[v] = true;
        forEachEdge(v, (u, w) -> {
            if (!visited[u]) dfs(visitor, u, visited);
        });
    }

    /**
     * DFS search on graph and performs some operation defined in visitor
     * on each vertex during traveling.
     * @param visitor the visitor functor.
     */
    public void dfs(VertexVisitor visitor) {
        int n = getVertexCount();
        boolean[] visited = new boolean[n];

        for (int i = 0; i < n; i++) {
            if (!visited[i]) {
                dfs(visitor, i, visited);
            }
        }
    }

    /**
     * Topological sort digraph by breadth-first search of graph.
     *
     * @return the vertices in the topological order.
     */
    public int[] bfsort() {
        if (!digraph) {
            throw new UnsupportedOperationException("Topological sort cannot be applied on undirected graph.");
        }

        int n = getVertexCount();
        int[] in = new int[n];
        int[] ts = new int[n];
        for (int i = 0; i < n; i++) {
            ts[i] = -1;
            forEachEdge(i, (j, w) -> in[j]++);
        }

        Queue queue = new LinkedList<>();
        for (int i = 0; i < n; i++) {
            if (in[i] == 0) {
                queue.offer(i);
            }
        }

        for (int i = 0; !queue.isEmpty(); i++) {
            int u = queue.poll();
            ts[i] = u;
            forEachEdge(u, (v, w) -> {
                if (--in[v] == 0) queue.offer(v);
            });
        }

        return ts;
    }

    /**
     * Breadth-first search connected components of graph.
     * @param v the start vertex.
     * @param cc the array to store the connected component id of each vertex.
     * @param id the current component id.
     */
    private void bfcc(int v, int[] cc, int id) {
        cc[v] = id;
        Queue queue = new LinkedList<>();
        queue.offer(v);
        while (!queue.isEmpty()) {
            int u = queue.poll();
            forEachEdge(u, (i, w) -> {
                if (cc[i] == -1) {
                    queue.offer(i);
                    cc[i] = id;
                }
            });
        }
    }

    /**
     * Returns the connected components by breadth-first search.
     *
     * @return a two-dimensional array of which each row is the vertices
     *         in the same connected component.
     */
    public int[][] bfcc() {
        if (digraph) {
            throw new UnsupportedOperationException("Connected components algorithm cannot be applied on digraph");
        }

        int n = getVertexCount();
        int[] cc = new int[n];
        Arrays.fill(cc, -1);

        int id = 0;
        for (int i = 0; i < n; i++) {
            if (cc[i] == -1) {
                bfcc(i, cc, id++);
            }
        }

        return connectedComponents(id, cc);
    }

    /**
     * Breadth-first search of graph.
     * @param visitor the visitor functor.
     * @param v the start vertex.
     * @param visited the flag if vertex has been visited.
     * @param queue a queue of vertices to visit.
     */
    private void bfs(VertexVisitor visitor, int v, boolean[] visited, Queue queue) {
        visitor.accept(v);
        visited[v] = true;
        queue.offer(v);
        while (!queue.isEmpty()) {
            int u = queue.poll();
            forEachEdge(u, (i, w) -> {
                if (!visited[i]) {
                    visitor.accept(i);
                    queue.offer(i);
                    visited[i] = true;
                }
            });
        }
    }

    /**
     * BFS search on graph and performs some operation defined in visitor
     * on each vertex during traveling.
     * @param visitor the visitor functor.
     */
    public void bfs(VertexVisitor visitor) {
        int n = getVertexCount();
        boolean[] visited = new boolean[n];

        Queue queue = new LinkedList<>();
        for (int i = 0; i < n; i++) {
            if (!visited[i]) {
                bfs(visitor, i, visited, queue);
            }
        }
    }

    /**
     * Calculate the shortest path from a source to all other vertices in the
     * graph by Dijkstra algorithm.
     *
     * @param s the source vertex.
     * @return The distance to all vertices from the source.
     */
    public double[] dijkstra(int s) {
        return dijkstra(s, true);
    }
    
    /**
     * Calculate the shortest path from a source to all other vertices in the
     * graph by Dijkstra algorithm.
     * @param s The source vertex.
     * @param weighted Ignore edge weights if false.
     * @return The distance to all vertices from the source. If weighted is false,
     *         it is the length of the shortest path to other vertices.
     */
    public double[] dijkstra(int s, boolean weighted) {
        int n = getVertexCount();
        double[] wt = new double[n];
        Arrays.fill(wt, Double.POSITIVE_INFINITY);

        PriorityQueue queue = new PriorityQueue(wt);
        for (int v = 0; v < n; v++) {
            queue.insert(v);
        }

        wt[s] = 0.0;
        queue.lower(s);

        while (!queue.isEmpty()) {
            int v = queue.poll();
            if (!Double.isInfinite(wt[v])) {
                forEachEdge(v, (u, weight) -> {
                    double p = wt[v] + (weighted ? weight : 1);
                    if (p < wt[u]) {
                        wt[u] = p;
                        queue.lower(u);
                    }
                });
            }
        }
        
        return wt;
    }

    /**
     * Calculates the all pair shortest-path by Dijkstra algorithm.
     *
     * @return the length of shortest-path between vertices.
     */
    public double[][] dijkstra() {
        int n = getVertexCount();
        double[][] wt = new double[n][];
        for (int i = 0; i < n; i++) {
            wt[i] = dijkstra(i);
        }
        return wt;
    }

    /**
     * Returns the minimum spanning tree (MST) for a weighted undirected
     * graph by Prim's algorithm. MST is a subset of the edges that forms
     * a tree that includes every vertex, where the total weight of all
     * the edges in the tree is minimized.
     * @param mst an output container to hold edges in the MST.
     * @return the cost of minimum spanning tree.
     */
    public double prim(List mst) {
        if (digraph) {
            throw new UnsupportedOperationException("Prim's algorithm cannot be applied on a digraph.");
        }

        int n = getVertexCount();
        if (n < 2) {
            throw new UnsupportedOperationException("Cannot construct MST with fewer than 2 vertices.");
        }

        // Tracks whether a node is included in the MST
        boolean[] inMST = new boolean[n];

        // Stores the minimum edge weight to add a node to the MST
        double[] minEdgeWeight = new double[n];
        Arrays.fill(minEdgeWeight, Double.MAX_VALUE);

        // Stores the parent of each node in the MST
        int[] parent = new int[n];
        Arrays.fill(parent, -1);

        // Total weight of the MST
        double totalWeight = 0.0;

        // Start the MST from node 0
        minEdgeWeight[0] = 0.0;

        // Iterate to include all nodes in the MST
        for (int i = 0; i < n; i++) {
            // Find the vertex with the smallest edge weight not yet included in the MST
            int u = -1;
            double minWeight = Double.MAX_VALUE;

            for (int v = 0; v < n; v++) {
                if (!inMST[v] && minEdgeWeight[v] < minWeight) {
                    minWeight = minEdgeWeight[v];
                    u = v;
                }
            }

            if (u == -1) {
                throw new RuntimeException("Failed to construct MST");
            }

            // Include this vertex in the MST
            inMST[u] = true;
            totalWeight += minWeight;

            // Update the edge weights for the remaining vertices
            final int p = u;
            forEachEdge(u, (v, weight) -> {
                if (!inMST[v]) {
                    if (weight < minEdgeWeight[v]) {
                        minEdgeWeight[v] = weight;
                        parent[v] = p; // Update the parent for this vertex
                    }
                }
            });
        }

        if (mst != null) {
            for (int v = 1; v < n; v++) {
                if (parent[v] != -1) {
                    int u = parent[v];
                    mst.add(new Edge(v, u, minEdgeWeight[v]));
                }
            }
        }
        return totalWeight;
    }

    /**
     * Returns the distance of path.
     * @param path the path.
     * @return the distance.
     */
    public double getPathDistance(int[] path) {
        double distance = 0.0;
        for (int i = 0; i < path.length-1; i++) {
            distance += getDistance(path[i], path[i+1]);
        }
        return distance;
    }

    /**
     * Replace edges path[i]->path[i+1] and path[j]->path[j+1]
     * with path[i]->path[j] and path[i+1]->path[j+1]
     */
    private void swapEdges(int[] path, int i, int j) {
        i += 1;
        while (i < j) {
            Sort.swap(path, i, j);
            i++;
            j--;
        }
    }

    /**
     * A search node in TSP branch and bound algorithm.
     */
    private record TspNode(int[] path, double lowerBound, double cost) implements Comparable {
        /**
         * Returns the level of search tree, i.e. the length of partial path.
         * @return the level of search tree.
         */
        public int level() {
            return path.length;
        }

        @Override
        public int compareTo(TspNode o) {
            return Double.compare(lowerBound, o.lowerBound);
        }
    }

    /**
     * Returns the MST cost of vertices not in the path.
     * @param inPath the flag if a node is in the partial path.
     * @return the MST cost.
     */
    private double mstLowerBound(boolean[] inPath) {
        int n = getVertexCount();

        // Tracks whether a node is included in the MST
        boolean[] inMST = new boolean[n];
        // Stores the minimum edge weight to add a node to the MST
        double[] minEdgeWeight = new double[n];
        Arrays.fill(minEdgeWeight, Double.MAX_VALUE);

        // Total weight of the MST
        double totalWeight = 0.0;

        // Find the start node
        for (int i = 0; i < n; i++) {
            if (!inPath[i]) {
                minEdgeWeight[i] = 0.0;
                break;
            }
        }

        // Iterate to include all nodes in the MST
        for (int i = 0; i < n; i++) {
            // Find the vertex with the smallest edge weight not yet included in the MST
            int u = -1;
            double minWeight = Double.MAX_VALUE;

            for (int v = 0; v < n; v++) {
                if (!inMST[v] && !inPath[v] && minEdgeWeight[v] < minWeight) {
                    minWeight = minEdgeWeight[v];
                    u = v;
                }
            }

            if (u == -1) break; // All reachable nodes are visited

            // Include this vertex in the MST
            inMST[u] = true;
            totalWeight += minWeight;

            // Update the edge weights for the remaining vertices
            final int p = u;
            forEachEdge(u, (v, weight) -> {
                if (!inMST[v] && !inPath[v]) {
                    minEdgeWeight[v] = Math.min(minEdgeWeight[v], weight);
                }
            });
        }

        // add the edge back to node 0
        forEachEdge(0, (v, weight) -> {
            if (inMST[v]) {
                minEdgeWeight[0] = Math.min(minEdgeWeight[0], weight);
            }
        });

        return totalWeight + minEdgeWeight[0];
    }

    /**
     * Returns the optimal travelling salesman problem (TSP) tour with
     * branch and bound algorithm.
     * @return the optimal TSP tour.
     */
    public int[] tsp() {
        int n = getVertexCount();
        if (n < 2) {
            throw new UnsupportedOperationException("Cannot construct TSP with fewer than 2 vertices.");
        }

        // Initialize the best cost with nearest insertion
        int[] tour = nearestInsertion();
        double bestCost = getPathDistance(tour);

        // Push the initial node into the priority queue
        java.util.PriorityQueue pq = new java.util.PriorityQueue<>();
        pq.offer(new TspNode(new int[1], 0.0, 0.0));

        // Perform Branch and Bound with Best-First Search
        while (!pq.isEmpty()) {
            var current = pq.poll();

            // Skip nodes with bounds worse than the current best solution
            if (current.lowerBound >= bestCost) continue;

            // If we reach the last level, check the complete path
            if (current.level() == n) {
                double cost = current.cost + getDistance(current.path[n-1], 0); // Return to start
                if (cost < bestCost) {
                    bestCost = cost;
                    System.arraycopy(current.path, 0, tour, 0, n);
                }
                continue;
            }

            boolean[] inPath = new boolean[n];
            for (var node : current.path) {
                inPath[node] = true;
            }
            // The cost of checking all remaining branches is cheaper compared
            // estimating MST lower bound if there are fewer than 5 nodes.
            double mst = (n - current.level() < 5) ? 0 : mstLowerBound(inPath);

            // Explore all possible next nodes
            final double currentBest = bestCost;
            forEachEdge(current.path[current.level() - 1], (v, weight) -> {
                if (inPath[v]) return;

                double nextCost = current.cost + weight;
                double nextLowerBound = nextCost + mst;

                // Prune branches with higher bounds
                if (nextLowerBound < currentBest) {
                    int[] nextPath = Arrays.copyOfRange(current.path, 0, current.level() + 1);
                    nextPath[current.level()] = v;
                    pq.offer(new TspNode(nextPath, nextLowerBound, nextCost));
                }
            });
        }

        logger.info("Branch and bound TSP cost = {}", bestCost);
        return tour;
    }

    /**
     * Returns the optimal TSP tour with Held-Karp algorithm.
     * It works on graph up to 31 vertices.
     * @return the optimal TSP tour.
     */
    public int[] heldKarp() {
        int n = getVertexCount();
        if (n < 2) {
            throw new UnsupportedOperationException("Cannot construct TSP with fewer than 2 vertices.");
        }

        if (n > 31) {
            throw new UnsupportedOperationException("Held-Karp cannot run with more than 31 vertices.");
        }

        // DP table: dp[mask][i] stores the shortest path to visit all nodes in mask
        // ending at node i
        int p = 1 << n;
        double[][] dp = new double[p][n];
        for (var row : dp) {
            Arrays.fill(row, Double.POSITIVE_INFINITY);
        }

        // Base case: cost of starting at node 0 and visiting only node 0 is 0
        dp[1][0] = 0.0;

        // Iterate through all possible subsets of nodes (represented by bitmasks)
        for (int mask = 1; mask < p; mask++) {
            for (int u = 0; u < n; u++) {
                if ((mask & (1 << u)) == 0) continue; // u is not in the subset represented by mask
                for (int v = 0; v < n; v++) {
                    if ((mask & (1 << v)) != 0) continue; // v is already in the subset
                    int nextMask = mask | (1 << v); // Add v to the subset
                    dp[nextMask][v] = Math.min(dp[nextMask][v], dp[mask][u] + getDistance(u, v));
                }
            }
        }

        // Reconstruct the optimal tour by backtracking
        int endMask = p - 1;
        int lastNode = 0;
        double bestCost = Double.POSITIVE_INFINITY;
        // Find the last node of the optimal tour (minimum cost returning to node 0)
        for (int u = 1; u < n; u++) {
            double cost = dp[endMask][u] + getDistance(u, 0);
            if (cost < bestCost) {
                bestCost = cost;
                lastNode = u;
            }
        }

        // Backtracking to find the tour
        int mask = endMask;
        int index = 1; // The tour always starts with node 0.
        int[] tour = new int[n+1];
        while (mask != 0) {
            tour[index++] = lastNode;
            int prevMask = mask ^ (1 << lastNode);
            for (int u = 0; u < n; u++) {
                if (dp[mask][lastNode] == dp[prevMask][u] + getDistance(u, lastNode)) {
                    lastNode = u;
                    break;
                }
            }
            mask = prevMask;
        }

        logger.info("Held-Karp TSP cost = {}", bestCost);
        MathEx.reverse(tour);
        return tour;
    }

    /**
     * Returns the insertion position that causes minimum increase of TSP tour distance.
     * @param node the node to insert.
     * @param tour the tour path.
     * @param length the length of tour path.
     * @return the insertion position.
     */
    private int getInsertPosition(int node, int[] tour, int length) {
        int insertIndex = 1;
        double minIncrease = Double.MAX_VALUE;
        for (int i = 0; i < length; i++) {
            int node1 = tour[i];
            int node2 = tour[(i+1) % length];
            double increase = getDistance(node1, node) + getDistance(node, node2) - getDistance(node1, node2);
            if (increase < minIncrease) {
                minIncrease = increase;
                insertIndex = i + 1;
            }
        }
        return insertIndex;
    }

    /**
     * Returns the approximate solution to TSP with the
     * nearest insertion heuristic.
     * @return the approximate solution to TSP.
     */
    public int[] nearestInsertion() {
        int n = getVertexCount();
        if (n < 2) {
            throw new UnsupportedOperationException("Cannot construct TSP with fewer than 2 vertices.");
        }

        int[] tour = new int[n+1];
        double[] dist = new double[n];
        boolean[] visited = new boolean[n];
        visited[0] = true;
        Arrays.fill(dist, Double.POSITIVE_INFINITY);

        forEachEdge(0, (i, weight) -> dist[i] = weight);

        int nearestNode = MathEx.whichMin(dist);
        tour[1] = nearestNode;
        visited[nearestNode] = true;

        for (int length = 2; length < n; length++) {
            forEachEdge(nearestNode, (i, weight) -> {
                if (!visited[i]) {
                    dist[i] = Math.min(dist[i], weight);
                }
            });

            nearestNode = -1;
            double minDistance = Double.POSITIVE_INFINITY;
            for (int i = 0; i < n; i++) {
                if (!visited[i]) {
                    if (dist[i] < minDistance) {
                        minDistance = dist[i];
                        nearestNode = i;
                    }
                }
            }

            // insert at the position that minimizes the increase in tour length
            int insertIndex = getInsertPosition(nearestNode, tour, length);
            System.arraycopy(tour, insertIndex, tour, insertIndex+1, length-insertIndex);
            tour[insertIndex] = nearestNode;
            visited[nearestNode] = true;
        }

        return tour;
    }

    /**
     * Returns the approximate solution to TSP with the
     * farthest insertion heuristic.
     * @return the approximate solution to TSP.
     */
    public int[] farthestInsertion() {
        int n = getVertexCount();
        if (n < 2) {
            throw new UnsupportedOperationException("Cannot construct TSP with fewer than 2 vertices.");
        }

        int[] tour = new int[n+1];
        double[] dist = new double[n];
        boolean[] visited = new boolean[n];
        visited[0] = true;

        forEachEdge(0, (i, weight) -> dist[i] = weight);

        int farthestNode = MathEx.whichMax(dist);
        tour[1] = farthestNode;
        visited[farthestNode] = true;

        for (int length = 2; length < n; length++) {
            forEachEdge(farthestNode, (i, weight) -> {
                if (!visited[i]) {
                    dist[i] = Math.max(dist[i], weight);
                }
            });

            farthestNode = -1;
            double maxDistance = Double.NEGATIVE_INFINITY;
            for (int i = 0; i < n; i++) {
                if (!visited[i]) {
                    if (dist[i] > maxDistance) {
                        maxDistance = dist[i];
                        farthestNode = i;
                    }
                }
            }

            // insert at the position that minimizes the increase in tour length
            int insertIndex = getInsertPosition(farthestNode, tour, length);
            System.arraycopy(tour, insertIndex, tour, insertIndex+1, length-insertIndex);
            tour[insertIndex] = farthestNode;
            visited[farthestNode] = true;
        }

        return tour;
    }

    /**
     * Returns the approximate solution to TSP with the
     * arbitrary insertion heuristic.
     * @return the approximate solution to TSP.
     */
    public int[] arbitraryInsertion() {
        int n = getVertexCount();
        if (n < 2) {
            throw new UnsupportedOperationException("Cannot construct TSP with fewer than 2 vertices.");
        }

        int[] tour = new int[n+1];
        tour[1] = 1;

        for (int v = 2; v < n; v++) {
            // insert at the position that minimizes the increase in tour length
            int insertIndex = getInsertPosition(v, tour, v);
            System.arraycopy(tour, insertIndex, tour, insertIndex+1, v-insertIndex);
            tour[insertIndex] = v;
        }

        return tour;
    }

    /**
     * Improves an existing TSP tour with the 2-opt heuristic. The method
     * reconnects pairs of non-adjacent edges until no more pairs can be
     * swapped to further improve the solution.
     * @param tour an existing TSP tour. It may be revised with a better tour
     *             of lower cost.
     * @param maxIter the maximum number of iterations of the outer loop.
     * @return the improved tour cost.
     */
    public double opt2(int[] tour, int maxIter) {
        int n = getVertexCount();
        if (tour.length != n+1) {
            throw new IllegalArgumentException("Invalid tour length: " + tour.length);
        }

        double cost = getPathDistance(tour);
        boolean improved = true;
        for (int iter = 0; improved && iter < maxIter; iter++) {
            improved = false;
            for (int i = 0; i < n - 2; i++) {
                for (int j = i + 2; j < n; j++) {
                    double d1 = getWeight(tour[i], tour[j]);
                    double d2 = getWeight(tour[i+1], tour[(j+1)%n]);
                    if (d1 != 0 && d2 != 0) {
                        double delta = d1 + d2 - getWeight(tour[i], tour[i+1]) - getWeight(tour[j], tour[(j+1)%n]);

                        // If the length of the path is reduced, do a 2-opt swap
                        if (delta < 0) {
                            swapEdges(tour, i, j);
                            cost += delta;
                            improved = true;
                        }
                    }
                }
            }
        }
        return cost;
    }

    /**
     * Returns the approximate solution to TSP with Christofides algorithm.
     * @return the approximate solution to TSP.
     */
    public int[] christofides() {
        int n = getVertexCount();
        if (n < 2) {
            throw new UnsupportedOperationException("Cannot construct TSP with fewer than 2 vertices.");
        }

        // Step 1: Find a Minimum Spanning Tree (MST)
        List mst = new ArrayList<>();
        prim(mst);

        // Step 2: Find vertices with odd degree in MST
        int[] degree = new int[n];
        for (var edge : mst) {
            degree[edge.u()]++;
            degree[edge.v()]++;
        }

        List oddDegreeVertices = new ArrayList<>();
        for (int i = 0; i < n; ++i) {
            if (degree[i] % 2 != 0) {
                oddDegreeVertices.add(i);
            }
        }

        // Step 3: Perform Minimum Weight Perfect Matching on odd-degree vertices
        int m = oddDegreeVertices.size();
        List matching = new ArrayList<>();
        boolean[] matched = new boolean[m];

        for (int i = 0; i < m; i++) {
            if (matched[i]) continue;

            int bestPartner = -1;
            double minWeight = Double.POSITIVE_INFINITY;

            for (int j = i + 1; j < m; j++) {
                if (matched[j]) continue;

                double weight = getDistance(oddDegreeVertices.get(i), oddDegreeVertices.get(j));
                if (weight < minWeight) {
                    minWeight = weight;
                    bestPartner = j;
                }
            }

            if (bestPartner != -1) {
                matching.add(new Edge(oddDegreeVertices.get(i), oddDegreeVertices.get(bestPartner)));
                matched[i] = true;
                matched[bestPartner] = true;
            }
        }

        // Step 4: Combine MST edges and matching edges to form a multigraph
        int[][] multigraph = new int[n][matching.size() + n - 1];
        int index = 0;
        for (var edge : mst) {
            multigraph[edge.u()][index] = edge.v();
            multigraph[edge.v()][index++] = edge.u();
        }
        for (var edge : matching) {
            multigraph[edge.u()][index] = edge.v();
            multigraph[edge.v()][index++] = edge.u();
        }
        assert index == matching.size() + n - 1;

        // Step 5: Find an Eulerian Circuit
        List circuit = new ArrayList<>();
        Stack stack = new Stack<>();
        boolean[][] visited = new boolean[n][n];

        stack.push(0);
        while (!stack.isEmpty()) {
            int u = stack.peek();

            boolean found = false;
            for (var v : multigraph[u]) {
                if (!visited[u][v]) {
                    stack.push(v);
                    visited[u][v] = true;
                    visited[v][u] = true;
                    found = true;
                    break;
                }
            }

            if (!found) {
                circuit.add(u);
                stack.pop();
            }
        }

        // Step 6: Shortcut the Eulerian Circuit to form a Hamiltonian Path
        boolean[] inPath = new boolean[n];
        int[] tour = new int[n+1];

        index = 0;
        for (var v : circuit) {
            if (!inPath[v]) {
                tour[index++] = v;
                inPath[v] = true;
            }
        }
        assert index == n;

        return tour;
    }
}




© 2015 - 2025 Weber Informatics LLC | Privacy Policy