/* -*- C++ -*-

  *

  * This file is a part of LEMON, a generic C++ optimization library

  *

  * Copyright (C) 2003-2007

  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

  * (Egervary Research Group on Combinatorial Optimization, EGRES).

  *

  * Permission to use, modify and distribute this software is granted

  * provided that this copyright notice appears in all copies. For

  * precise terms see the accompanying LICENSE file.

  *

  * This software is provided "AS IS" with no warranty of any kind,

  * express or implied, and with no claim as to its suitability for any

  * purpose.

  *

  */

 ///\ingroup tools

 ///\file

 ///\brief Special plane graph generator.

 ///

 ///Graph generator application for various types of plane graphs.

 ///

 ///\verbatim

 /// Usage:

 ///   ./tools/lgf-gen [-2con|-tree|-tsp|-tsp2|-dela] [-disc|-square|-gauss]

 ///      [-rand|-seed int] [--help|-h|-help] [-area num] [-cities int] [-dir]

 ///      [-eps] [-g int] [-n int] [prefix]

 /// Where:

 ///   [prefix]

 ///      Prefix of the output files. Default is 'lgf-gen-out'

 ///   --help|-h|-help

 ///      Print a short help message

 ///   -2con

 ///      Create a two connected planar graph

 ///   -area num

 ///      Full relative area of the cities (default is 1)

 ///   -cities int

 ///      Number of cities (default is 1)

 ///   -dela

 ///      Delaunay triangulation graph

 ///   -dir

 ///      Directed graph is generated (each edges are replaced by two directed ones)

 ///   -disc

 ///      Nodes are evenly distributed on a unit disc (default)

 ///   -eps

 ///      Also generate .eps output (prefix.eps)

 ///   -g int

 ///      Girth parameter (default is 10)

 ///   -gauss

 ///      Nodes are located according to a two-dim gauss distribution

 ///   -n int

 ///      Number of nodes (default is 100)

 ///   -rand

 ///      Use time seed for random number generator

 ///   -seed int

 ///      Random seed

 ///   -square

 ///      Nodes are evenly distributed on a unit square

 ///   -tree

 ///      Create a min. cost spanning tree

 ///   -tsp

 ///      Create a TSP tour

 ///   -tsp2

 ///      Create a TSP tour (tree based)

 ///\endverbatim

 /// \image html plane_tree.png

 /// \image latex plane_tree.eps "Eucledian spanning tree" width=\textwidth

 ///

 #include <lemon/list_graph.h>

 #include <lemon/graph_utils.h>

 #include <lemon/random.h>

 #include <lemon/dim2.h>

 #include <lemon/bfs.h>

 #include <lemon/counter.h>

 #include <lemon/suurballe.h>

 #include <lemon/graph_to_eps.h>

 #include <lemon/graph_writer.h>

 #include <lemon/arg_parser.h>

 #include <lemon/euler.h>

 #include <cmath>

 #include <algorithm>

 #include <lemon/kruskal.h>

 #include <lemon/time_measure.h>

 using namespace lemon;

 typedef dim2::Point<double> Point;

 UGRAPH_TYPEDEFS(ListUGraph);

 bool progress=true;

 int N;

 // int girth;

 ListUGraph g;

 std::vector<Node> nodes;

 ListUGraph::NodeMap<Point> coords(g);

 double totalLen(){

   double tlen=0;

   for(UEdgeIt e(g);e!=INVALID;++e)

     tlen+=sqrt((coords[g.source(e)]-coords[g.target(e)]).normSquare());

   return tlen;

 }

 int tsp_impr_num=0;

 const double EPSILON=1e-8; 

 bool tsp_improve(Node u, Node v)

 {

   double luv=std::sqrt((coords[v]-coords[u]).normSquare());

   Node u2=u;

   Node v2=v;

   do {

     Node n;

     for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e);

     u2=v2;

     v2=n;

     if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON>

        std::sqrt((coords[u]-coords[u2]).normSquare())+

        std::sqrt((coords[v]-coords[v2]).normSquare()))

       {

  	g.erase(findUEdge(g,u,v));

  	g.erase(findUEdge(g,u2,v2));

 	g.addEdge(u2,u);

 	g.addEdge(v,v2);

 	tsp_impr_num++;

 	return true;

       }

   } while(v2!=u);

   return false;

 }

 bool tsp_improve(Node u)

 {

   for(IncEdgeIt e(g,u);e!=INVALID;++e)

     if(tsp_improve(u,g.runningNode(e))) return true;

   return false;

 }

 void tsp_improve()

 {

   bool b;

   do {

     b=false;

     for(NodeIt n(g);n!=INVALID;++n)

       if(tsp_improve(n)) b=true;

   } while(b);

 }

 void tsp()

 {

   for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);

   tsp_improve();

 }

 class Line

 {

 public:

   Point a;

   Point b;

   Line(Point _a,Point _b) :a(_a),b(_b) {}

   Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}

   Line(const Edge &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}

   Line(const UEdge &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}

 };

 inline std::ostream& operator<<(std::ostream &os, const Line &l)

 {

   os << l.a << "->" << l.b;

   return os;

 }

 bool cross(Line a, Line b) 

 {

   Point ao=rot90(a.b-a.a);

   Point bo=rot90(b.b-b.a);

   return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 &&

     (bo*(a.a-b.a))*(bo*(a.b-b.a))<0;

 }

 struct Pedge

 {

   Node a;

   Node b;

   double len;

 };

 bool pedgeLess(Pedge a,Pedge b)

 {

   return a.len<b.len;

 }

 std::vector<UEdge> edges;

 namespace _delaunay_bits {

   struct Part {

     int prev, curr, next;

     Part(int p, int c, int n) : prev(p), curr(c), next(n) {} 

   };

   inline std::ostream& operator<<(std::ostream& os, const Part& part) {

     os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';

     return os;

   }

   inline double circle_point(const Point& p, const Point& q, const Point& r) {

     double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);

     if (a == 0) return std::numeric_limits<double>::quiet_NaN();

     double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +

       (q.x * q.x + q.y * q.y) * (r.y - p.y) +

       (r.x * r.x + r.y * r.y) * (p.y - q.y);

     double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +

       (q.x * q.x + q.y * q.y) * (r.x - p.x) +

       (r.x * r.x + r.y * r.y) * (p.x - q.x);

     double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +

       (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +

       (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);

     return d / (2 * a) + sqrt((d * d + e * e) / (4 * a * a) + f / a);

   }

   inline bool circle_form(const Point& p, const Point& q, const Point& r) {

     return rot90(q - p) * (r - q) < 0.0;

   }

   inline double intersection(const Point& p, const Point& q, double sx) {

     const double epsilon = 1e-8;

     if (p.x == q.x) return (p.y + q.y) / 2.0;

     if (sx < p.x + epsilon) return p.y;

     if (sx < q.x + epsilon) return q.y;

     double a = q.x - p.x;

     double b = (q.x - sx) * p.y - (p.x - sx) * q.y;    

     double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();

     return (b - sqrt(d)) / a;

   }

   struct YLess {

     YLess(const std::vector<Point>& points, double& sweep) 

       : _points(points), _sweep(sweep) {}

     bool operator()(const Part& l, const Part& r) const {

       const double epsilon = 1e-8;

       //      std::cerr << l << " vs " << r << std::endl;

       double lbx = l.prev != -1 ?

 	intersection(_points[l.prev], _points[l.curr], _sweep) :

 	- std::numeric_limits<double>::infinity();

       double rbx = r.prev != -1 ?

 	intersection(_points[r.prev], _points[r.curr], _sweep) :

 	- std::numeric_limits<double>::infinity();

       double lex = l.next != -1 ?

 	intersection(_points[l.curr], _points[l.next], _sweep) :

 	std::numeric_limits<double>::infinity();

       double rex = r.next != -1 ?

 	intersection(_points[r.curr], _points[r.next], _sweep) :

 	std::numeric_limits<double>::infinity();

       if (lbx > lex) std::swap(lbx, lex);

       if (rbx > rex) std::swap(rbx, rex);

       if (lex < epsilon + rex && lbx + epsilon < rex) return true;

       if (rex < epsilon + lex && rbx + epsilon < lex) return false;

       return lex < rex;

     }

     const std::vector<Point>& _points;

     double& _sweep;

   };

   struct BeachIt;

   typedef std::multimap<double, BeachIt> SpikeHeap;

   typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;

   struct BeachIt {

     Beach::iterator it;

     BeachIt(Beach::iterator iter) : it(iter) {}

   };

 }

 inline void delaunay() {

   Counter cnt("Number of edges added: ");

   using namespace _delaunay_bits;

   typedef _delaunay_bits::Part Part;

   typedef std::vector<std::pair<double, int> > SiteHeap;

   std::vector<Point> points;

   std::vector<Node> nodes;

   for (NodeIt it(g); it != INVALID; ++it) {

     nodes.push_back(it);

     points.push_back(coords[it]);

   }

   SiteHeap siteheap(points.size());

   double sweep;

   for (int i = 0; i < int(siteheap.size()); ++i) {

     siteheap[i] = std::make_pair(points[i].x, i);

   }

   std::sort(siteheap.begin(), siteheap.end());

   sweep = siteheap.front().first;

   YLess yless(points, sweep);

   Beach beach(yless);

   SpikeHeap spikeheap;

   std::set<std::pair<int, int> > edges;

   int siteindex = 0;

   {

     SiteHeap front;

     while (siteindex < int(siteheap.size()) &&

 	   siteheap[0].first == siteheap[siteindex].first) {

       front.push_back(std::make_pair(points[siteheap[siteindex].second].y,

 				     siteheap[siteindex].second));

       ++siteindex;

     }

     std::sort(front.begin(), front.end());

     for (int i = 0; i < int(front.size()); ++i) {

       int prev = (i == 0 ? -1 : front[i - 1].second);

       int curr = front[i].second;

       int next = (i + 1 == int(front.size()) ? -1 : front[i + 1].second);

       beach.insert(std::make_pair(Part(prev, curr, next), 

 				  spikeheap.end()));      

     }

   }

   while (siteindex < int(points.size()) || !spikeheap.empty()) {

     SpikeHeap::iterator spit = spikeheap.begin();

     if (siteindex < int(points.size()) && 

 	(spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {

       int site = siteheap[siteindex].second;

       sweep = siteheap[siteindex].first;

       Beach::iterator bit = beach.upper_bound(Part(site, site, site));

       if (bit->second != spikeheap.end()) {

 	spikeheap.erase(bit->second);	

       }

       int prev = bit->first.prev;

       int curr = bit->first.curr;

       int next = bit->first.next;

       beach.erase(bit);

       SpikeHeap::iterator pit = spikeheap.end();

       if (prev != -1 && 

 	  circle_form(points[prev], points[curr], points[site])) {

 	double x = circle_point(points[prev], points[curr], points[site]);

 	pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));

 	pit->second.it =

 	  beach.insert(std::make_pair(Part(prev, curr, site), pit));

       } else {

 	beach.insert(std::make_pair(Part(prev, curr, site), pit));

       }

       beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));

       SpikeHeap::iterator nit = spikeheap.end();

       if (next != -1 && 

 	  circle_form(points[site], points[curr],points[next])) {

 	double x = circle_point(points[site], points[curr], points[next]);

 	nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));

 	nit->second.it =

 	  beach.insert(std::make_pair(Part(site, curr, next), nit));

       } else {

 	beach.insert(std::make_pair(Part(site, curr, next), nit));

       }

       ++siteindex;

     } else {

       sweep = spit->first;      

       Beach::iterator bit = spit->second.it;

       int prev = bit->first.prev;

       int curr = bit->first.curr;

       int next = bit->first.next;

       {

 	std::pair<int, int> edge;

 	edge = prev < curr ? 

 	  std::make_pair(prev, curr) : std::make_pair(curr, prev);

 	if (edges.find(edge) == edges.end()) {

 	  edges.insert(edge);

 	  g.addEdge(nodes[prev], nodes[curr]);

 	  ++cnt;

 	}

 	edge = curr < next ? 

 	  std::make_pair(curr, next) : std::make_pair(next, curr);

 	if (edges.find(edge) == edges.end()) {

 	  edges.insert(edge);

 	  g.addEdge(nodes[curr], nodes[next]);

 	  ++cnt;

 	}

       }

       Beach::iterator pbit = bit; --pbit;

       int ppv = pbit->first.prev;

       Beach::iterator nbit = bit; ++nbit;

       int nnt = nbit->first.next;

       if (bit->second != spikeheap.end()) spikeheap.erase(bit->second);

       if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second);

       if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second);

       beach.erase(nbit);

       beach.erase(bit);

       beach.erase(pbit);

       SpikeHeap::iterator pit = spikeheap.end();

       if (ppv != -1 && ppv != next && 

 	  circle_form(points[ppv], points[prev], points[next])) {

 	double x = circle_point(points[ppv], points[prev], points[next]);

 	if (x < sweep) x = sweep;

 	pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));

 	pit->second.it =

 	  beach.insert(std::make_pair(Part(ppv, prev, next), pit));

       } else {

 	beach.insert(std::make_pair(Part(ppv, prev, next), pit));

       }

       SpikeHeap::iterator nit = spikeheap.end();

       if (nnt != -1 && prev != nnt &&

 	  circle_form(points[prev], points[next], points[nnt])) {

 	double x = circle_point(points[prev], points[next], points[nnt]);

 	if (x < sweep) x = sweep;

 	nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));

 	nit->second.it =

 	  beach.insert(std::make_pair(Part(prev, next, nnt), nit));

       } else {

 	beach.insert(std::make_pair(Part(prev, next, nnt), nit));

       }

     }

   }

   for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {

     int curr = it->first.curr;

     int next = it->first.next;

     if (next == -1) continue;

     std::pair<int, int> edge;

     edge = curr < next ? 

       std::make_pair(curr, next) : std::make_pair(next, curr);

     if (edges.find(edge) == edges.end()) {

       edges.insert(edge);

       g.addEdge(nodes[curr], nodes[next]);

       ++cnt;

     }

   }

 }

 void sparse(int d) 

 {

   Counter cnt("Number of edges removed: ");

   Bfs<ListUGraph> bfs(g);

   for(std::vector<UEdge>::reverse_iterator ei=edges.rbegin();

       ei!=edges.rend();++ei)

     {

       Node a=g.source(*ei);

       Node b=g.target(*ei);

       g.erase(*ei);

       bfs.run(a,b);

       if(bfs.predEdge(b)==INVALID || bfs.dist(b)>d)

 	g.addEdge(a,b);

       else cnt++;

     }

 }

 void sparse2(int d) 

 {

   Counter cnt("Number of edges removed: ");

   for(std::vector<UEdge>::reverse_iterator ei=edges.rbegin();

       ei!=edges.rend();++ei)

     {

       Node a=g.source(*ei);

       Node b=g.target(*ei);

       g.erase(*ei);

       ConstMap<Edge,int> cegy(1);

       Suurballe<ListUGraph,ConstMap<Edge,int> > sur(g,cegy,a,b);

       int k=sur.run(2);

       if(k<2 || sur.totalLength()>d)

 	g.addEdge(a,b);

       else cnt++;

 //       else std::cout << "Remove edge " << g.id(a) << "-" << g.id(b) << '\n';

     }

 }

 void sparseTriangle(int d)

 {

   Counter cnt("Number of edges added: ");

   std::vector<Pedge> pedges;

   for(NodeIt n(g);n!=INVALID;++n) 

     for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)

       {

 	Pedge p;

 	p.a=n;

 	p.b=m;

 	p.len=(coords[m]-coords[n]).normSquare();

 	pedges.push_back(p);

       }

   std::sort(pedges.begin(),pedges.end(),pedgeLess);

   for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)

     {

       Line li(pi->a,pi->b);

       UEdgeIt e(g);

       for(;e!=INVALID && !cross(e,li);++e) ;

       UEdge ne;

       if(e==INVALID) {

 	ConstMap<Edge,int> cegy(1);

 	Suurballe<ListUGraph,ConstMap<Edge,int> >

 	  sur(g,cegy,pi->a,pi->b);

 	int k=sur.run(2);

 	if(k<2 || sur.totalLength()>d)

 	  {

 	    ne=g.addEdge(pi->a,pi->b);

 	    edges.push_back(ne);

 	    cnt++;

 	  }

       }

     }

 }

 template <typename UGraph, typename CoordMap>

 class LengthSquareMap {

 public:

   typedef typename UGraph::UEdge Key;

   typedef typename CoordMap::Value::Value Value;

   LengthSquareMap(const UGraph& ugraph, const CoordMap& coords)

     : _ugraph(ugraph), _coords(coords) {}

   Value operator[](const Key& key) const {

     return (_coords[_ugraph.target(key)] -

 	    _coords[_ugraph.source(key)]).normSquare();

   }

 private:

   const UGraph& _ugraph;

   const CoordMap& _coords;

 };

 void minTree() {

   std::vector<Pedge> pedges;

   Timer T;

   std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";

   delaunay();

   std::cout << T.realTime() << "s: Calculating spanning tree...\n";

   LengthSquareMap<ListUGraph, ListUGraph::NodeMap<Point> > ls(g, coords);

   ListUGraph::UEdgeMap<bool> tree(g);

   kruskal(g, ls, tree);

   std::cout << T.realTime() << "s: Removing non tree edges...\n";

   std::vector<UEdge> remove;

   for (UEdgeIt e(g); e != INVALID; ++e) {

     if (!tree[e]) remove.push_back(e);

   }

   for(int i = 0; i < int(remove.size()); ++i) {

     g.erase(remove[i]);

   }

   std::cout << T.realTime() << "s: Done\n";

 }

 void tsp2() 

 {

   std::cout << "Find a tree..." << std::endl;

   minTree();

   std::cout << "Total edge length (tree) : " << totalLen() << std::endl;

   std::cout << "Make it Euler..." << std::endl;

   {

     std::vector<Node> leafs;

     for(NodeIt n(g);n!=INVALID;++n)

       if(countIncEdges(g,n)%2==1) leafs.push_back(n);

 //    for(unsigned int i=0;i<leafs.size();i+=2)

 //       g.addEdge(leafs[i],leafs[i+1]);

     std::vector<Pedge> pedges;

     for(unsigned int i=0;i<leafs.size()-1;i++)

       for(unsigned int j=i+1;j<leafs.size();j++)

 	{

 	  Node n=leafs[i];

 	  Node m=leafs[j];

 	  Pedge p;

 	  p.a=n;

 	  p.b=m;

 	  p.len=(coords[m]-coords[n]).normSquare();

 	  pedges.push_back(p);

 	}

     std::sort(pedges.begin(),pedges.end(),pedgeLess);

     for(unsigned int i=0;i<pedges.size();i++)

       if(countIncEdges(g,pedges[i].a)%2 &&

 	 countIncEdges(g,pedges[i].b)%2)

 	g.addEdge(pedges[i].a,pedges[i].b);

   }

   for(NodeIt n(g);n!=INVALID;++n)

     if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )

       std::cout << "GEBASZ!!!" << std::endl;

   for(UEdgeIt e(g);e!=INVALID;++e)

     if(g.source(e)==g.target(e))

       std::cout << "LOOP GEBASZ!!!" << std::endl;

   std::cout << "Number of edges : " << countUEdges(g) << std::endl;

   std::cout << "Total edge length (euler) : " << totalLen() << std::endl;

   ListUGraph::UEdgeMap<Edge> enext(g);

   {

     UEulerIt<ListUGraph> e(g);

     Edge eo=e;

     Edge ef=e;

 //     std::cout << "Tour edge: " << g.id(UEdge(e)) << std::endl;      

     for(++e;e!=INVALID;++e)

       {

 // 	std::cout << "Tour edge: " << g.id(UEdge(e)) << std::endl;      

 	enext[eo]=e;

 	eo=e;

       }

     enext[eo]=ef;

   }

   std::cout << "Creating a tour from that..." << std::endl;

   int nnum = countNodes(g);

   int ednum = countUEdges(g);

   for(Edge p=enext[UEdgeIt(g)];ednum>nnum;p=enext[p]) 

     {

 //       std::cout << "Checking edge " << g.id(p) << std::endl;      

       Edge e=enext[p];

       Edge f=enext[e];

       Node n2=g.source(f);

       Node n1=g.oppositeNode(n2,e);

       Node n3=g.oppositeNode(n2,f);

       if(countIncEdges(g,n2)>2)

 	{

 // 	  std::cout << "Remove an Edge" << std::endl;

 	  Edge ff=enext[f];

 	  g.erase(e);

 	  g.erase(f);

 	  if(n1!=n3)

 	    {

 	      Edge ne=g.direct(g.addEdge(n1,n3),n1);

 	      enext[p]=ne;

 	      enext[ne]=ff;

 	      ednum--;

 	    }

 	  else {

 	    enext[p]=ff;

 	    ednum-=2;

 	  }

 	}

     }

   std::cout << "Total edge length (tour) : " << totalLen() << std::endl;

   std::cout << "2-opt the tour..." << std::endl;

   tsp_improve();

   std::cout << "Total edge length (2-opt tour) : " << totalLen() << std::endl;

 }

 int main(int argc,const char **argv) 

 {

   ArgParser ap(argc,argv);

 //   bool eps;

   bool disc_d, square_d, gauss_d;

 //   bool tsp_a,two_a,tree_a;

   int num_of_cities=1;

   double area=1;

   N=100;

 //   girth=10;

   std::string ndist("disc");

   ap.refOption("n", "Number of nodes (default is 100)", N)

     .intOption("g", "Girth parameter (default is 10)", 10)

     .refOption("cities", "Number of cities (default is 1)", num_of_cities)

     .refOption("area", "Full relative area of the cities (default is 1)", area)

     .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d)

     .optionGroup("dist", "disc")

     .refOption("square", "Nodes are evenly distributed on a unit square", square_d)

     .optionGroup("dist", "square")

     .refOption("gauss",

 	    "Nodes are located according to a two-dim gauss distribution",

 	    gauss_d)

     .optionGroup("dist", "gauss")

 //     .mandatoryGroup("dist")

     .onlyOneGroup("dist")

     .boolOption("eps", "Also generate .eps output (prefix.eps)")

     .boolOption("dir", "Directed graph is generated (each edges are replaced by two directed ones)")

     .boolOption("2con", "Create a two connected planar graph")

     .optionGroup("alg","2con")

     .boolOption("tree", "Create a min. cost spanning tree")

     .optionGroup("alg","tree")

     .boolOption("tsp", "Create a TSP tour")

     .optionGroup("alg","tsp")

     .boolOption("tsp2", "Create a TSP tour (tree based)")

     .optionGroup("alg","tsp2")

     .boolOption("dela", "Delaunay triangulation graph")

     .optionGroup("alg","dela")

     .onlyOneGroup("alg")

     .boolOption("rand", "Use time seed for random number generator")

     .optionGroup("rand", "rand")

     .intOption("seed", "Random seed", -1)

     .optionGroup("rand", "seed")

     .onlyOneGroup("rand")

     .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")

     .run();

   if (ap["rand"]) {

     int seed = time(0);

     std::cout << "Random number seed: " << seed << std::endl;

     rnd = Random(seed);

   }

   if (ap.given("seed")) {

     int seed = ap["seed"];

     std::cout << "Random number seed: " << seed << std::endl;

     rnd = Random(seed);

   }

   std::string prefix;

   switch(ap.files().size()) 

     {

     case 0:

       prefix="lgf-gen-out";

       break;

     case 1:

       prefix=ap.files()[0];

       break;

     default:

       std::cerr << "\nAt most one prefix can be given\n\n";

       exit(1);

     }

   double sum_sizes=0;

   std::vector<double> sizes;

   std::vector<double> cum_sizes;

   for(int s=0;s<num_of_cities;s++) 

     {

       // 	sum_sizes+=rnd.exponential();

       double d=rnd();

       sum_sizes+=d;

       sizes.push_back(d);

       cum_sizes.push_back(sum_sizes);

     }

   int i=0;

   for(int s=0;s<num_of_cities;s++) 

     {

       Point center=(num_of_cities==1?Point(0,0):rnd.disc());

       if(gauss_d)

 	for(;i<N*(cum_sizes[s]/sum_sizes);i++) {

 	  Node n=g.addNode();

 	  nodes.push_back(n);

 	  coords[n]=center+rnd.gauss2()*area*

 	    std::sqrt(sizes[s]/sum_sizes);

 	}

       else if(square_d)

 	for(;i<N*(cum_sizes[s]/sum_sizes);i++) {

 	  Node n=g.addNode();

 	  nodes.push_back(n);

 	  coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*

 	    std::sqrt(sizes[s]/sum_sizes);

 	}

       else if(disc_d || true)

 	for(;i<N*(cum_sizes[s]/sum_sizes);i++) {

 	  Node n=g.addNode();

 	  nodes.push_back(n);

 	  coords[n]=center+rnd.disc()*area*

 	    std::sqrt(sizes[s]/sum_sizes);

 	}

     }

 //   for (ListUGraph::NodeIt n(g); n != INVALID; ++n) {

 //     std::cerr << coords[n] << std::endl;

 //   }

   if(ap["tsp"]) {

     tsp();

     std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;

   }

   if(ap["tsp2"]) {

     tsp2();

     std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;

   }

   else if(ap["2con"]) {

     std::cout << "Make triangles\n";

     //   triangle();

     sparseTriangle(ap["g"]);

     std::cout << "Make it sparser\n";

     sparse2(ap["g"]);

   }

   else if(ap["tree"]) {

     minTree();

   }

   else if(ap["dela"]) {

     delaunay();

   }

   std::cout << "Number of nodes    : " << countNodes(g) << std::endl;

   std::cout << "Number of edges    : " << countUEdges(g) << std::endl;

   double tlen=0;

   for(UEdgeIt e(g);e!=INVALID;++e)

     tlen+=sqrt((coords[g.source(e)]-coords[g.target(e)]).normSquare());

   std::cout << "Total edge length  : " << tlen << std::endl;

   if(ap["eps"])

     graphToEps(g,prefix+".eps").scaleToA4().

       scale(600).nodeScale(.2).edgeWidthScale(.001).preScale(false).

       coords(coords).run();

   if(ap["dir"])

     GraphWriter<ListUGraph>(prefix+".lgf",g).

       writeNodeMap("coordinates_x",scaleMap(xMap(coords),600)).

       writeNodeMap("coordinates_y",scaleMap(yMap(coords),600)).

       run();

   else UGraphWriter<ListUGraph>(prefix+".lgf",g).

 	 writeNodeMap("coordinates_x",scaleMap(xMap(coords),600)).

 	 writeNodeMap("coordinates_y",scaleMap(yMap(coords),600)).

 	 run();

 }