[Boost-users] tsp 2-approximation using triangle inequality

From: Matyas W Egyhazy (matt.egyhazy_at_[hidden])
Date: 2006-07-27 21:04:45

Hello,

I recently implemented an approximation for the traveling salesperson
using the boost graph library.  The interface as well as the
implementation are not polished and they require improvement, mostly due
to the work being part of a school homework assignment.  If anyone is
interested, I would like outside comments before spending (or wasting,
if no one is interested) effort to fix it up.  I have attached the .h,a
driver .cpp, and a sample input file.  The driver program will produce
output files that are viewable in a neato dot renderer.  I have tested
this only on VC++ 2005 (8.0) although I cannot think of any specific
reason why it would not compile on other compilers.

The approximation algorithm runs in polynomial time and solves the TSP
problem with the worst possible solution producing a path twice as long
as the optimal solution.  More information can be found regarding the
algorithm/heuristic on page 1028 in "Introduction to Algorithms" by
Cormen, Leiserson, Rivest, and Stein.  A slightly different description
can be found on Wikipedia:
http://en.wikipedia.org/wiki/Traveling_salesperson_problem#Triangle_inequality_and_the_Christofides_algorithm 

    My implementation matches the Introduction to Algorithms description
as I use a preorder walk of the MST without bothering to create an Euler
tour.  I believe the Euler tour part of the algorithm is only used to
prove the limit on the length of the produced solution.

Regards,

Matt


//tsp_2_approximation
//Takes a sequence of simple_points and performs 2-approximation triangle inequality to
//generate a solution for the traveling salesperson problem.
//Outputs the graphs in dot format to passed ostream objects
//Also returns the MST generated as part of the approximation and a predecessor map
//for a tour
#include <algorithm>  
#include <vector>
#include <map>
#include <string>
#include <iostream>
#include <complex>
#include <boost/graph/graph_as_tree.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/prim_minimum_spanning_tree.hpp>
#include <boost/graph/simple_point.hpp>  
#include <boost/graph/graphviz.hpp>
#include <boost/graph/iteration_macros.hpp>
#include <boost/graph/graph_traits.hpp>

template<typename Pt>
std::istream& operator >> (std::istream& is, const boost::simple_point<Pt>& p)
{
    return is >> boost::lexical_cast<std::string>(p.x) >> std::string(", ") >> 
        boost::lexical_cast<std::string>(p.y);
}

template<typename Pt>
std::ostream& operator << (std::ostream& os, const boost::simple_point<Pt>& p)
{
    return os << boost::lexical_cast<std::string>(p.x) << std::string(",") <<
        boost::lexical_cast<std::string>(p.y) << "!";
}

