/* -*- 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.

 *

 */

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

}