tools/lgf-gen.cc
author Alpar Juttner <alpar@cs.elte.hu>
Thu, 23 Apr 2009 10:44:35 +0100
changeset 659 0c8e5c688440
parent 631 33c6b6e755cd
child 663 f2d6d3446adf
permissions -rw-r--r--
Fix usage of sqrt() (#268)
     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 /// \code
    27 ///   lgf-gen --help
    28 /// \endcode
    29 /// for more info 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         spikeheap.erase(bit->second);
   333       }
   334 
   335       int prev = bit->first.prev;
   336       int curr = bit->first.curr;
   337       int next = bit->first.next;
   338 
   339       beach.erase(bit);
   340 
   341       SpikeHeap::iterator pit = spikeheap.end();
   342       if (prev != -1 &&
   343           circle_form(points[prev], points[curr], points[site])) {
   344         double x = circle_point(points[prev], points[curr], points[site]);
   345         pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
   346         pit->second.it =
   347           beach.insert(std::make_pair(Part(prev, curr, site), pit));
   348       } else {
   349         beach.insert(std::make_pair(Part(prev, curr, site), pit));
   350       }
   351 
   352       beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
   353 
   354       SpikeHeap::iterator nit = spikeheap.end();
   355       if (next != -1 &&
   356           circle_form(points[site], points[curr],points[next])) {
   357         double x = circle_point(points[site], points[curr], points[next]);
   358         nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
   359         nit->second.it =
   360           beach.insert(std::make_pair(Part(site, curr, next), nit));
   361       } else {
   362         beach.insert(std::make_pair(Part(site, curr, next), nit));
   363       }
   364 
   365       ++siteindex;
   366     } else {
   367       sweep = spit->first;
   368 
   369       Beach::iterator bit = spit->second.it;
   370 
   371       int prev = bit->first.prev;
   372       int curr = bit->first.curr;
   373       int next = bit->first.next;
   374 
   375       {
   376         std::pair<int, int> arc;
   377 
   378         arc = prev < curr ?
   379           std::make_pair(prev, curr) : std::make_pair(curr, prev);
   380 
   381         if (arcs.find(arc) == arcs.end()) {
   382           arcs.insert(arc);
   383           g.addEdge(nodes[prev], nodes[curr]);
   384           ++cnt;
   385         }
   386 
   387         arc = curr < next ?
   388           std::make_pair(curr, next) : std::make_pair(next, curr);
   389 
   390         if (arcs.find(arc) == arcs.end()) {
   391           arcs.insert(arc);
   392           g.addEdge(nodes[curr], nodes[next]);
   393           ++cnt;
   394         }
   395       }
   396 
   397       Beach::iterator pbit = bit; --pbit;
   398       int ppv = pbit->first.prev;
   399       Beach::iterator nbit = bit; ++nbit;
   400       int nnt = nbit->first.next;
   401 
   402       if (bit->second != spikeheap.end()) spikeheap.erase(bit->second);
   403       if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second);
   404       if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second);
   405 
   406       beach.erase(nbit);
   407       beach.erase(bit);
   408       beach.erase(pbit);
   409 
   410       SpikeHeap::iterator pit = spikeheap.end();
   411       if (ppv != -1 && ppv != next &&
   412           circle_form(points[ppv], points[prev], points[next])) {
   413         double x = circle_point(points[ppv], points[prev], points[next]);
   414         if (x < sweep) x = sweep;
   415         pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
   416         pit->second.it =
   417           beach.insert(std::make_pair(Part(ppv, prev, next), pit));
   418       } else {
   419         beach.insert(std::make_pair(Part(ppv, prev, next), pit));
   420       }
   421 
   422       SpikeHeap::iterator nit = spikeheap.end();
   423       if (nnt != -1 && prev != nnt &&
   424           circle_form(points[prev], points[next], points[nnt])) {
   425         double x = circle_point(points[prev], points[next], points[nnt]);
   426         if (x < sweep) x = sweep;
   427         nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
   428         nit->second.it =
   429           beach.insert(std::make_pair(Part(prev, next, nnt), nit));
   430       } else {
   431         beach.insert(std::make_pair(Part(prev, next, nnt), nit));
   432       }
   433 
   434     }
   435   }
   436 
   437   for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
   438     int curr = it->first.curr;
   439     int next = it->first.next;
   440 
   441     if (next == -1) continue;
   442 
   443     std::pair<int, int> arc;
   444 
   445     arc = curr < next ?
   446       std::make_pair(curr, next) : std::make_pair(next, curr);
   447 
   448     if (arcs.find(arc) == arcs.end()) {
   449       arcs.insert(arc);
   450       g.addEdge(nodes[curr], nodes[next]);
   451       ++cnt;
   452     }
   453   }
   454 }
   455 
   456 void sparse(int d)
   457 {
   458   Counter cnt("Number of arcs removed: ");
   459   Bfs<ListGraph> bfs(g);
   460   for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
   461       ei!=arcs.rend();++ei)
   462     {
   463       Node a=g.u(*ei);
   464       Node b=g.v(*ei);
   465       g.erase(*ei);
   466       bfs.run(a,b);
   467       if(bfs.predArc(b)==INVALID || bfs.dist(b)>d)
   468         g.addEdge(a,b);
   469       else cnt++;
   470     }
   471 }
   472 
   473 void sparse2(int d)
   474 {
   475   Counter cnt("Number of arcs removed: ");
   476   for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
   477       ei!=arcs.rend();++ei)
   478     {
   479       Node a=g.u(*ei);
   480       Node b=g.v(*ei);
   481       g.erase(*ei);
   482       ConstMap<Arc,int> cegy(1);
   483       Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy,a,b);
   484       int k=sur.run(2);
   485       if(k<2 || sur.totalLength()>d)
   486         g.addEdge(a,b);
   487       else cnt++;
   488 //       else std::cout << "Remove arc " << g.id(a) << "-" << g.id(b) << '\n';
   489     }
   490 }
   491 
   492 void sparseTriangle(int d)
   493 {
   494   Counter cnt("Number of arcs added: ");
   495   std::vector<Parc> pedges;
   496   for(NodeIt n(g);n!=INVALID;++n)
   497     for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
   498       {
   499         Parc p;
   500         p.a=n;
   501         p.b=m;
   502         p.len=(coords[m]-coords[n]).normSquare();
   503         pedges.push_back(p);
   504       }
   505   std::sort(pedges.begin(),pedges.end(),pedgeLess);
   506   for(std::vector<Parc>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
   507     {
   508       Line li(pi->a,pi->b);
   509       EdgeIt e(g);
   510       for(;e!=INVALID && !cross(e,li);++e) ;
   511       Edge ne;
   512       if(e==INVALID) {
   513         ConstMap<Arc,int> cegy(1);
   514         Suurballe<ListGraph,ConstMap<Arc,int> >
   515           sur(g,cegy,pi->a,pi->b);
   516         int k=sur.run(2);
   517         if(k<2 || sur.totalLength()>d)
   518           {
   519             ne=g.addEdge(pi->a,pi->b);
   520             arcs.push_back(ne);
   521             cnt++;
   522           }
   523       }
   524     }
   525 }
   526 
   527 template <typename Graph, typename CoordMap>
   528 class LengthSquareMap {
   529 public:
   530   typedef typename Graph::Edge Key;
   531   typedef typename CoordMap::Value::Value Value;
   532 
   533   LengthSquareMap(const Graph& graph, const CoordMap& coords)
   534     : _graph(graph), _coords(coords) {}
   535 
   536   Value operator[](const Key& key) const {
   537     return (_coords[_graph.v(key)] -
   538             _coords[_graph.u(key)]).normSquare();
   539   }
   540 
   541 private:
   542 
   543   const Graph& _graph;
   544   const CoordMap& _coords;
   545 };
   546 
   547 void minTree() {
   548   std::vector<Parc> pedges;
   549   Timer T;
   550   std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
   551   delaunay();
   552   std::cout << T.realTime() << "s: Calculating spanning tree...\n";
   553   LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords);
   554   ListGraph::EdgeMap<bool> tree(g);
   555   kruskal(g, ls, tree);
   556   std::cout << T.realTime() << "s: Removing non tree arcs...\n";
   557   std::vector<Edge> remove;
   558   for (EdgeIt e(g); e != INVALID; ++e) {
   559     if (!tree[e]) remove.push_back(e);
   560   }
   561   for(int i = 0; i < int(remove.size()); ++i) {
   562     g.erase(remove[i]);
   563   }
   564   std::cout << T.realTime() << "s: Done\n";
   565 }
   566 
   567 void tsp2()
   568 {
   569   std::cout << "Find a tree..." << std::endl;
   570 
   571   minTree();
   572 
   573   std::cout << "Total arc length (tree) : " << totalLen() << std::endl;
   574 
   575   std::cout << "Make it Euler..." << std::endl;
   576 
   577   {
   578     std::vector<Node> leafs;
   579     for(NodeIt n(g);n!=INVALID;++n)
   580       if(countIncEdges(g,n)%2==1) leafs.push_back(n);
   581 
   582 //    for(unsigned int i=0;i<leafs.size();i+=2)
   583 //       g.addArc(leafs[i],leafs[i+1]);
   584 
   585     std::vector<Parc> pedges;
   586     for(unsigned int i=0;i<leafs.size()-1;i++)
   587       for(unsigned int j=i+1;j<leafs.size();j++)
   588         {
   589           Node n=leafs[i];
   590           Node m=leafs[j];
   591           Parc p;
   592           p.a=n;
   593           p.b=m;
   594           p.len=(coords[m]-coords[n]).normSquare();
   595           pedges.push_back(p);
   596         }
   597     std::sort(pedges.begin(),pedges.end(),pedgeLess);
   598     for(unsigned int i=0;i<pedges.size();i++)
   599       if(countIncEdges(g,pedges[i].a)%2 &&
   600          countIncEdges(g,pedges[i].b)%2)
   601         g.addEdge(pedges[i].a,pedges[i].b);
   602   }
   603 
   604   for(NodeIt n(g);n!=INVALID;++n)
   605     if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
   606       std::cout << "GEBASZ!!!" << std::endl;
   607 
   608   for(EdgeIt e(g);e!=INVALID;++e)
   609     if(g.u(e)==g.v(e))
   610       std::cout << "LOOP GEBASZ!!!" << std::endl;
   611 
   612   std::cout << "Number of arcs : " << countEdges(g) << std::endl;
   613 
   614   std::cout << "Total arc length (euler) : " << totalLen() << std::endl;
   615 
   616   ListGraph::EdgeMap<Arc> enext(g);
   617   {
   618     EulerIt<ListGraph> e(g);
   619     Arc eo=e;
   620     Arc ef=e;
   621 //     std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
   622     for(++e;e!=INVALID;++e)
   623       {
   624 //         std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
   625         enext[eo]=e;
   626         eo=e;
   627       }
   628     enext[eo]=ef;
   629   }
   630 
   631   std::cout << "Creating a tour from that..." << std::endl;
   632 
   633   int nnum = countNodes(g);
   634   int ednum = countEdges(g);
   635 
   636   for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p])
   637     {
   638 //       std::cout << "Checking arc " << g.id(p) << std::endl;
   639       Arc e=enext[p];
   640       Arc f=enext[e];
   641       Node n2=g.source(f);
   642       Node n1=g.oppositeNode(n2,e);
   643       Node n3=g.oppositeNode(n2,f);
   644       if(countIncEdges(g,n2)>2)
   645         {
   646 //           std::cout << "Remove an Arc" << std::endl;
   647           Arc ff=enext[f];
   648           g.erase(e);
   649           g.erase(f);
   650           if(n1!=n3)
   651             {
   652               Arc ne=g.direct(g.addEdge(n1,n3),n1);
   653               enext[p]=ne;
   654               enext[ne]=ff;
   655               ednum--;
   656             }
   657           else {
   658             enext[p]=ff;
   659             ednum-=2;
   660           }
   661         }
   662     }
   663 
   664   std::cout << "Total arc length (tour) : " << totalLen() << std::endl;
   665 
   666   std::cout << "2-opt the tour..." << std::endl;
   667 
   668   tsp_improve();
   669 
   670   std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl;
   671 }
   672 
   673 
   674 int main(int argc,const char **argv)
   675 {
   676   ArgParser ap(argc,argv);
   677 
   678 //   bool eps;
   679   bool disc_d, square_d, gauss_d;
   680 //   bool tsp_a,two_a,tree_a;
   681   int num_of_cities=1;
   682   double area=1;
   683   N=100;
   684 //   girth=10;
   685   std::string ndist("disc");
   686   ap.refOption("n", "Number of nodes (default is 100)", N)
   687     .intOption("g", "Girth parameter (default is 10)", 10)
   688     .refOption("cities", "Number of cities (default is 1)", num_of_cities)
   689     .refOption("area", "Full relative area of the cities (default is 1)", area)
   690     .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d)
   691     .optionGroup("dist", "disc")
   692     .refOption("square", "Nodes are evenly distributed on a unit square", square_d)
   693     .optionGroup("dist", "square")
   694     .refOption("gauss",
   695             "Nodes are located according to a two-dim gauss distribution",
   696             gauss_d)
   697     .optionGroup("dist", "gauss")
   698 //     .mandatoryGroup("dist")
   699     .onlyOneGroup("dist")
   700     .boolOption("eps", "Also generate .eps output (prefix.eps)")
   701     .boolOption("nonodes", "Draw the edges only in the generated .eps")
   702     .boolOption("dir", "Directed digraph is generated (each arcs are replaced by two directed ones)")
   703     .boolOption("2con", "Create a two connected planar digraph")
   704     .optionGroup("alg","2con")
   705     .boolOption("tree", "Create a min. cost spanning tree")
   706     .optionGroup("alg","tree")
   707     .boolOption("tsp", "Create a TSP tour")
   708     .optionGroup("alg","tsp")
   709     .boolOption("tsp2", "Create a TSP tour (tree based)")
   710     .optionGroup("alg","tsp2")
   711     .boolOption("dela", "Delaunay triangulation digraph")
   712     .optionGroup("alg","dela")
   713     .onlyOneGroup("alg")
   714     .boolOption("rand", "Use time seed for random number generator")
   715     .optionGroup("rand", "rand")
   716     .intOption("seed", "Random seed", -1)
   717     .optionGroup("rand", "seed")
   718     .onlyOneGroup("rand")
   719     .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
   720     .run();
   721 
   722   if (ap["rand"]) {
   723     int seed = time(0);
   724     std::cout << "Random number seed: " << seed << std::endl;
   725     rnd = Random(seed);
   726   }
   727   if (ap.given("seed")) {
   728     int seed = ap["seed"];
   729     std::cout << "Random number seed: " << seed << std::endl;
   730     rnd = Random(seed);
   731   }
   732 
   733   std::string prefix;
   734   switch(ap.files().size())
   735     {
   736     case 0:
   737       prefix="lgf-gen-out";
   738       break;
   739     case 1:
   740       prefix=ap.files()[0];
   741       break;
   742     default:
   743       std::cerr << "\nAt most one prefix can be given\n\n";
   744       exit(1);
   745     }
   746 
   747   double sum_sizes=0;
   748   std::vector<double> sizes;
   749   std::vector<double> cum_sizes;
   750   for(int s=0;s<num_of_cities;s++)
   751     {
   752       //         sum_sizes+=rnd.exponential();
   753       double d=rnd();
   754       sum_sizes+=d;
   755       sizes.push_back(d);
   756       cum_sizes.push_back(sum_sizes);
   757     }
   758   int i=0;
   759   for(int s=0;s<num_of_cities;s++)
   760     {
   761       Point center=(num_of_cities==1?Point(0,0):rnd.disc());
   762       if(gauss_d)
   763         for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
   764           Node n=g.addNode();
   765           nodes.push_back(n);
   766           coords[n]=center+rnd.gauss2()*area*
   767             std::sqrt(sizes[s]/sum_sizes);
   768         }
   769       else if(square_d)
   770         for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
   771           Node n=g.addNode();
   772           nodes.push_back(n);
   773           coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
   774             std::sqrt(sizes[s]/sum_sizes);
   775         }
   776       else if(disc_d || true)
   777         for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
   778           Node n=g.addNode();
   779           nodes.push_back(n);
   780           coords[n]=center+rnd.disc()*area*
   781             std::sqrt(sizes[s]/sum_sizes);
   782         }
   783     }
   784 
   785 //   for (ListGraph::NodeIt n(g); n != INVALID; ++n) {
   786 //     std::cerr << coords[n] << std::endl;
   787 //   }
   788 
   789   if(ap["tsp"]) {
   790     tsp();
   791     std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
   792   }
   793   if(ap["tsp2"]) {
   794     tsp2();
   795     std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
   796   }
   797   else if(ap["2con"]) {
   798     std::cout << "Make triangles\n";
   799     //   triangle();
   800     sparseTriangle(ap["g"]);
   801     std::cout << "Make it sparser\n";
   802     sparse2(ap["g"]);
   803   }
   804   else if(ap["tree"]) {
   805     minTree();
   806   }
   807   else if(ap["dela"]) {
   808     delaunay();
   809   }
   810 
   811 
   812   std::cout << "Number of nodes    : " << countNodes(g) << std::endl;
   813   std::cout << "Number of arcs    : " << countEdges(g) << std::endl;
   814   double tlen=0;
   815   for(EdgeIt e(g);e!=INVALID;++e)
   816     tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
   817   std::cout << "Total arc length  : " << tlen << std::endl;
   818 
   819   if(ap["eps"])
   820     graphToEps(g,prefix+".eps").scaleToA4().
   821       scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false).
   822       coords(coords).hideNodes(ap.given("nonodes")).run();
   823 
   824   if(ap["dir"])
   825     DigraphWriter<ListGraph>(g,prefix+".lgf").
   826       nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
   827       nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
   828       run();
   829   else GraphWriter<ListGraph>(g,prefix+".lgf").
   830          nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
   831          nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
   832          run();
   833 }
   834