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