tools/lgf-gen.cc
author Akos Ladanyi <ladanyi@tmit.bme.hu>
Mon, 27 Apr 2009 18:03:18 +0100
changeset 635 89705c452130
parent 616 f2d6d3446adf
child 654 a312f84d86c6
permissions -rw-r--r--
Add CPLEX_ROOT_DIR variable to FindCPLEX (#277)
     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);
   484       int k=sur.run(a,b,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> > sur(g,cegy);
   515         int k=sur.run(pi->a,pi->b,2);
   516         if(k<2 || sur.totalLength()>d)
   517           {
   518             ne=g.addEdge(pi->a,pi->b);
   519             arcs.push_back(ne);
   520             cnt++;
   521           }
   522       }
   523     }
   524 }
   525 
   526 template <typename Graph, typename CoordMap>
   527 class LengthSquareMap {
   528 public:
   529   typedef typename Graph::Edge Key;
   530   typedef typename CoordMap::Value::Value Value;
   531 
   532   LengthSquareMap(const Graph& graph, const CoordMap& coords)
   533     : _graph(graph), _coords(coords) {}
   534 
   535   Value operator[](const Key& key) const {
   536     return (_coords[_graph.v(key)] -
   537             _coords[_graph.u(key)]).normSquare();
   538   }
   539 
   540 private:
   541 
   542   const Graph& _graph;
   543   const CoordMap& _coords;
   544 };
   545 
   546 void minTree() {
   547   std::vector<Parc> pedges;
   548   Timer T;
   549   std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
   550   delaunay();
   551   std::cout << T.realTime() << "s: Calculating spanning tree...\n";
   552   LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords);
   553   ListGraph::EdgeMap<bool> tree(g);
   554   kruskal(g, ls, tree);
   555   std::cout << T.realTime() << "s: Removing non tree arcs...\n";
   556   std::vector<Edge> remove;
   557   for (EdgeIt e(g); e != INVALID; ++e) {
   558     if (!tree[e]) remove.push_back(e);
   559   }
   560   for(int i = 0; i < int(remove.size()); ++i) {
   561     g.erase(remove[i]);
   562   }
   563   std::cout << T.realTime() << "s: Done\n";
   564 }
   565 
   566 void tsp2()
   567 {
   568   std::cout << "Find a tree..." << std::endl;
   569 
   570   minTree();
   571 
   572   std::cout << "Total arc length (tree) : " << totalLen() << std::endl;
   573 
   574   std::cout << "Make it Euler..." << std::endl;
   575 
   576   {
   577     std::vector<Node> leafs;
   578     for(NodeIt n(g);n!=INVALID;++n)
   579       if(countIncEdges(g,n)%2==1) leafs.push_back(n);
   580 
   581 //    for(unsigned int i=0;i<leafs.size();i+=2)
   582 //       g.addArc(leafs[i],leafs[i+1]);
   583 
   584     std::vector<Parc> pedges;
   585     for(unsigned int i=0;i<leafs.size()-1;i++)
   586       for(unsigned int j=i+1;j<leafs.size();j++)
   587         {
   588           Node n=leafs[i];
   589           Node m=leafs[j];
   590           Parc p;
   591           p.a=n;
   592           p.b=m;
   593           p.len=(coords[m]-coords[n]).normSquare();
   594           pedges.push_back(p);
   595         }
   596     std::sort(pedges.begin(),pedges.end(),pedgeLess);
   597     for(unsigned int i=0;i<pedges.size();i++)
   598       if(countIncEdges(g,pedges[i].a)%2 &&
   599          countIncEdges(g,pedges[i].b)%2)
   600         g.addEdge(pedges[i].a,pedges[i].b);
   601   }
   602 
   603   for(NodeIt n(g);n!=INVALID;++n)
   604     if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
   605       std::cout << "GEBASZ!!!" << std::endl;
   606 
   607   for(EdgeIt e(g);e!=INVALID;++e)
   608     if(g.u(e)==g.v(e))
   609       std::cout << "LOOP GEBASZ!!!" << std::endl;
   610 
   611   std::cout << "Number of arcs : " << countEdges(g) << std::endl;
   612 
   613   std::cout << "Total arc length (euler) : " << totalLen() << std::endl;
   614 
   615   ListGraph::EdgeMap<Arc> enext(g);
   616   {
   617     EulerIt<ListGraph> e(g);
   618     Arc eo=e;
   619     Arc ef=e;
   620 //     std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
   621     for(++e;e!=INVALID;++e)
   622       {
   623 //         std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
   624         enext[eo]=e;
   625         eo=e;
   626       }
   627     enext[eo]=ef;
   628   }
   629 
   630   std::cout << "Creating a tour from that..." << std::endl;
   631 
   632   int nnum = countNodes(g);
   633   int ednum = countEdges(g);
   634 
   635   for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p])
   636     {
   637 //       std::cout << "Checking arc " << g.id(p) << std::endl;
   638       Arc e=enext[p];
   639       Arc f=enext[e];
   640       Node n2=g.source(f);
   641       Node n1=g.oppositeNode(n2,e);
   642       Node n3=g.oppositeNode(n2,f);
   643       if(countIncEdges(g,n2)>2)
   644         {
   645 //           std::cout << "Remove an Arc" << std::endl;
   646           Arc ff=enext[f];
   647           g.erase(e);
   648           g.erase(f);
   649           if(n1!=n3)
   650             {
   651               Arc ne=g.direct(g.addEdge(n1,n3),n1);
   652               enext[p]=ne;
   653               enext[ne]=ff;
   654               ednum--;
   655             }
   656           else {
   657             enext[p]=ff;
   658             ednum-=2;
   659           }
   660         }
   661     }
   662 
   663   std::cout << "Total arc length (tour) : " << totalLen() << std::endl;
   664 
   665   std::cout << "2-opt the tour..." << std::endl;
   666 
   667   tsp_improve();
   668 
   669   std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl;
   670 }
   671 
   672 
   673 int main(int argc,const char **argv)
   674 {
   675   ArgParser ap(argc,argv);
   676 
   677 //   bool eps;
   678   bool disc_d, square_d, gauss_d;
   679 //   bool tsp_a,two_a,tree_a;
   680   int num_of_cities=1;
   681   double area=1;
   682   N=100;
   683 //   girth=10;
   684   std::string ndist("disc");
   685   ap.refOption("n", "Number of nodes (default is 100)", N)
   686     .intOption("g", "Girth parameter (default is 10)", 10)
   687     .refOption("cities", "Number of cities (default is 1)", num_of_cities)
   688     .refOption("area", "Full relative area of the cities (default is 1)", area)
   689     .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d)
   690     .optionGroup("dist", "disc")
   691     .refOption("square", "Nodes are evenly distributed on a unit square", square_d)
   692     .optionGroup("dist", "square")
   693     .refOption("gauss",
   694             "Nodes are located according to a two-dim gauss distribution",
   695             gauss_d)
   696     .optionGroup("dist", "gauss")
   697 //     .mandatoryGroup("dist")
   698     .onlyOneGroup("dist")
   699     .boolOption("eps", "Also generate .eps output (prefix.eps)")
   700     .boolOption("nonodes", "Draw the edges only in the generated .eps")
   701     .boolOption("dir", "Directed digraph is generated (each arcs are replaced by two directed ones)")
   702     .boolOption("2con", "Create a two connected planar digraph")
   703     .optionGroup("alg","2con")
   704     .boolOption("tree", "Create a min. cost spanning tree")
   705     .optionGroup("alg","tree")
   706     .boolOption("tsp", "Create a TSP tour")
   707     .optionGroup("alg","tsp")
   708     .boolOption("tsp2", "Create a TSP tour (tree based)")
   709     .optionGroup("alg","tsp2")
   710     .boolOption("dela", "Delaunay triangulation digraph")
   711     .optionGroup("alg","dela")
   712     .onlyOneGroup("alg")
   713     .boolOption("rand", "Use time seed for random number generator")
   714     .optionGroup("rand", "rand")
   715     .intOption("seed", "Random seed", -1)
   716     .optionGroup("rand", "seed")
   717     .onlyOneGroup("rand")
   718     .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
   719     .run();
   720 
   721   if (ap["rand"]) {
   722     int seed = int(time(0));
   723     std::cout << "Random number seed: " << seed << std::endl;
   724     rnd = Random(seed);
   725   }
   726   if (ap.given("seed")) {
   727     int seed = ap["seed"];
   728     std::cout << "Random number seed: " << seed << std::endl;
   729     rnd = Random(seed);
   730   }
   731 
   732   std::string prefix;
   733   switch(ap.files().size())
   734     {
   735     case 0:
   736       prefix="lgf-gen-out";
   737       break;
   738     case 1:
   739       prefix=ap.files()[0];
   740       break;
   741     default:
   742       std::cerr << "\nAt most one prefix can be given\n\n";
   743       exit(1);
   744     }
   745 
   746   double sum_sizes=0;
   747   std::vector<double> sizes;
   748   std::vector<double> cum_sizes;
   749   for(int s=0;s<num_of_cities;s++)
   750     {
   751       //         sum_sizes+=rnd.exponential();
   752       double d=rnd();
   753       sum_sizes+=d;
   754       sizes.push_back(d);
   755       cum_sizes.push_back(sum_sizes);
   756     }
   757   int i=0;
   758   for(int s=0;s<num_of_cities;s++)
   759     {
   760       Point center=(num_of_cities==1?Point(0,0):rnd.disc());
   761       if(gauss_d)
   762         for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
   763           Node n=g.addNode();
   764           nodes.push_back(n);
   765           coords[n]=center+rnd.gauss2()*area*
   766             std::sqrt(sizes[s]/sum_sizes);
   767         }
   768       else if(square_d)
   769         for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
   770           Node n=g.addNode();
   771           nodes.push_back(n);
   772           coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
   773             std::sqrt(sizes[s]/sum_sizes);
   774         }
   775       else if(disc_d || true)
   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.disc()*area*
   780             std::sqrt(sizes[s]/sum_sizes);
   781         }
   782     }
   783 
   784 //   for (ListGraph::NodeIt n(g); n != INVALID; ++n) {
   785 //     std::cerr << coords[n] << std::endl;
   786 //   }
   787 
   788   if(ap["tsp"]) {
   789     tsp();
   790     std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
   791   }
   792   if(ap["tsp2"]) {
   793     tsp2();
   794     std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
   795   }
   796   else if(ap["2con"]) {
   797     std::cout << "Make triangles\n";
   798     //   triangle();
   799     sparseTriangle(ap["g"]);
   800     std::cout << "Make it sparser\n";
   801     sparse2(ap["g"]);
   802   }
   803   else if(ap["tree"]) {
   804     minTree();
   805   }
   806   else if(ap["dela"]) {
   807     delaunay();
   808   }
   809 
   810 
   811   std::cout << "Number of nodes    : " << countNodes(g) << std::endl;
   812   std::cout << "Number of arcs    : " << countEdges(g) << std::endl;
   813   double tlen=0;
   814   for(EdgeIt e(g);e!=INVALID;++e)
   815     tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
   816   std::cout << "Total arc length  : " << tlen << std::endl;
   817 
   818   if(ap["eps"])
   819     graphToEps(g,prefix+".eps").scaleToA4().
   820       scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false).
   821       coords(coords).hideNodes(ap.given("nonodes")).run();
   822 
   823   if(ap["dir"])
   824     DigraphWriter<ListGraph>(g,prefix+".lgf").
   825       nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
   826       nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
   827       run();
   828   else GraphWriter<ListGraph>(g,prefix+".lgf").
   829          nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
   830          nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
   831          run();
   832 }
   833