tools/lgf-gen.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Thu, 01 Nov 2018 11:27:05 +0100
changeset 1197 f179aa1045a4
parent 654 a312f84d86c6
permissions -rw-r--r--
Suppress unused typdef warnings (#615)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2009
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 /// \ingroup tools
    20 /// \file
    21 /// \brief Special plane graph generator.
    22 ///
    23 /// Graph generator application for various types of plane graphs.
    24 ///
    25 /// See
    26 /// \code
    27 ///   lgf-gen --help
    28 /// \endcode
    29 /// for more information on the usage.
    30 
    31 #include <algorithm>
    32 #include <set>
    33 #include <ctime>
    34 #include <lemon/list_graph.h>
    35 #include <lemon/random.h>
    36 #include <lemon/dim2.h>
    37 #include <lemon/bfs.h>
    38 #include <lemon/counter.h>
    39 #include <lemon/suurballe.h>
    40 #include <lemon/graph_to_eps.h>
    41 #include <lemon/lgf_writer.h>
    42 #include <lemon/arg_parser.h>
    43 #include <lemon/euler.h>
    44 #include <lemon/math.h>
    45 #include <lemon/kruskal.h>
    46 #include <lemon/time_measure.h>
    47 
    48 using namespace lemon;
    49 
    50 typedef dim2::Point<double> Point;
    51 
    52 GRAPH_TYPEDEFS(ListGraph);
    53 
    54 bool progress=true;
    55 
    56 int N;
    57 // int girth;
    58 
    59 ListGraph g;
    60 
    61 std::vector<Node> nodes;
    62 ListGraph::NodeMap<Point> coords(g);
    63 
    64 
    65 double totalLen(){
    66   double tlen=0;
    67   for(EdgeIt e(g);e!=INVALID;++e)
    68     tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
    69   return tlen;
    70 }
    71 
    72 int tsp_impr_num=0;
    73 
    74 const double EPSILON=1e-8;
    75 bool tsp_improve(Node u, Node v)
    76 {
    77   double luv=std::sqrt((coords[v]-coords[u]).normSquare());
    78   Node u2=u;
    79   Node v2=v;
    80   do {
    81     Node n;
    82     for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e) { }
    83     u2=v2;
    84     v2=n;
    85     if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON>
    86        std::sqrt((coords[u]-coords[u2]).normSquare())+
    87        std::sqrt((coords[v]-coords[v2]).normSquare()))
    88       {
    89          g.erase(findEdge(g,u,v));
    90          g.erase(findEdge(g,u2,v2));
    91         g.addEdge(u2,u);
    92         g.addEdge(v,v2);
    93         tsp_impr_num++;
    94         return true;
    95       }
    96   } while(v2!=u);
    97   return false;
    98 }
    99 
   100 bool tsp_improve(Node u)
   101 {
   102   for(IncEdgeIt e(g,u);e!=INVALID;++e)
   103     if(tsp_improve(u,g.runningNode(e))) return true;
   104   return false;
   105 }
   106 
   107 void tsp_improve()
   108 {
   109   bool b;
   110   do {
   111     b=false;
   112     for(NodeIt n(g);n!=INVALID;++n)
   113       if(tsp_improve(n)) b=true;
   114   } while(b);
   115 }
   116 
   117 void tsp()
   118 {
   119   for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);
   120   tsp_improve();
   121 }
   122 
   123 class Line
   124 {
   125 public:
   126   Point a;
   127   Point b;
   128   Line(Point _a,Point _b) :a(_a),b(_b) {}
   129   Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}
   130   Line(const Arc &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
   131   Line(const Edge &e) : a(coords[g.u(e)]),b(coords[g.v(e)]) {}
   132 };
   133 
   134 inline std::ostream& operator<<(std::ostream &os, const Line &l)
   135 {
   136   os << l.a << "->" << l.b;
   137   return os;
   138 }
   139 
   140 bool cross(Line a, Line b)
   141 {
   142   Point ao=rot90(a.b-a.a);
   143   Point bo=rot90(b.b-b.a);
   144   return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 &&
   145     (bo*(a.a-b.a))*(bo*(a.b-b.a))<0;
   146 }
   147 
   148 struct Parc
   149 {
   150   Node a;
   151   Node b;
   152   double len;
   153 };
   154 
   155 bool pedgeLess(Parc a,Parc b)
   156 {
   157   return a.len<b.len;
   158 }
   159 
   160 std::vector<Edge> arcs;
   161 
   162 namespace _delaunay_bits {
   163 
   164   struct Part {
   165     int prev, curr, next;
   166 
   167     Part(int p, int c, int n) : prev(p), curr(c), next(n) {}
   168   };
   169 
   170   inline std::ostream& operator<<(std::ostream& os, const Part& part) {
   171     os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';
   172     return os;
   173   }
   174 
   175   inline double circle_point(const Point& p, const Point& q, const Point& r) {
   176     double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);
   177     if (a == 0) return std::numeric_limits<double>::quiet_NaN();
   178 
   179     double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +
   180       (q.x * q.x + q.y * q.y) * (r.y - p.y) +
   181       (r.x * r.x + r.y * r.y) * (p.y - q.y);
   182 
   183     double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +
   184       (q.x * q.x + q.y * q.y) * (r.x - p.x) +
   185       (r.x * r.x + r.y * r.y) * (p.x - q.x);
   186 
   187     double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +
   188       (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +
   189       (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);
   190 
   191     return d / (2 * a) + std::sqrt((d * d + e * e) / (4 * a * a) + f / a);
   192   }
   193 
   194   inline bool circle_form(const Point& p, const Point& q, const Point& r) {
   195     return rot90(q - p) * (r - q) < 0.0;
   196   }
   197 
   198   inline double intersection(const Point& p, const Point& q, double sx) {
   199     const double epsilon = 1e-8;
   200 
   201     if (p.x == q.x) return (p.y + q.y) / 2.0;
   202 
   203     if (sx < p.x + epsilon) return p.y;
   204     if (sx < q.x + epsilon) return q.y;
   205 
   206     double a = q.x - p.x;
   207     double b = (q.x - sx) * p.y - (p.x - sx) * q.y;
   208     double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();
   209     return (b - std::sqrt(d)) / a;
   210   }
   211 
   212   struct YLess {
   213 
   214 
   215     YLess(const std::vector<Point>& points, double& sweep)
   216       : _points(points), _sweep(sweep) {}
   217 
   218     bool operator()(const Part& l, const Part& r) const {
   219       const double epsilon = 1e-8;
   220 
   221       //      std::cerr << l << " vs " << r << std::endl;
   222       double lbx = l.prev != -1 ?
   223         intersection(_points[l.prev], _points[l.curr], _sweep) :
   224         - std::numeric_limits<double>::infinity();
   225       double rbx = r.prev != -1 ?
   226         intersection(_points[r.prev], _points[r.curr], _sweep) :
   227         - std::numeric_limits<double>::infinity();
   228       double lex = l.next != -1 ?
   229         intersection(_points[l.curr], _points[l.next], _sweep) :
   230         std::numeric_limits<double>::infinity();
   231       double rex = r.next != -1 ?
   232         intersection(_points[r.curr], _points[r.next], _sweep) :
   233         std::numeric_limits<double>::infinity();
   234 
   235       if (lbx > lex) std::swap(lbx, lex);
   236       if (rbx > rex) std::swap(rbx, rex);
   237 
   238       if (lex < epsilon + rex && lbx + epsilon < rex) return true;
   239       if (rex < epsilon + lex && rbx + epsilon < lex) return false;
   240       return lex < rex;
   241     }
   242 
   243     const std::vector<Point>& _points;
   244     double& _sweep;
   245   };
   246 
   247   struct BeachIt;
   248 
   249   typedef std::multimap<double, BeachIt*> SpikeHeap;
   250 
   251   typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;
   252 
   253   struct BeachIt {
   254     Beach::iterator it;
   255 
   256     BeachIt(Beach::iterator iter) : it(iter) {}
   257   };
   258 
   259 }
   260 
   261 inline void delaunay() {
   262   Counter cnt("Number of arcs added: ");
   263 
   264   using namespace _delaunay_bits;
   265 
   266   typedef _delaunay_bits::Part Part;
   267   typedef std::vector<std::pair<double, int> > SiteHeap;
   268 
   269 
   270   std::vector<Point> points;
   271   std::vector<Node> nodes;
   272 
   273   for (NodeIt it(g); it != INVALID; ++it) {
   274     nodes.push_back(it);
   275     points.push_back(coords[it]);
   276   }
   277 
   278   SiteHeap siteheap(points.size());
   279 
   280   double sweep;
   281 
   282 
   283   for (int i = 0; i < int(siteheap.size()); ++i) {
   284     siteheap[i] = std::make_pair(points[i].x, i);
   285   }
   286 
   287   std::sort(siteheap.begin(), siteheap.end());
   288   sweep = siteheap.front().first;
   289 
   290   YLess yless(points, sweep);
   291   Beach beach(yless);
   292 
   293   SpikeHeap spikeheap;
   294 
   295   std::set<std::pair<int, int> > arcs;
   296 
   297   int siteindex = 0;
   298   {
   299     SiteHeap front;
   300 
   301     while (siteindex < int(siteheap.size()) &&
   302            siteheap[0].first == siteheap[siteindex].first) {
   303       front.push_back(std::make_pair(points[siteheap[siteindex].second].y,
   304                                      siteheap[siteindex].second));
   305       ++siteindex;
   306     }
   307 
   308     std::sort(front.begin(), front.end());
   309 
   310     for (int i = 0; i < int(front.size()); ++i) {
   311       int prev = (i == 0 ? -1 : front[i - 1].second);
   312       int curr = front[i].second;
   313       int next = (i + 1 == int(front.size()) ? -1 : front[i + 1].second);
   314 
   315       beach.insert(std::make_pair(Part(prev, curr, next),
   316                                   spikeheap.end()));
   317     }
   318   }
   319 
   320   while (siteindex < int(points.size()) || !spikeheap.empty()) {
   321 
   322     SpikeHeap::iterator spit = spikeheap.begin();
   323 
   324     if (siteindex < int(points.size()) &&
   325         (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {
   326       int site = siteheap[siteindex].second;
   327       sweep = siteheap[siteindex].first;
   328 
   329       Beach::iterator bit = beach.upper_bound(Part(site, site, site));
   330 
   331       if (bit->second != spikeheap.end()) {
   332         delete bit->second->second;
   333         spikeheap.erase(bit->second);
   334       }
   335 
   336       int prev = bit->first.prev;
   337       int curr = bit->first.curr;
   338       int next = bit->first.next;
   339 
   340       beach.erase(bit);
   341 
   342       SpikeHeap::iterator pit = spikeheap.end();
   343       if (prev != -1 &&
   344           circle_form(points[prev], points[curr], points[site])) {
   345         double x = circle_point(points[prev], points[curr], points[site]);
   346         pit = spikeheap.insert(std::make_pair(x, new BeachIt(beach.end())));
   347         pit->second->it =
   348           beach.insert(std::make_pair(Part(prev, curr, site), pit));
   349       } else {
   350         beach.insert(std::make_pair(Part(prev, curr, site), pit));
   351       }
   352 
   353       beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
   354 
   355       SpikeHeap::iterator nit = spikeheap.end();
   356       if (next != -1 &&
   357           circle_form(points[site], points[curr],points[next])) {
   358         double x = circle_point(points[site], points[curr], points[next]);
   359         nit = spikeheap.insert(std::make_pair(x, new BeachIt(beach.end())));
   360         nit->second->it =
   361           beach.insert(std::make_pair(Part(site, curr, next), nit));
   362       } else {
   363         beach.insert(std::make_pair(Part(site, curr, next), nit));
   364       }
   365 
   366       ++siteindex;
   367     } else {
   368       sweep = spit->first;
   369 
   370       Beach::iterator bit = spit->second->it;
   371 
   372       int prev = bit->first.prev;
   373       int curr = bit->first.curr;
   374       int next = bit->first.next;
   375 
   376       {
   377         std::pair<int, int> arc;
   378 
   379         arc = prev < curr ?
   380           std::make_pair(prev, curr) : std::make_pair(curr, prev);
   381 
   382         if (arcs.find(arc) == arcs.end()) {
   383           arcs.insert(arc);
   384           g.addEdge(nodes[prev], nodes[curr]);
   385           ++cnt;
   386         }
   387 
   388         arc = curr < next ?
   389           std::make_pair(curr, next) : std::make_pair(next, curr);
   390 
   391         if (arcs.find(arc) == arcs.end()) {
   392           arcs.insert(arc);
   393           g.addEdge(nodes[curr], nodes[next]);
   394           ++cnt;
   395         }
   396       }
   397 
   398       Beach::iterator pbit = bit; --pbit;
   399       int ppv = pbit->first.prev;
   400       Beach::iterator nbit = bit; ++nbit;
   401       int nnt = nbit->first.next;
   402 
   403       if (bit->second != spikeheap.end())
   404         {
   405           delete bit->second->second;
   406           spikeheap.erase(bit->second);
   407         }
   408       if (pbit->second != spikeheap.end())
   409         {
   410           delete pbit->second->second;
   411           spikeheap.erase(pbit->second);
   412         }
   413       if (nbit->second != spikeheap.end())
   414         {
   415           delete nbit->second->second;
   416           spikeheap.erase(nbit->second);
   417         }
   418       
   419       beach.erase(nbit);
   420       beach.erase(bit);
   421       beach.erase(pbit);
   422 
   423       SpikeHeap::iterator pit = spikeheap.end();
   424       if (ppv != -1 && ppv != next &&
   425           circle_form(points[ppv], points[prev], points[next])) {
   426         double x = circle_point(points[ppv], points[prev], points[next]);
   427         if (x < sweep) x = sweep;
   428         pit = spikeheap.insert(std::make_pair(x, new BeachIt(beach.end())));
   429         pit->second->it =
   430           beach.insert(std::make_pair(Part(ppv, prev, next), pit));
   431       } else {
   432         beach.insert(std::make_pair(Part(ppv, prev, next), pit));
   433       }
   434 
   435       SpikeHeap::iterator nit = spikeheap.end();
   436       if (nnt != -1 && prev != nnt &&
   437           circle_form(points[prev], points[next], points[nnt])) {
   438         double x = circle_point(points[prev], points[next], points[nnt]);
   439         if (x < sweep) x = sweep;
   440         nit = spikeheap.insert(std::make_pair(x, new BeachIt(beach.end())));
   441         nit->second->it =
   442           beach.insert(std::make_pair(Part(prev, next, nnt), nit));
   443       } else {
   444         beach.insert(std::make_pair(Part(prev, next, nnt), nit));
   445       }
   446 
   447     }
   448   }
   449 
   450   for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
   451     int curr = it->first.curr;
   452     int next = it->first.next;
   453 
   454     if (next == -1) continue;
   455 
   456     std::pair<int, int> arc;
   457 
   458     arc = curr < next ?
   459       std::make_pair(curr, next) : std::make_pair(next, curr);
   460 
   461     if (arcs.find(arc) == arcs.end()) {
   462       arcs.insert(arc);
   463       g.addEdge(nodes[curr], nodes[next]);
   464       ++cnt;
   465     }
   466   }
   467 }
   468 
   469 void sparse(int d)
   470 {
   471   Counter cnt("Number of arcs removed: ");
   472   Bfs<ListGraph> bfs(g);
   473   for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
   474       ei!=arcs.rend();++ei)
   475     {
   476       Node a=g.u(*ei);
   477       Node b=g.v(*ei);
   478       g.erase(*ei);
   479       bfs.run(a,b);
   480       if(bfs.predArc(b)==INVALID || bfs.dist(b)>d)
   481         g.addEdge(a,b);
   482       else cnt++;
   483     }
   484 }
   485 
   486 void sparse2(int d)
   487 {
   488   Counter cnt("Number of arcs removed: ");
   489   for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
   490       ei!=arcs.rend();++ei)
   491     {
   492       Node a=g.u(*ei);
   493       Node b=g.v(*ei);
   494       g.erase(*ei);
   495       ConstMap<Arc,int> cegy(1);
   496       Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy);
   497       int k=sur.run(a,b,2);
   498       if(k<2 || sur.totalLength()>d)
   499         g.addEdge(a,b);
   500       else cnt++;
   501 //       else std::cout << "Remove arc " << g.id(a) << "-" << g.id(b) << '\n';
   502     }
   503 }
   504 
   505 void sparseTriangle(int d)
   506 {
   507   Counter cnt("Number of arcs added: ");
   508   std::vector<Parc> pedges;
   509   for(NodeIt n(g);n!=INVALID;++n)
   510     for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
   511       {
   512         Parc p;
   513         p.a=n;
   514         p.b=m;
   515         p.len=(coords[m]-coords[n]).normSquare();
   516         pedges.push_back(p);
   517       }
   518   std::sort(pedges.begin(),pedges.end(),pedgeLess);
   519   for(std::vector<Parc>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
   520     {
   521       Line li(pi->a,pi->b);
   522       EdgeIt e(g);
   523       for(;e!=INVALID && !cross(e,li);++e) ;
   524       Edge ne;
   525       if(e==INVALID) {
   526         ConstMap<Arc,int> cegy(1);
   527         Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy);
   528         int k=sur.run(pi->a,pi->b,2);
   529         if(k<2 || sur.totalLength()>d)
   530           {
   531             ne=g.addEdge(pi->a,pi->b);
   532             arcs.push_back(ne);
   533             cnt++;
   534           }
   535       }
   536     }
   537 }
   538 
   539 template <typename Graph, typename CoordMap>
   540 class LengthSquareMap {
   541 public:
   542   typedef typename Graph::Edge Key;
   543   typedef typename CoordMap::Value::Value Value;
   544 
   545   LengthSquareMap(const Graph& graph, const CoordMap& coords)
   546     : _graph(graph), _coords(coords) {}
   547 
   548   Value operator[](const Key& key) const {
   549     return (_coords[_graph.v(key)] -
   550             _coords[_graph.u(key)]).normSquare();
   551   }
   552 
   553 private:
   554 
   555   const Graph& _graph;
   556   const CoordMap& _coords;
   557 };
   558 
   559 void minTree() {
   560   std::vector<Parc> pedges;
   561   Timer T;
   562   std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
   563   delaunay();
   564   std::cout << T.realTime() << "s: Calculating spanning tree...\n";
   565   LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords);
   566   ListGraph::EdgeMap<bool> tree(g);
   567   kruskal(g, ls, tree);
   568   std::cout << T.realTime() << "s: Removing non tree arcs...\n";
   569   std::vector<Edge> remove;
   570   for (EdgeIt e(g); e != INVALID; ++e) {
   571     if (!tree[e]) remove.push_back(e);
   572   }
   573   for(int i = 0; i < int(remove.size()); ++i) {
   574     g.erase(remove[i]);
   575   }
   576   std::cout << T.realTime() << "s: Done\n";
   577 }
   578 
   579 void tsp2()
   580 {
   581   std::cout << "Find a tree..." << std::endl;
   582 
   583   minTree();
   584 
   585   std::cout << "Total arc length (tree) : " << totalLen() << std::endl;
   586 
   587   std::cout << "Make it Euler..." << std::endl;
   588 
   589   {
   590     std::vector<Node> leafs;
   591     for(NodeIt n(g);n!=INVALID;++n)
   592       if(countIncEdges(g,n)%2==1) leafs.push_back(n);
   593 
   594 //    for(unsigned int i=0;i<leafs.size();i+=2)
   595 //       g.addArc(leafs[i],leafs[i+1]);
   596 
   597     std::vector<Parc> pedges;
   598     for(unsigned int i=0;i<leafs.size()-1;i++)
   599       for(unsigned int j=i+1;j<leafs.size();j++)
   600         {
   601           Node n=leafs[i];
   602           Node m=leafs[j];
   603           Parc p;
   604           p.a=n;
   605           p.b=m;
   606           p.len=(coords[m]-coords[n]).normSquare();
   607           pedges.push_back(p);
   608         }
   609     std::sort(pedges.begin(),pedges.end(),pedgeLess);
   610     for(unsigned int i=0;i<pedges.size();i++)
   611       if(countIncEdges(g,pedges[i].a)%2 &&
   612          countIncEdges(g,pedges[i].b)%2)
   613         g.addEdge(pedges[i].a,pedges[i].b);
   614   }
   615 
   616   for(NodeIt n(g);n!=INVALID;++n)
   617     if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
   618       std::cout << "GEBASZ!!!" << std::endl;
   619 
   620   for(EdgeIt e(g);e!=INVALID;++e)
   621     if(g.u(e)==g.v(e))
   622       std::cout << "LOOP GEBASZ!!!" << std::endl;
   623 
   624   std::cout << "Number of arcs : " << countEdges(g) << std::endl;
   625 
   626   std::cout << "Total arc length (euler) : " << totalLen() << std::endl;
   627 
   628   ListGraph::EdgeMap<Arc> enext(g);
   629   {
   630     EulerIt<ListGraph> e(g);
   631     Arc eo=e;
   632     Arc ef=e;
   633 //     std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
   634     for(++e;e!=INVALID;++e)
   635       {
   636 //         std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
   637         enext[eo]=e;
   638         eo=e;
   639       }
   640     enext[eo]=ef;
   641   }
   642 
   643   std::cout << "Creating a tour from that..." << std::endl;
   644 
   645   int nnum = countNodes(g);
   646   int ednum = countEdges(g);
   647 
   648   for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p])
   649     {
   650 //       std::cout << "Checking arc " << g.id(p) << std::endl;
   651       Arc e=enext[p];
   652       Arc f=enext[e];
   653       Node n2=g.source(f);
   654       Node n1=g.oppositeNode(n2,e);
   655       Node n3=g.oppositeNode(n2,f);
   656       if(countIncEdges(g,n2)>2)
   657         {
   658 //           std::cout << "Remove an Arc" << std::endl;
   659           Arc ff=enext[f];
   660           g.erase(e);
   661           g.erase(f);
   662           if(n1!=n3)
   663             {
   664               Arc ne=g.direct(g.addEdge(n1,n3),n1);
   665               enext[p]=ne;
   666               enext[ne]=ff;
   667               ednum--;
   668             }
   669           else {
   670             enext[p]=ff;
   671             ednum-=2;
   672           }
   673         }
   674     }
   675 
   676   std::cout << "Total arc length (tour) : " << totalLen() << std::endl;
   677 
   678   std::cout << "2-opt the tour..." << std::endl;
   679 
   680   tsp_improve();
   681 
   682   std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl;
   683 }
   684 
   685 
   686 int main(int argc,const char **argv)
   687 {
   688   ArgParser ap(argc,argv);
   689 
   690 //   bool eps;
   691   bool disc_d, square_d, gauss_d;
   692 //   bool tsp_a,two_a,tree_a;
   693   int num_of_cities=1;
   694   double area=1;
   695   N=100;
   696 //   girth=10;
   697   std::string ndist("disc");
   698   ap.refOption("n", "Number of nodes (default is 100)", N)
   699     .intOption("g", "Girth parameter (default is 10)", 10)
   700     .refOption("cities", "Number of cities (default is 1)", num_of_cities)
   701     .refOption("area", "Full relative area of the cities (default is 1)", area)
   702     .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",
   703                disc_d)
   704     .optionGroup("dist", "disc")
   705     .refOption("square", "Nodes are evenly distributed on a unit square",
   706                square_d)
   707     .optionGroup("dist", "square")
   708     .refOption("gauss", "Nodes are located according to a two-dim Gauss "
   709                "distribution", gauss_d)
   710     .optionGroup("dist", "gauss")
   711     .onlyOneGroup("dist")
   712     .boolOption("eps", "Also generate .eps output (<prefix>.eps)")
   713     .boolOption("nonodes", "Draw only the edges in the generated .eps output")
   714     .boolOption("dir", "Directed graph is generated (each edge is replaced by "
   715                 "two directed arcs)")
   716     .boolOption("2con", "Create a two connected planar graph")
   717     .optionGroup("alg","2con")
   718     .boolOption("tree", "Create a min. cost spanning tree")
   719     .optionGroup("alg","tree")
   720     .boolOption("tsp", "Create a TSP tour")
   721     .optionGroup("alg","tsp")
   722     .boolOption("tsp2", "Create a TSP tour (tree based)")
   723     .optionGroup("alg","tsp2")
   724     .boolOption("dela", "Delaunay triangulation graph")
   725     .optionGroup("alg","dela")
   726     .onlyOneGroup("alg")
   727     .boolOption("rand", "Use time seed for random number generator")
   728     .optionGroup("rand", "rand")
   729     .intOption("seed", "Random seed", -1)
   730     .optionGroup("rand", "seed")
   731     .onlyOneGroup("rand")
   732     .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
   733     .run();
   734 
   735   if (ap["rand"]) {
   736     int seed = int(time(0));
   737     std::cout << "Random number seed: " << seed << std::endl;
   738     rnd = Random(seed);
   739   }
   740   if (ap.given("seed")) {
   741     int seed = ap["seed"];
   742     std::cout << "Random number seed: " << seed << std::endl;
   743     rnd = Random(seed);
   744   }
   745 
   746   std::string prefix;
   747   switch(ap.files().size())
   748     {
   749     case 0:
   750       prefix="lgf-gen-out";
   751       break;
   752     case 1:
   753       prefix=ap.files()[0];
   754       break;
   755     default:
   756       std::cerr << "\nAt most one prefix can be given\n\n";
   757       exit(1);
   758     }
   759 
   760   double sum_sizes=0;
   761   std::vector<double> sizes;
   762   std::vector<double> cum_sizes;
   763   for(int s=0;s<num_of_cities;s++)
   764     {
   765       //         sum_sizes+=rnd.exponential();
   766       double d=rnd();
   767       sum_sizes+=d;
   768       sizes.push_back(d);
   769       cum_sizes.push_back(sum_sizes);
   770     }
   771   int i=0;
   772   for(int s=0;s<num_of_cities;s++)
   773     {
   774       Point center=(num_of_cities==1?Point(0,0):rnd.disc());
   775       if(gauss_d)
   776         for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
   777           Node n=g.addNode();
   778           nodes.push_back(n);
   779           coords[n]=center+rnd.gauss2()*area*
   780             std::sqrt(sizes[s]/sum_sizes);
   781         }
   782       else if(square_d)
   783         for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
   784           Node n=g.addNode();
   785           nodes.push_back(n);
   786           coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
   787             std::sqrt(sizes[s]/sum_sizes);
   788         }
   789       else if(disc_d || true)
   790         for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
   791           Node n=g.addNode();
   792           nodes.push_back(n);
   793           coords[n]=center+rnd.disc()*area*
   794             std::sqrt(sizes[s]/sum_sizes);
   795         }
   796     }
   797 
   798 //   for (ListGraph::NodeIt n(g); n != INVALID; ++n) {
   799 //     std::cerr << coords[n] << std::endl;
   800 //   }
   801 
   802   if(ap["tsp"]) {
   803     tsp();
   804     std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
   805   }
   806   if(ap["tsp2"]) {
   807     tsp2();
   808     std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
   809   }
   810   else if(ap["2con"]) {
   811     std::cout << "Make triangles\n";
   812     //   triangle();
   813     sparseTriangle(ap["g"]);
   814     std::cout << "Make it sparser\n";
   815     sparse2(ap["g"]);
   816   }
   817   else if(ap["tree"]) {
   818     minTree();
   819   }
   820   else if(ap["dela"]) {
   821     delaunay();
   822   }
   823 
   824 
   825   std::cout << "Number of nodes    : " << countNodes(g) << std::endl;
   826   std::cout << "Number of arcs    : " << countEdges(g) << std::endl;
   827   double tlen=0;
   828   for(EdgeIt e(g);e!=INVALID;++e)
   829     tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
   830   std::cout << "Total arc length  : " << tlen << std::endl;
   831 
   832   if(ap["eps"])
   833     graphToEps(g,prefix+".eps").scaleToA4().
   834       scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false).
   835       coords(coords).hideNodes(ap.given("nonodes")).run();
   836 
   837   if(ap["dir"])
   838     DigraphWriter<ListGraph>(g,prefix+".lgf").
   839       nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
   840       nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
   841       run();
   842   else GraphWriter<ListGraph>(g,prefix+".lgf").
   843          nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
   844          nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
   845          run();
   846 }
   847