lemon-0.x: tools/lgf-gen.cc@8450951a8e2d

     1 #include <lemon/list_graph.h>

     2 #include <lemon/graph_utils.h>

     3 #include <lemon/random.h>

     4 #include <lemon/dim2.h>

     5 #include <lemon/bfs.h>

     6 #include <lemon/counter.h>

     7 #include <lemon/suurballe.h>

     8 #include <lemon/graph_to_eps.h>

     9 #include <lemon/graph_writer.h>

    10 #include <lemon/arg_parser.h>

    11 #include <cmath>

    12 #include <algorithm>

    13 #include <lemon/unionfind.h>

    15 using namespace lemon;

    17 typedef dim2::Point<double> Point;

    19 UGRAPH_TYPEDEFS(ListUGraph);

    21 int N;

    22 int girth;

    24 ListUGraph g;

    26 std::vector<Node> nodes;

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

    29 int tsp_impr_num=0;

    31 const double EPSILON=1e-8;

    32 bool tsp_improve(Node u, Node v)

    33 {

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

    35   Node u2=u;

    36   Node v2=v;

    37   do {

    38     Node n;

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

    40     u2=v2;

    41     v2=n;

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

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

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

    45       {

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

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

    48 	g.addEdge(u2,u);

    49 	g.addEdge(v,v2);

    50 	tsp_impr_num++;

    51 	return true;

    52       }

    53   } while(v2!=u);

    54   return false;

    55 }

    57 bool tsp_improve(Node u)

    58 {

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

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

    61   return false;

    62 }

    64 void tsp_improve()

    65 {

    66   bool b;

    67   do {

    68     b=false;

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

    70       if(tsp_improve(n)) b=true;

    71   } while(b);

    72 }

    74 void tsp()

    75 {

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

    77   tsp_improve();

    78 }

    80 class Line

    81 {

    82 public:

    83   Point a;

    84   Point b;

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

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

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

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

    89 };

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

    92 {

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

    94   return os;

    95 }

    97 bool cross(Line a, Line b)

    98 {

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

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

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

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

   103 }

   105 struct Pedge

   106 {

   107   Node a;

   108   Node b;

   109   double len;

   110 };

   112 bool pedgeLess(Pedge a,Pedge b)

   113 {

   114   return a.len<b.len;

   115 }

   117 std::vector<UEdge> edges;

   119 void triangle()

   120 {

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

   122   std::vector<Pedge> pedges;

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

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

   125       {

   126 	Pedge p;

   127 	p.a=n;

   128 	p.b=m;

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

   130 	pedges.push_back(p);

   131       }

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

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

   134     {

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

   136       UEdgeIt e(g);

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

   138       UEdge ne;

   139       if(e==INVALID) {

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

   141 	edges.push_back(ne);

   142 	cnt++;

   143       }

   144     }

   145 }

   147 void sparse(int d)

   148 {

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

   150   Bfs<ListUGraph> bfs(g);

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

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

   153     {

   154       Node a=g.source(*ei);

   155       Node b=g.target(*ei);

   156       g.erase(*ei);

   157       bfs.run(a,b);

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

   159 	g.addEdge(a,b);

   160       else cnt++;

   161     }

   162 }

   164 void sparse2(int d)

   165 {

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

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

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

   169     {

   170       Node a=g.source(*ei);

   171       Node b=g.target(*ei);

   172       g.erase(*ei);

   173       ConstMap<Edge,int> cegy(1);

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

   175       int k=sur.run(2);

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

   177 	g.addEdge(a,b);

   178       else cnt++;

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

   180     }

   181 }

   183 void sparseTriangle(int d)

   184 {

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

   186   std::vector<Pedge> pedges;

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

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

   189       {

   190 	Pedge p;

   191 	p.a=n;

   192 	p.b=m;

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

   194 	pedges.push_back(p);

   195       }

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

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

   198     {

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

   200       UEdgeIt e(g);

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

   202       UEdge ne;

   203       if(e==INVALID) {

   204 	ConstMap<Edge,int> cegy(1);

   205 	Suurballe<ListUGraph,ConstMap<Edge,int> >

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

   207 	int k=sur.run(2);

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

   209 	  {

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

   211 	    edges.push_back(ne);

   212 	    cnt++;

   213 	  }

   214       }

   215     }

   216 }

   218 void minTree() {

   219   std::vector<Pedge> pedges;

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

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

   222       {

   223 	Pedge p;

   224 	p.a=n;

   225 	p.b=m;

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

   227 	pedges.push_back(p);

   228       }

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

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

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

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

   233     uf.insert(it);

   235   int en=0;

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

   237     {

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

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

   240 	en++;

   241 	if(en>=N-1) return;

   242       }

   243     }

   244 }

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

   249 {

   250   ArgParser ap(argc,argv);

   252   bool eps;

   253   bool disc_d, square_d, gauss_d;

   254   bool tsp_a,two_a,tree_a;

   255   int num_of_cities=1;

   256   double area=1;

   257   N=100;

   258   girth=10;

   259   std::string ndist("disc");

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

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

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

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

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

   265     .optionGroup("dist", "disc")

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

   267     .optionGroup("dist", "square")

   268     .option("gauss",

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

   270 	    gauss_d)

   271     .optionGroup("dist", "gauss")

   272 //     .mandatoryGroup("dist")

   273     .onlyOneGroup("dist")

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

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

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

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

   278     .optionGroup("alg","tree")

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

   280     .optionGroup("alg","tsp")

   281     .onlyOneGroup("alg")

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

   283     .run();

   285   std::string prefix;

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

   287     {

   288     case 0:

   289       prefix="lgf-gen-out";

   290       break;

   291     case 1:

   292       prefix=ap.files()[0];

   293       break;

   294     default:

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

   296       exit(1);

   297     }

   299   double sum_sizes=0;

   300   std::vector<double> sizes;

   301   std::vector<double> cum_sizes;

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

   303     {

   304       // 	sum_sizes+=rnd.exponential();

   305       double d=rnd();

   306       sum_sizes+=d;

   307       sizes.push_back(d);

   308       cum_sizes.push_back(sum_sizes);

   309     }

   310   int i=0;

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

   312     {

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

   314       if(gauss_d)

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

   316 	  Node n=g.addNode();

   317 	  nodes.push_back(n);

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

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

   320 	}

   321       else if(square_d)

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

   323 	  Node n=g.addNode();

   324 	  nodes.push_back(n);

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

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

   327 	}

   328       else if(disc_d || true)

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

   330 	  Node n=g.addNode();

   331 	  nodes.push_back(n);

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

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

   334 	}

   335     }

   337   if(tsp_a) {

   338     tsp();

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

   340   }

   341   else if(two_a) {

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

   343     //   triangle();

   344     sparseTriangle(girth);

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

   346     sparse2(girth);

   347   }

   348   else if(tree_a) {

   349     minTree();

   350   }

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

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

   355   double tlen=0;

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

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

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

   359   if(eps)

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

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

   362       coords(coords).run();

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

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

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

   367     run();

   368 }

author	alpar
	Sat, 03 Mar 2007 16:04:50 +0000
changeset 2390	8450951a8e2d
child 2391	14a343be7a5a
permissions	-rw-r--r--