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