//g is returned as the MST
//preds is the predecessor map for the TSP approximation tour
//the PositionMap is a coordinate map of the vertices represented as simple_points
//the ostreams are optional, if they are provided the dot graphs are printed to them
template<typename Graph, typename PositionMap, typename PredecessorMap>
void tsp_2_approximation(Graph& g, const PositionMap& pos, PredecessorMap preds, 
                         std::ostream& pts_out = std::ostream(),
                         std::ostream& mst_out = std::ostream(), 
                         std::ostream& tour_out = std::ostream())
{
    using namespace boost;
    using namespace std;

    typedef graph_traits<Graph>::vertex_iterator VertexItr;
    typedef graph_traits<Graph> GraphTraits;
    typedef property_map<Graph, vertex_index_t>::type IndexMap;
    typedef graph_traits<Graph>::vertex_descriptor Vertex;
    typedef iterator_property_map<vector<Vertex>::iterator,
        property_map<Graph, vertex_index_t>::type> ParentMap;
    typedef graph_as_tree<Graph, ParentMap> Tree;
    typedef tree_traits<Tree>::node_descriptor Node;
    typedef graph_traits <Graph>::edge_descriptor Edge;
    typedef vector<simple_point<double> > PositionVec;
    typedef iterator_property_map<PositionVec::iterator, 
        property_map<Graph, vertex_index_t>::type>
        PositionMap;

    int cnt(0);
    property_map<Graph, vertex_name_t>::type name_map(get(vertex_name, g));

    //setup name map for the vertices
    IndexMap index(get(vertex_index, g));
    for (pair<VertexItr, VertexItr> vp(vertices(g)); vp.first != vp.second; 
        ++vp.first)
    {
        cnt++;
        name_map[*(vp.first)] = lexical_cast<string>(index[*(vp.first)]);
    }

    //print point layout to dot format
    if (pts_out)
    {
        dynamic_properties dp;
        dp.property("pos", pos);
        dp.property("id", get(vertex_name, g));
        write_graphviz(pts_out, g, dp, string("id"));
    }

    property_map<Graph, edge_weight_t>::type weight_map(get(edge_weight, g));
    graph_traits<Graph>::edge_iterator ei, ei_end;

    for (int q(0); q < cnt; ++q)
    {
        for (int i(q); i < cnt; ++i)
        {
            if (q != i)
            {
                //Calculate distance
                double weight(sqrt(pow(static_cast<double>(pos[q].x - 
                    pos[i].x), 2.0) + pow(static_cast<double>(pos[i].y - 
                    pos[q].y), 2.0)));

                graph_traits<Graph>::edge_descriptor e; 
                bool inserted;
                tie(e, inserted) = add_edge(q, i, g);
                weight_map[e] = weight;
            }
        }
    }

    index = get(vertex_index, g);

    //extract the mst using prim's minimum spanning tree algorithm
    vector<Vertex> p(num_vertices(g));
    prim_minimum_spanning_tree(g, &p[0]);

    Graph mst(cnt);

    vector<Vertex> parent(num_vertices(mst));
    Tree t(mst, 0, make_iterator_property_map(parent.begin(), 
        get(vertex_index, mst)));

    property_map<Graph, edge_weight_t>::type 
        mst_weight_map(get(edge_weight, mst)); 

    //build mst graph using the predecessor map from prim's mst algorithm
    for (tie(ei, ei_end) = edges(g); ei != ei_end; ++ei) 
    {
        Edge e(*ei);
        Vertex u(source(e, g));
        Vertex v(target(e, g));

        if (p[v] == u || p[u] == v)
        {
            add_edge(u, v, mst);
        }
    }

    name_map = get(vertex_name, mst);
    for (tie(ei, ei_end) = edges(mst); ei != ei_end; ++ei) 
    {
        Edge e(*ei);
        Vertex u(source(e, g));
        Vertex v(target(e, g));
        double weight(sqrt(pow(static_cast<double>(pos[u].x - 
            pos[v].x), 2.0) + pow(static_cast<double>(pos[u].y - 
            pos[v].y), 2.0)));

        mst_weight_map[e] = weight;

        preds[v] = u;
        name_map[u] = lexical_cast<string>(u);
        name_map[v] = lexical_cast<string>(v);

    }

    //print mst in dot format
    if (mst_out)
    {
        dynamic_properties dp;
        dp.property("pos", pos);
        dp.property("label", get(edge_weight, mst));
        dp.property("id", get(vertex_name, mst));
        write_graphviz(mst_out, mst, dp, string("id"));
    }

    //create tour using a preorder traversal of the mst
    property_map<Graph, edge_weight_t>::type 
        tour_weight_map(get(edge_weight, mst)); 

    PreorderTraverser<Node, Tree> vis;

    traverse_tree_ref(0, t, vis);

    const vector<Node>& v(vis.getPath());

    Graph tour(cnt);

    //print out tour in dot format
    for (vector<Node>::const_iterator itr(v.begin());
        itr != v.end() - 1; ++itr)
    {
        Node s(*itr);
        Node d(*(itr + 1));
        add_edge(s, d, tour);
    }

    add_edge(*(v.end() - 1), *(v.begin()), tour);

    name_map = get(vertex_name, tour);
    for (tie(ei, ei_end) = edges(tour); ei != ei_end; ++ei) 
    {
        Edge e(*ei);
        Vertex u(source(e, g)); 
        Vertex v(target(e, g));
        double weight(sqrt(pow(static_cast<double>(pos[u].x - pos[v].x), 2.0) +
            pow(static_cast<double>(pos[u].y - pos[v].y), 2.0)));

        tour_weight_map[e] = weight;

        name_map[u] = lexical_cast<string>(u);
        name_map[v] = lexical_cast<string>(v);

    }

    if (tour_out)
    {
        dynamic_properties dp;
        dp.property("pos", pos);
        dp.property("label", get(edge_weight, tour));
        dp.property("id", get(vertex_name, tour));
        write_graphviz(tour_out, tour, dp, string("id"));
    }

    g = mst;
}

