#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 <cmath>

 #include <algorithm>

 #include <lemon/unionfind.h>

 using namespace lemon;

 typedef dim2::Point<double> Point;

 UGRAPH_TYPEDEFS(ListUGraph);

 int N;

 int girth;

 ListUGraph g;

 std::vector<Node> nodes;

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

 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;

 void triangle()

 {

   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) {

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

 	edges.push_back(ne);

 	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++;

 	  }

       }

     }

 }

 void minTree() {

   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);

   ListUGraph::NodeMap<int> comp(g);

   UnionFind<ListUGraph::NodeMap<int> > uf(comp);

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

     uf.insert(it);

   int en=0;

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

     {

       if ( uf.join(pi->a,pi->b) ) {

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

 	en++;

 	if(en>=N-1) return;

       }

     }

 }

 int main(int argc,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.option("n", "Number of nodes (default is 100)", N)

     .option("g", "Girth parameter (default is 10)", girth)

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

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

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

     .optionGroup("dist", "disc")

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

     .optionGroup("dist", "square")

     .option("gauss",

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

 	    gauss_d)

     .optionGroup("dist", "gauss")

 //     .mandatoryGroup("dist")

     .onlyOneGroup("dist")

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

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

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

     .option("tree", "Create a min. cost spanning tree",tree_a)

     .optionGroup("alg","tree")

     .option("tsp", "Create a TSP tour",tsp_a)

     .optionGroup("alg","tsp")

     .onlyOneGroup("alg")

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

     .run();

   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);

 	}

     }

   if(tsp_a) {

     tsp();

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

   }

   else if(two_a) {

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

     //   triangle();

     sparseTriangle(girth);

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

     sparse2(girth);

   }

   else if(tree_a) {

     minTree();

   }

   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(eps)

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

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

       coords(coords).run();

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

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

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

     run();

 }