libs/graph/test/metric_tsp_approx.cpp

   1 //=======================================================================
   2 // Copyright 2008
   3 // Author: Matyas W Egyhazy
   4 //
   5 // Distributed under the Boost Software License, Version 1.0. (See
   6 // accompanying file LICENSE_1_0.txt or copy at
   7 // http://www.boost.org/LICENSE_1_0.txt)
   8 //=======================================================================
   9
  10 #include <iostream>
  11 #include <vector>
  12 #include <fstream>
  13 #include <set>
  14 #include <ctime>
  15
  16 #include <boost/assert.hpp>
  17 #include <boost/lexical_cast.hpp>
  18 #include <boost/random.hpp>
  19 #include <boost/timer.hpp>
  20 #include <boost/integer_traits.hpp>
  21 #include <boost/graph/adjacency_matrix.hpp>
  22 #include <boost/graph/adjacency_list.hpp>
  23 #include <boost/graph/simple_point.hpp>
  24 #include <boost/graph/metric_tsp_approx.hpp>
  25 #include <boost/graph/graphviz.hpp>
  26
  27 // TODO: Integrate this into the test system a little better. We need to run
  28 // the test with some kind of input file.
  29
  30 template<typename PointType>
  31 struct cmpPnt
  32 {
  33     bool operator()(const boost::simple_point<PointType>& l,
  34                     const boost::simple_point<PointType>& r) const
  35     { return (l.x > r.x); }
  36 };
  37
  38 //add edges to the graph (for each node connect it to all other nodes)
  39 template<typename VertexListGraph, typename PointContainer,
  40     typename WeightMap, typename VertexIndexMap>
  41 void connectAllEuclidean(VertexListGraph& g,
  42                         const PointContainer& points,
  43                         WeightMap wmap,            // Property maps passed by value
  44                         VertexIndexMap vmap,       // Property maps passed by value
  45                         int /*sz*/)
  46 {
  47     using namespace boost;
  48     using namespace std;
  49     typedef typename graph_traits<VertexListGraph>::edge_descriptor Edge;
  50     typedef typename graph_traits<VertexListGraph>::vertex_iterator VItr;
  51
  52     Edge e;
  53     bool inserted;
  54
  55     pair<VItr, VItr> verts(vertices(g));
  56     for (VItr src(verts.first); src != verts.second; src++)
  57     {
  58         for (VItr dest(src); dest != verts.second; dest++)
  59         {
  60             if (dest != src)
  61             {
  62                 double weight(sqrt(pow(
  63                     static_cast<double>(points[vmap[*src]].x -
  64                         points[vmap[*dest]].x), 2.0) +
  65                     pow(static_cast<double>(points[vmap[*dest]].y -
  66                         points[vmap[*src]].y), 2.0)));
  67
  68                 boost::tie(e, inserted) = add_edge(*src, *dest, g);
  69
  70                 wmap[e] = weight;
  71             }
  72
  73         }
  74
  75     }
  76 }
  77
  78 // Create a randomly generated point
  79 // scatter time execution
  80 void testScalability(unsigned numpts)
  81 {
  82     using namespace boost;
  83     using namespace std;
  84
  85     typedef adjacency_matrix<undirectedS, no_property,
  86         property <edge_weight_t, double,
  87         property<edge_index_t, int> > > Graph;
  88     typedef graph_traits<Graph>::vertex_descriptor Vertex;
  89     typedef property_map<Graph, edge_weight_t>::type WeightMap;
  90     typedef set<simple_point<double>, cmpPnt<double> > PointSet;
  91     typedef vector< Vertex > Container;
  92
  93     boost::mt19937 rng(time(0));
  94     uniform_real<> range(0.01, (numpts * 2));
  95     variate_generator<boost::mt19937&, uniform_real<> >
  96         pnt_gen(rng, range);
  97
  98     PointSet points;
  99     simple_point<double> pnt;
 100
 101     while (points.size() < numpts)
 102     {
 103         pnt.x = pnt_gen();
 104         pnt.y = pnt_gen();
 105         points.insert(pnt);
 106     }
 107
 108     Graph g(numpts);
 109     WeightMap weight_map(get(edge_weight, g));
 110     vector<simple_point<double> > point_vec(points.begin(), points.end());
 111
 112     connectAllEuclidean(g, point_vec, weight_map, get(vertex_index, g), numpts);
 113
 114     Container c;
 115     timer t;
 116     double len = 0.0;
 117
 118     // Run the TSP approx, creating the visitor on the fly.
 119     metric_tsp_approx(g, make_tsp_tour_len_visitor(g, back_inserter(c), len, weight_map));
 120
 121     cout << "Number of points: " << num_vertices(g) << endl;
 122     cout << "Number of edges: " << num_edges(g) << endl;
 123     cout << "Length of tour: " << len << endl;
 124     cout << "Elapsed: " << t.elapsed() << endl;
 125 }
 126
 127 template <typename PositionVec>
 128 void checkAdjList(PositionVec v)
 129 {
 130     using namespace std;
 131     using namespace boost;
 132
 133     typedef adjacency_list<listS, listS, undirectedS> Graph;
 134     typedef graph_traits<Graph>::vertex_descriptor Vertex;
 135     typedef graph_traits <Graph>::edge_descriptor Edge;
 136     typedef vector<Vertex> Container;
 137     typedef map<Vertex, std::size_t> VertexIndexMap;
 138     typedef map<Edge, double> EdgeWeightMap;
 139     typedef associative_property_map<VertexIndexMap> VPropertyMap;
 140     typedef associative_property_map<EdgeWeightMap> EWeightPropertyMap;
 141     typedef graph_traits<Graph>::vertex_iterator VItr;
 142
 143     Container c;
 144     EdgeWeightMap w_map;
 145     VertexIndexMap v_map;
 146     VPropertyMap v_pmap(v_map);
 147     EWeightPropertyMap w_pmap(w_map);
 148
 149     Graph g(v.size());
 150
 151     //create vertex index map
 152     VItr vi, ve;
 153     int idx(0);
 154     for (boost::tie(vi, ve) = vertices(g); vi != ve; ++vi)
 155     {
 156         Vertex v(*vi);
 157         v_pmap[v] = idx;
 158         idx++;
 159     }
 160
 161     connectAllEuclidean(g, v, w_pmap,
 162         v_pmap, v.size());
 163
 164     metric_tsp_approx_from_vertex(g,
 165         *vertices(g).first,
 166         w_pmap,
 167         v_pmap,
 168         tsp_tour_visitor<back_insert_iterator<Container > >
 169         (back_inserter(c)));
 170
 171     cout << "adj_list" << endl;
 172     for (Container::iterator itr = c.begin(); itr != c.end(); ++itr) {
 173         cout << v_map[*itr] << " ";
 174     }
 175     cout << endl << endl;
 176
 177     c.clear();
 178 }
 179
 180 static void usage()
 181 {
 182     using namespace std;
 183     cerr << "To run this program properly please place a "
 184          << "file called graph.txt"
 185          << endl << "into the current working directory." << endl
 186          << "Each line of this file should be a coordinate specifying the"
 187          << endl << "location of a vertex" << endl
 188          << "For example: " << endl << "1,2" << endl << "20,4" << endl
 189          << "15,7" << endl << endl;
 190 }
 191
 192 int main(int argc, char* argv[])
 193 {
 194    using namespace boost;
 195    using namespace std;
 196
 197     typedef vector<simple_point<double> > PositionVec;
 198     typedef adjacency_matrix<undirectedS, no_property,
 199         property <edge_weight_t, double> > Graph;
 200     typedef graph_traits<Graph>::vertex_descriptor Vertex;
 201     typedef vector<Vertex> Container;
 202     typedef property_map<Graph, edge_weight_t>::type WeightMap;
 203     typedef property_map<Graph, vertex_index_t>::type VertexMap;
 204
 205     // Make sure that the the we can parse the given file.
 206     if(argc < 2) {
 207         usage();
 208         // return -1;
 209         return 0;
 210     }
 211
 212     // Open the graph file, failing if one isn't given on the command line.
 213     ifstream fin(argv[1]);
 214     if (!fin)
 215     {
 216         usage();
 217         // return -1;
 218         return 0;
 219     }
 220
 221    string line;
 222    PositionVec position_vec;
 223
 224    int n(0);
 225    while (getline(fin, line))
 226    {
 227        simple_point<double> vertex;
 228
 229        size_t idx(line.find(","));
 230        string xStr(line.substr(0, idx));
 231        string yStr(line.substr(idx + 1, line.size() - idx));
 232
 233        vertex.x = lexical_cast<double>(xStr);
 234        vertex.y = lexical_cast<double>(yStr);
 235
 236        position_vec.push_back(vertex);
 237        n++;
 238    }
 239
 240    fin.close();
 241
 242    Container c;
 243    Graph g(position_vec.size());
 244    WeightMap weight_map(get(edge_weight, g));
 245    VertexMap v_map = get(vertex_index, g);
 246
 247    connectAllEuclidean(g, position_vec, weight_map, v_map, n);
 248
 249    metric_tsp_approx_tour(g, back_inserter(c));
 250
 251    for (vector<Vertex>::iterator itr = c.begin(); itr != c.end(); ++itr)
 252    {
 253        cout << *itr << " ";
 254    }
 255    cout << endl << endl;
 256
 257    c.clear();
 258
 259    checkAdjList(position_vec);
 260
 261    metric_tsp_approx_from_vertex(g, *vertices(g).first,
 262        get(edge_weight, g), get(vertex_index, g),
 263        tsp_tour_visitor<back_insert_iterator<vector<Vertex> > >
 264        (back_inserter(c)));
 265
 266    for (vector<Vertex>::iterator itr = c.begin(); itr != c.end(); ++itr)
 267    {
 268        cout << *itr << " ";
 269    }
 270    cout << endl << endl;
 271
 272    c.clear();
 273
 274    double len(0.0);
 275    try {
 276        metric_tsp_approx(g, make_tsp_tour_len_visitor(g, back_inserter(c), len, weight_map));
 277    }
 278    catch (const bad_graph& e) {
 279        cerr << "bad_graph: " << e.what() << endl;
 280        return -1;
 281    }
 282
 283    cout << "Number of points: " << num_vertices(g) << endl;
 284    cout << "Number of edges: " << num_edges(g) << endl;
 285    cout << "Length of Tour: " << len << endl;
 286
 287    int cnt(0);
 288    pair<Vertex,Vertex> triangleEdge;
 289    for (vector<Vertex>::iterator itr = c.begin(); itr != c.end();
 290        ++itr, ++cnt)
 291    {
 292        cout << *itr << " ";
 293
 294        if (cnt == 2)
 295        {
 296            triangleEdge.first = *itr;
 297        }
 298        if (cnt == 3)
 299        {
 300            triangleEdge.second = *itr;
 301        }
 302    }
 303    cout << endl << endl;
 304    c.clear();
 305
 306    testScalability(1000);
 307
 308    // if the graph is not fully connected then some of the
 309    // assumed triangle-inequality edges may not exist
 310    remove_edge(edge(triangleEdge.first, triangleEdge.second, g).first, g);
 311
 312     // Make sure that we can actually trap incomplete graphs.
 313     bool caught = false;
 314     try {
 315         double len = 0.0;
 316         metric_tsp_approx(g, make_tsp_tour_len_visitor(g, back_inserter(c), len, weight_map));
 317     }
 318     catch (const bad_graph& e) { caught = true; }
 319     BOOST_ASSERT(caught);
 320
 321    return 0;
 322 }