tools/lgf-gen.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Tue, 24 Mar 2009 00:18:25 +0100
changeset 604 8c3112a66878
parent 523 d9e43511d11c
child 570 ab6da8cf5ab2
permissions -rw-r--r--
Use XTI implementation instead of ATI in NetworkSimplex (#234)

XTI (eXtended Threaded Index) is an imporved version of the widely
known ATI (Augmented Threaded Index) method for storing and updating
the spanning tree structure in Network Simplex algorithms.

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