//Recursive function that traverses a tree with a custom visttor
template <class Tree, class TreeVisitor>
void traverse_tree_ref(typename boost::tree_traits<Tree>::node_descriptor v,
                       Tree& t, TreeVisitor& visitor)
{
    visitor.preorder(v, t);
    typename boost::tree_traits<Tree>::children_iterator i, end;
    tie(i, end) = children(v, t);
    if (i != end) 
    {
        traverse_tree_ref(*i++, t, visitor);
        visitor.inorder(v, t);
        while (i != end)
        {
            traverse_tree_ref(*i++, t, visitor);
        }
    } 
    else
    {
        visitor.inorder(v, t);
    }

    visitor.postorder(v, t);
}

//Keeps track of a preorder traversal of a tree
template<typename Node, typename Tree> class PreorderTraverser 
{
private:
    std::vector<Node> path_;
public:
    void preorder(Node n, Tree& t) 
    {
        path_.push_back(n);
    }

    void inorder(Node n, Tree& t) {}

    void postorder(Node, Tree& t) {}

    const std::vector<Node>& getPath() 
    { 
        return path_;
    }
};


//Main test harness for the tsp_2_approximation solution generator
#include <iostream>
#include <utility>                                  
#include <vector>
#include <fstream>
#include <boost/lexical_cast.hpp>
#include "tsp_2_approximation.h"

int main(int,char*[])
{
    using namespace boost;
    using namespace std;

    typedef adjacency_list<vecS, vecS, directedS, 
        property <vertex_name_t, string>,
        property <edge_weight_t, double > > Graph;

    typedef vector<simple_point<double> > PositionVec;

    typedef iterator_property_map<PositionVec::iterator, 
        property_map<Graph, vertex_index_t>::type> PositionMap;

    ifstream fin("graph.txt");
    if (!fin)
    {
        cerr << "To run this program properly please place a file called graph.txt" 
            << endl << "into the current working directory." << endl
            << "Each line of this file should be a coordinate specifying the" 
            << endl << "location of a vertex" << endl
            << "For example: " << endl << "1,2" << endl << "20,4" << endl 
            << "15,7" << endl << endl << 
            "The output of the program will be three DOT formatted files in the" 
            << endl << "current working directory: " << endl
            << "pts.dot - Only displays the points as described in graph.txt" << endl << 
            "mst.dot - Displays an MST for the points in graph.txt" << endl << 
            "tour.dot - Displays the TSP 2-approximation triangle inequality tour " 
            << endl << "for the points in graph.txt" 
            << endl << endl << "These files should be viewed with a 'neato'" <<  
            "DOT graph viewer" << endl << "as it can handle forced coordinates" << endl;

        return -1;

    }

    ofstream ptsOut("pts.dot");
    ofstream mstOut("mst.dot");
    ofstream tourOut("tour.dot");
    string line;
    PositionVec position_vec;

    int n(0);
    while (getline(fin, line))
    {
        simple_point<double> vertex;

        size_t idx(line.find(","));
        string xStr(line.substr(0, idx));
        string yStr(line.substr(idx + 1, line.size() - idx));

        vertex.x = lexical_cast<double>(xStr);
        vertex.y = lexical_cast<double>(yStr);

        position_vec.push_back(vertex);
        n++;
    }

    fin.close();

    Graph g(n);
    PositionMap pos(position_vec.begin(), get(vertex_index, g));
    vector<graph_traits<Graph>::vertex_descriptor> preds(n);
    tsp_2_approximation(g, pos, &preds.front(), ptsOut, mstOut, tourOut);

    return 0;
}

1,2
2,4
10,20
45,20
15,1
17,20
20,25
30,20