alpar@2391: /* -*- C++ -*-
alpar@2391:  *
alpar@2391:  * This file is a part of LEMON, a generic C++ optimization library
alpar@2391:  *
alpar@2391:  * Copyright (C) 2003-2007
alpar@2391:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@2391:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@2391:  *
alpar@2391:  * Permission to use, modify and distribute this software is granted
alpar@2391:  * provided that this copyright notice appears in all copies. For
alpar@2391:  * precise terms see the accompanying LICENSE file.
alpar@2391:  *
alpar@2391:  * This software is provided "AS IS" with no warranty of any kind,
alpar@2391:  * express or implied, and with no claim as to its suitability for any
alpar@2391:  * purpose.
alpar@2391:  *
alpar@2391:  */
alpar@2391: 
alpar@2390: #include <lemon/list_graph.h>
alpar@2390: #include <lemon/graph_utils.h>
alpar@2390: #include <lemon/random.h>
alpar@2390: #include <lemon/dim2.h>
alpar@2390: #include <lemon/bfs.h>
alpar@2390: #include <lemon/counter.h>
alpar@2390: #include <lemon/suurballe.h>
alpar@2390: #include <lemon/graph_to_eps.h>
alpar@2390: #include <lemon/graph_writer.h>
alpar@2390: #include <lemon/arg_parser.h>
alpar@2390: #include <cmath>
alpar@2390: #include <algorithm>
alpar@2390: #include <lemon/unionfind.h>
alpar@2402: #include <lemon/time_measure.h>
alpar@2390: 
alpar@2390: using namespace lemon;
alpar@2390: 
alpar@2390: typedef dim2::Point<double> Point;
alpar@2390: 
alpar@2390: UGRAPH_TYPEDEFS(ListUGraph);
alpar@2390: 
alpar@2402: bool progress=true;
alpar@2402: 
alpar@2390: int N;
alpar@2402: // int girth;
alpar@2390: 
alpar@2390: ListUGraph g;
alpar@2390: 
alpar@2390: std::vector<Node> nodes;
alpar@2390: ListUGraph::NodeMap<Point> coords(g);
alpar@2390: 
alpar@2390: int tsp_impr_num=0;
alpar@2390: 
alpar@2390: const double EPSILON=1e-8; 
alpar@2390: bool tsp_improve(Node u, Node v)
alpar@2390: {
alpar@2390:   double luv=std::sqrt((coords[v]-coords[u]).normSquare());
alpar@2390:   Node u2=u;
alpar@2390:   Node v2=v;
alpar@2390:   do {
alpar@2390:     Node n;
alpar@2390:     for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e);
alpar@2390:     u2=v2;
alpar@2390:     v2=n;
alpar@2390:     if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON>
alpar@2390:        std::sqrt((coords[u]-coords[u2]).normSquare())+
alpar@2390:        std::sqrt((coords[v]-coords[v2]).normSquare()))
alpar@2390:       {
alpar@2390:  	g.erase(findUEdge(g,u,v));
alpar@2390:  	g.erase(findUEdge(g,u2,v2));
alpar@2390: 	g.addEdge(u2,u);
alpar@2390: 	g.addEdge(v,v2);
alpar@2390: 	tsp_impr_num++;
alpar@2390: 	return true;
alpar@2390:       }
alpar@2390:   } while(v2!=u);
alpar@2390:   return false;
alpar@2390: }
alpar@2390: 
alpar@2390: bool tsp_improve(Node u)
alpar@2390: {
alpar@2390:   for(IncEdgeIt e(g,u);e!=INVALID;++e)
alpar@2390:     if(tsp_improve(u,g.runningNode(e))) return true;
alpar@2390:   return false;
alpar@2390: }
alpar@2390: 
alpar@2390: void tsp_improve()
alpar@2390: {
alpar@2390:   bool b;
alpar@2390:   do {
alpar@2390:     b=false;
alpar@2390:     for(NodeIt n(g);n!=INVALID;++n)
alpar@2390:       if(tsp_improve(n)) b=true;
alpar@2390:   } while(b);
alpar@2390: }
alpar@2390: 
alpar@2390: void tsp()
alpar@2390: {
alpar@2390:   for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);
alpar@2390:   tsp_improve();
alpar@2390: }
alpar@2390: 
alpar@2390: class Line
alpar@2390: {
alpar@2390: public:
alpar@2390:   Point a;
alpar@2390:   Point b;
alpar@2390:   Line(Point _a,Point _b) :a(_a),b(_b) {}
alpar@2390:   Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}
alpar@2390:   Line(const Edge &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
alpar@2390:   Line(const UEdge &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
alpar@2390: };
alpar@2390:   
alpar@2390: inline std::ostream& operator<<(std::ostream &os, const Line &l)
alpar@2390: {
alpar@2390:   os << l.a << "->" << l.b;
alpar@2390:   return os;
alpar@2390: }
alpar@2390: 
alpar@2390: bool cross(Line a, Line b) 
alpar@2390: {
alpar@2390:   Point ao=rot90(a.b-a.a);
alpar@2390:   Point bo=rot90(b.b-b.a);
alpar@2390:   return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 &&
alpar@2390:     (bo*(a.a-b.a))*(bo*(a.b-b.a))<0;
alpar@2390: }
alpar@2390: 
alpar@2390: struct Pedge
alpar@2390: {
alpar@2390:   Node a;
alpar@2390:   Node b;
alpar@2390:   double len;
alpar@2390: };
alpar@2390: 
alpar@2390: bool pedgeLess(Pedge a,Pedge b)
alpar@2390: {
alpar@2390:   return a.len<b.len;
alpar@2390: }
alpar@2390: 
alpar@2390: std::vector<UEdge> edges;
alpar@2390: 
alpar@2390: void triangle()
alpar@2390: {
alpar@2390:   Counter cnt("Number of edges added: ");
alpar@2390:   std::vector<Pedge> pedges;
alpar@2390:   for(NodeIt n(g);n!=INVALID;++n) 
alpar@2390:     for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
alpar@2390:       {
alpar@2390: 	Pedge p;
alpar@2390: 	p.a=n;
alpar@2390: 	p.b=m;
alpar@2390: 	p.len=(coords[m]-coords[n]).normSquare();
alpar@2390: 	pedges.push_back(p);
alpar@2390:       }
alpar@2390:   std::sort(pedges.begin(),pedges.end(),pedgeLess);
alpar@2390:   for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
alpar@2390:     {
alpar@2390:       Line li(pi->a,pi->b);
alpar@2390:       UEdgeIt e(g);
alpar@2390:       for(;e!=INVALID && !cross(e,li);++e) ;
alpar@2390:       UEdge ne;
alpar@2390:       if(e==INVALID) {
alpar@2390: 	ne=g.addEdge(pi->a,pi->b);
alpar@2390: 	edges.push_back(ne);
alpar@2390: 	cnt++;
alpar@2390:       }
alpar@2390:     }
alpar@2390: }
alpar@2390: 
alpar@2390: void sparse(int d) 
alpar@2390: {
alpar@2390:   Counter cnt("Number of edges removed: ");
alpar@2390:   Bfs<ListUGraph> bfs(g);
alpar@2390:   for(std::vector<UEdge>::reverse_iterator ei=edges.rbegin();
alpar@2390:       ei!=edges.rend();++ei)
alpar@2390:     {
alpar@2390:       Node a=g.source(*ei);
alpar@2390:       Node b=g.target(*ei);
alpar@2390:       g.erase(*ei);
alpar@2390:       bfs.run(a,b);
alpar@2390:       if(bfs.predEdge(b)==INVALID || bfs.dist(b)>d)
alpar@2390: 	g.addEdge(a,b);
alpar@2390:       else cnt++;
alpar@2390:     }
alpar@2390: }
alpar@2390: 
alpar@2390: void sparse2(int d) 
alpar@2390: {
alpar@2390:   Counter cnt("Number of edges removed: ");
alpar@2390:   for(std::vector<UEdge>::reverse_iterator ei=edges.rbegin();
alpar@2390:       ei!=edges.rend();++ei)
alpar@2390:     {
alpar@2390:       Node a=g.source(*ei);
alpar@2390:       Node b=g.target(*ei);
alpar@2390:       g.erase(*ei);
alpar@2390:       ConstMap<Edge,int> cegy(1);
alpar@2390:       Suurballe<ListUGraph,ConstMap<Edge,int> > sur(g,cegy,a,b);
alpar@2390:       int k=sur.run(2);
alpar@2390:       if(k<2 || sur.totalLength()>d)
alpar@2390: 	g.addEdge(a,b);
alpar@2390:       else cnt++;
alpar@2390: //       else std::cout << "Remove edge " << g.id(a) << "-" << g.id(b) << '\n';
alpar@2390:     }
alpar@2390: }
alpar@2390: 
alpar@2390: void sparseTriangle(int d)
alpar@2390: {
alpar@2390:   Counter cnt("Number of edges added: ");
alpar@2390:   std::vector<Pedge> pedges;
alpar@2390:   for(NodeIt n(g);n!=INVALID;++n) 
alpar@2390:     for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
alpar@2390:       {
alpar@2390: 	Pedge p;
alpar@2390: 	p.a=n;
alpar@2390: 	p.b=m;
alpar@2390: 	p.len=(coords[m]-coords[n]).normSquare();
alpar@2390: 	pedges.push_back(p);
alpar@2390:       }
alpar@2390:   std::sort(pedges.begin(),pedges.end(),pedgeLess);
alpar@2390:   for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
alpar@2390:     {
alpar@2390:       Line li(pi->a,pi->b);
alpar@2390:       UEdgeIt e(g);
alpar@2390:       for(;e!=INVALID && !cross(e,li);++e) ;
alpar@2390:       UEdge ne;
alpar@2390:       if(e==INVALID) {
alpar@2390: 	ConstMap<Edge,int> cegy(1);
alpar@2390: 	Suurballe<ListUGraph,ConstMap<Edge,int> >
alpar@2390: 	  sur(g,cegy,pi->a,pi->b);
alpar@2390: 	int k=sur.run(2);
alpar@2390: 	if(k<2 || sur.totalLength()>d)
alpar@2390: 	  {
alpar@2390: 	    ne=g.addEdge(pi->a,pi->b);
alpar@2390: 	    edges.push_back(ne);
alpar@2390: 	    cnt++;
alpar@2390: 	  }
alpar@2390:       }
alpar@2390:     }
alpar@2390: }
alpar@2390: 
alpar@2390: void minTree() {
alpar@2402:   int en=0;
alpar@2402:   int pr=0;
alpar@2390:   std::vector<Pedge> pedges;
alpar@2402:   Timer T;
alpar@2402:   std::cout << T.realTime() << "s: Setting up the edges...\n";
alpar@2390:   for(NodeIt n(g);n!=INVALID;++n) 
alpar@2402:     {
alpar@2402:       for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
alpar@2402: 	{
alpar@2402: 	  Pedge p;
alpar@2402: 	  p.a=n;
alpar@2402: 	  p.b=m;
alpar@2402: 	  p.len=(coords[m]-coords[n]).normSquare();
alpar@2402: 	  pedges.push_back(p);
alpar@2402: 	}
alpar@2402:       en++;
alpar@2402:       if(progress && en>=pr*double(N)/100) 
alpar@2402: 	{
alpar@2402: 	  std::cout << pr << "%  \r" << std::flush;
alpar@2402: 	  pr++;
alpar@2402: 	}
alpar@2402:     }
alpar@2402:   std::cout << T.realTime() << "s: Sorting the edges...\n";
alpar@2390:   std::sort(pedges.begin(),pedges.end(),pedgeLess);
alpar@2402:   std::cout << T.realTime() << "s: Creating the tree...\n";
alpar@2390:   ListUGraph::NodeMap<int> comp(g);
alpar@2390:   UnionFind<ListUGraph::NodeMap<int> > uf(comp);
alpar@2390:   for (NodeIt it(g); it != INVALID; ++it)
alpar@2402:     uf.insert(it);  
alpar@2390:   for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
alpar@2390:     {
alpar@2390:       if ( uf.join(pi->a,pi->b) ) {
alpar@2390: 	g.addEdge(pi->a,pi->b);
alpar@2402: 	if(en>=N-1)
alpar@2402: 	  break;
alpar@2390:       }
alpar@2390:     }
alpar@2402:   std::cout << T.realTime() << "s: Done\n";
alpar@2390: }
alpar@2390: 
alpar@2390: 
alpar@2390: 
deba@2410: int main(int argc,const char **argv) 
alpar@2390: {
alpar@2390:   ArgParser ap(argc,argv);
alpar@2390: 
alpar@2402: //   bool eps;
alpar@2390:   bool disc_d, square_d, gauss_d;
alpar@2402: //   bool tsp_a,two_a,tree_a;
alpar@2390:   int num_of_cities=1;
alpar@2390:   double area=1;
alpar@2390:   N=100;
alpar@2402: //   girth=10;
alpar@2390:   std::string ndist("disc");
alpar@2402:   ap.refOption("n", "Number of nodes (default is 100)", N)
alpar@2402:     .intOption("g", "Girth parameter (default is 10)", 10)
alpar@2402:     .refOption("cities", "Number of cities (default is 1)", num_of_cities)
alpar@2402:     .refOption("area", "Full relative area of the cities (default is 1)", area)
alpar@2402:     .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d)
alpar@2390:     .optionGroup("dist", "disc")
alpar@2402:     .refOption("square", "Nodes are evenly distributed on a unit square", square_d)
alpar@2390:     .optionGroup("dist", "square")
alpar@2402:     .refOption("gauss",
alpar@2390: 	    "Nodes are located according to a two-dim gauss distribution",
alpar@2390: 	    gauss_d)
alpar@2390:     .optionGroup("dist", "gauss")
alpar@2390: //     .mandatoryGroup("dist")
alpar@2390:     .onlyOneGroup("dist")
alpar@2402:     .boolOption("eps", "Also generate .eps output (prefix.eps)")
alpar@2402:     .boolOption("2con", "Create a two connected planar graph")
alpar@2390:     .optionGroup("alg","2con")
alpar@2402:     .boolOption("tree", "Create a min. cost spanning tree")
alpar@2390:     .optionGroup("alg","tree")
alpar@2402:     .boolOption("tsp", "Create a TSP tour")
alpar@2390:     .optionGroup("alg","tsp")
alpar@2390:     .onlyOneGroup("alg")
alpar@2390:     .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
alpar@2390:     .run();
alpar@2390:   
alpar@2390:   std::string prefix;
alpar@2390:   switch(ap.files().size()) 
alpar@2390:     {
alpar@2390:     case 0:
alpar@2390:       prefix="lgf-gen-out";
alpar@2390:       break;
alpar@2390:     case 1:
alpar@2390:       prefix=ap.files()[0];
alpar@2390:       break;
alpar@2390:     default:
alpar@2390:       std::cerr << "\nAt most one prefix can be given\n\n";
alpar@2390:       exit(1);
alpar@2390:     }
alpar@2390:   
alpar@2390:   double sum_sizes=0;
alpar@2390:   std::vector<double> sizes;
alpar@2390:   std::vector<double> cum_sizes;
alpar@2390:   for(int s=0;s<num_of_cities;s++) 
alpar@2390:     {
alpar@2390:       // 	sum_sizes+=rnd.exponential();
alpar@2390:       double d=rnd();
alpar@2390:       sum_sizes+=d;
alpar@2390:       sizes.push_back(d);
alpar@2390:       cum_sizes.push_back(sum_sizes);
alpar@2390:     }
alpar@2390:   int i=0;
alpar@2390:   for(int s=0;s<num_of_cities;s++) 
alpar@2390:     {
alpar@2390:       Point center=(num_of_cities==1?Point(0,0):rnd.disc());
alpar@2390:       if(gauss_d)
alpar@2390: 	for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
alpar@2390: 	  Node n=g.addNode();
alpar@2390: 	  nodes.push_back(n);
alpar@2390: 	  coords[n]=center+rnd.gauss2()*area*
alpar@2390: 	    std::sqrt(sizes[s]/sum_sizes);
alpar@2390: 	}
alpar@2390:       else if(square_d)
alpar@2390: 	for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
alpar@2390: 	  Node n=g.addNode();
alpar@2390: 	  nodes.push_back(n);
alpar@2390: 	  coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
alpar@2390: 	    std::sqrt(sizes[s]/sum_sizes);
alpar@2390: 	}
alpar@2390:       else if(disc_d || true)
alpar@2390: 	for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
alpar@2390: 	  Node n=g.addNode();
alpar@2390: 	  nodes.push_back(n);
alpar@2390: 	  coords[n]=center+rnd.disc()*area*
alpar@2390: 	    std::sqrt(sizes[s]/sum_sizes);
alpar@2390: 	}
alpar@2390:     }
alpar@2390:   
alpar@2402:   if(ap["tsp"]) {
alpar@2390:     tsp();
alpar@2390:     std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
alpar@2390:   }
alpar@2402:   else if(ap["2con"]) {
alpar@2390:     std::cout << "Make triangles\n";
alpar@2390:     //   triangle();
alpar@2402:     sparseTriangle(ap["g"]);
alpar@2390:     std::cout << "Make it sparser\n";
alpar@2402:     sparse2(ap["g"]);
alpar@2390:   }
alpar@2402:   else if(ap["tree"]) {
alpar@2390:     minTree();
alpar@2390:   }
alpar@2390:   
alpar@2390: 
alpar@2390:   std::cout << "Number of nodes    : " << countNodes(g) << std::endl;
alpar@2390:   std::cout << "Number of edges    : " << countUEdges(g) << std::endl;
alpar@2390:   double tlen=0;
alpar@2390:   for(UEdgeIt e(g);e!=INVALID;++e)
alpar@2390:     tlen+=sqrt((coords[g.source(e)]-coords[g.target(e)]).normSquare());
alpar@2390:   std::cout << "Total edge length  : " << tlen << std::endl;
alpar@2402:   if(ap["eps"])
alpar@2390:     graphToEps(g,prefix+".eps").
alpar@2390:       scale(600).nodeScale(.2).edgeWidthScale(.001).preScale(false).
alpar@2390:       coords(coords).run();
alpar@2390: 
alpar@2390:   UGraphWriter<ListUGraph>(prefix+".lgf",g).
alpar@2390:     writeNodeMap("coordinates_x",scaleMap(xMap(coords),600)).
alpar@2390:     writeNodeMap("coordinates_y",scaleMap(yMap(coords),600)).
alpar@2390:     run();
alpar@2390: }
alpar@2390: