1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/christofides_tsp.h Wed Oct 17 19:14:07 2018 +0200
1.3 @@ -0,0 +1,254 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2013
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_CHRISTOFIDES_TSP_H
1.23 +#define LEMON_CHRISTOFIDES_TSP_H
1.24 +
1.25 +/// \ingroup tsp
1.26 +/// \file
1.27 +/// \brief Christofides algorithm for symmetric TSP
1.28 +
1.29 +#include <lemon/full_graph.h>
1.30 +#include <lemon/smart_graph.h>
1.31 +#include <lemon/kruskal.h>
1.32 +#include <lemon/matching.h>
1.33 +#include <lemon/euler.h>
1.34 +
1.35 +namespace lemon {
1.36 +
1.37 + /// \ingroup tsp
1.38 + ///
1.39 + /// \brief Christofides algorithm for symmetric TSP.
1.40 + ///
1.41 + /// ChristofidesTsp implements Christofides' heuristic for solving
1.42 + /// symmetric \ref tsp "TSP".
1.43 + ///
1.44 + /// This a well-known approximation method for the TSP problem with
1.45 + /// metric cost function.
1.46 + /// It has a guaranteed approximation factor of 3/2 (i.e. it finds a tour
1.47 + /// whose total cost is at most 3/2 of the optimum), but it usually
1.48 + /// provides better solutions in practice.
1.49 + /// This implementation runs in O(n<sup>3</sup>log(n)) time.
1.50 + ///
1.51 + /// The algorithm starts with a \ref spantree "minimum cost spanning tree" and
1.52 + /// finds a \ref MaxWeightedPerfectMatching "minimum cost perfect matching"
1.53 + /// in the subgraph induced by the nodes that have odd degree in the
1.54 + /// spanning tree.
1.55 + /// Finally, it constructs the tour from the \ref EulerIt "Euler traversal"
1.56 + /// of the union of the spanning tree and the matching.
1.57 + /// During this last step, the algorithm simply skips the visited nodes
1.58 + /// (i.e. creates shortcuts) assuming that the triangle inequality holds
1.59 + /// for the cost function.
1.60 + ///
1.61 + /// \tparam CM Type of the cost map.
1.62 + ///
1.63 + /// \warning CM::Value must be a signed number type.
1.64 + template <typename CM>
1.65 + class ChristofidesTsp
1.66 + {
1.67 + public:
1.68 +
1.69 + /// Type of the cost map
1.70 + typedef CM CostMap;
1.71 + /// Type of the edge costs
1.72 + typedef typename CM::Value Cost;
1.73 +
1.74 + private:
1.75 +
1.76 + GRAPH_TYPEDEFS(FullGraph);
1.77 +
1.78 + const FullGraph &_gr;
1.79 + const CostMap &_cost;
1.80 + std::vector<Node> _path;
1.81 + Cost _sum;
1.82 +
1.83 + public:
1.84 +
1.85 + /// \brief Constructor
1.86 + ///
1.87 + /// Constructor.
1.88 + /// \param gr The \ref FullGraph "full graph" the algorithm runs on.
1.89 + /// \param cost The cost map.
1.90 + ChristofidesTsp(const FullGraph &gr, const CostMap &cost)
1.91 + : _gr(gr), _cost(cost) {}
1.92 +
1.93 + /// \name Execution Control
1.94 + /// @{
1.95 +
1.96 + /// \brief Runs the algorithm.
1.97 + ///
1.98 + /// This function runs the algorithm.
1.99 + ///
1.100 + /// \return The total cost of the found tour.
1.101 + Cost run() {
1.102 + _path.clear();
1.103 +
1.104 + if (_gr.nodeNum() == 0) return _sum = 0;
1.105 + else if (_gr.nodeNum() == 1) {
1.106 + _path.push_back(_gr(0));
1.107 + return _sum = 0;
1.108 + }
1.109 + else if (_gr.nodeNum() == 2) {
1.110 + _path.push_back(_gr(0));
1.111 + _path.push_back(_gr(1));
1.112 + return _sum = 2 * _cost[_gr.edge(_gr(0), _gr(1))];
1.113 + }
1.114 +
1.115 + // Compute min. cost spanning tree
1.116 + std::vector<Edge> tree;
1.117 + kruskal(_gr, _cost, std::back_inserter(tree));
1.118 +
1.119 + FullGraph::NodeMap<int> deg(_gr, 0);
1.120 + for (int i = 0; i != int(tree.size()); ++i) {
1.121 + Edge e = tree[i];
1.122 + ++deg[_gr.u(e)];
1.123 + ++deg[_gr.v(e)];
1.124 + }
1.125 +
1.126 + // Copy the induced subgraph of odd nodes
1.127 + std::vector<Node> odd_nodes;
1.128 + for (NodeIt u(_gr); u != INVALID; ++u) {
1.129 + if (deg[u] % 2 == 1) odd_nodes.push_back(u);
1.130 + }
1.131 +
1.132 + SmartGraph sgr;
1.133 + SmartGraph::EdgeMap<Cost> scost(sgr);
1.134 + for (int i = 0; i != int(odd_nodes.size()); ++i) {
1.135 + sgr.addNode();
1.136 + }
1.137 + for (int i = 0; i != int(odd_nodes.size()); ++i) {
1.138 + for (int j = 0; j != int(odd_nodes.size()); ++j) {
1.139 + if (j == i) continue;
1.140 + SmartGraph::Edge e =
1.141 + sgr.addEdge(sgr.nodeFromId(i), sgr.nodeFromId(j));
1.142 + scost[e] = -_cost[_gr.edge(odd_nodes[i], odd_nodes[j])];
1.143 + }
1.144 + }
1.145 +
1.146 + // Compute min. cost perfect matching
1.147 + MaxWeightedPerfectMatching<SmartGraph, SmartGraph::EdgeMap<Cost> >
1.148 + mwpm(sgr, scost);
1.149 + mwpm.run();
1.150 +
1.151 + for (SmartGraph::EdgeIt e(sgr); e != INVALID; ++e) {
1.152 + if (mwpm.matching(e)) {
1.153 + tree.push_back( _gr.edge(odd_nodes[sgr.id(sgr.u(e))],
1.154 + odd_nodes[sgr.id(sgr.v(e))]) );
1.155 + }
1.156 + }
1.157 +
1.158 + // Join the spanning tree and the matching
1.159 + sgr.clear();
1.160 + for (int i = 0; i != _gr.nodeNum(); ++i) {
1.161 + sgr.addNode();
1.162 + }
1.163 + for (int i = 0; i != int(tree.size()); ++i) {
1.164 + int ui = _gr.id(_gr.u(tree[i])),
1.165 + vi = _gr.id(_gr.v(tree[i]));
1.166 + sgr.addEdge(sgr.nodeFromId(ui), sgr.nodeFromId(vi));
1.167 + }
1.168 +
1.169 + // Compute the tour from the Euler traversal
1.170 + SmartGraph::NodeMap<bool> visited(sgr, false);
1.171 + for (EulerIt<SmartGraph> e(sgr); e != INVALID; ++e) {
1.172 + SmartGraph::Node n = sgr.target(e);
1.173 + if (!visited[n]) {
1.174 + _path.push_back(_gr(sgr.id(n)));
1.175 + visited[n] = true;
1.176 + }
1.177 + }
1.178 +
1.179 + _sum = _cost[_gr.edge(_path.back(), _path.front())];
1.180 + for (int i = 0; i < int(_path.size())-1; ++i) {
1.181 + _sum += _cost[_gr.edge(_path[i], _path[i+1])];
1.182 + }
1.183 +
1.184 + return _sum;
1.185 + }
1.186 +
1.187 + /// @}
1.188 +
1.189 + /// \name Query Functions
1.190 + /// @{
1.191 +
1.192 + /// \brief The total cost of the found tour.
1.193 + ///
1.194 + /// This function returns the total cost of the found tour.
1.195 + ///
1.196 + /// \pre run() must be called before using this function.
1.197 + Cost tourCost() const {
1.198 + return _sum;
1.199 + }
1.200 +
1.201 + /// \brief Returns a const reference to the node sequence of the
1.202 + /// found tour.
1.203 + ///
1.204 + /// This function returns a const reference to a vector
1.205 + /// that stores the node sequence of the found tour.
1.206 + ///
1.207 + /// \pre run() must be called before using this function.
1.208 + const std::vector<Node>& tourNodes() const {
1.209 + return _path;
1.210 + }
1.211 +
1.212 + /// \brief Gives back the node sequence of the found tour.
1.213 + ///
1.214 + /// This function copies the node sequence of the found tour into
1.215 + /// an STL container through the given output iterator. The
1.216 + /// <tt>value_type</tt> of the container must be <tt>FullGraph::Node</tt>.
1.217 + /// For example,
1.218 + /// \code
1.219 + /// std::vector<FullGraph::Node> nodes(countNodes(graph));
1.220 + /// tsp.tourNodes(nodes.begin());
1.221 + /// \endcode
1.222 + /// or
1.223 + /// \code
1.224 + /// std::list<FullGraph::Node> nodes;
1.225 + /// tsp.tourNodes(std::back_inserter(nodes));
1.226 + /// \endcode
1.227 + ///
1.228 + /// \pre run() must be called before using this function.
1.229 + template <typename Iterator>
1.230 + void tourNodes(Iterator out) const {
1.231 + std::copy(_path.begin(), _path.end(), out);
1.232 + }
1.233 +
1.234 + /// \brief Gives back the found tour as a path.
1.235 + ///
1.236 + /// This function copies the found tour as a list of arcs/edges into
1.237 + /// the given \ref lemon::concepts::Path "path structure".
1.238 + ///
1.239 + /// \pre run() must be called before using this function.
1.240 + template <typename Path>
1.241 + void tour(Path &path) const {
1.242 + path.clear();
1.243 + for (int i = 0; i < int(_path.size()) - 1; ++i) {
1.244 + path.addBack(_gr.arc(_path[i], _path[i+1]));
1.245 + }
1.246 + if (int(_path.size()) >= 2) {
1.247 + path.addBack(_gr.arc(_path.back(), _path.front()));
1.248 + }
1.249 + }
1.250 +
1.251 + /// @}
1.252 +
1.253 + };
1.254 +
1.255 +}; // namespace lemon
1.256 +
1.257 +#endif