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

     1 /* -*- C++ -*-

     2  *

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

     4  *

     5  * Copyright (C) 2003-2007

     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

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

     8  *

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

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

    11  * precise terms see the accompanying LICENSE file.

    12  *

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

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

    15  * purpose.

    16  *

    17  */

    19 #include <lemon/list_graph.h>

    20 #include <lemon/graph_utils.h>

    21 #include <lemon/random.h>

    22 #include <lemon/dim2.h>

    23 #include <lemon/bfs.h>

    24 #include <lemon/counter.h>

    25 #include <lemon/suurballe.h>

    26 #include <lemon/graph_to_eps.h>

    27 #include <lemon/graph_writer.h>

    28 #include <lemon/arg_parser.h>

    29 #include <cmath>

    30 #include <algorithm>

    31 #include <lemon/unionfind.h>

    32 #include <lemon/time_measure.h>

    34 using namespace lemon;

    36 typedef dim2::Point<double> Point;

    38 UGRAPH_TYPEDEFS(ListUGraph);

    40 bool progress=true;

    42 int N;

    43 // int girth;

    45 ListUGraph g;

    47 std::vector<Node> nodes;

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

    50 int tsp_impr_num=0;

    52 const double EPSILON=1e-8;

    53 bool tsp_improve(Node u, Node v)

    54 {

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

    56   Node u2=u;

    57   Node v2=v;

    58   do {

    59     Node n;

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

    61     u2=v2;

    62     v2=n;

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

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

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

    66       {

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

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

    69 	g.addEdge(u2,u);

    70 	g.addEdge(v,v2);

    71 	tsp_impr_num++;

    72 	return true;

    73       }

    74   } while(v2!=u);

    75   return false;

    76 }

    78 bool tsp_improve(Node u)

    79 {

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

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

    82   return false;

    83 }

    85 void tsp_improve()

    86 {

    87   bool b;

    88   do {

    89     b=false;

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

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

    92   } while(b);

    93 }

    95 void tsp()

    96 {

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

    98   tsp_improve();

    99 }

   101 class Line

   102 {

   103 public:

   104   Point a;

   105   Point b;

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

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

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

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

   110 };

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

   113 {

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

   115   return os;

   116 }

   118 bool cross(Line a, Line b)

   119 {

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

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

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

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

   124 }

   126 struct Pedge

   127 {

   128   Node a;

   129   Node b;

   130   double len;

   131 };

   133 bool pedgeLess(Pedge a,Pedge b)

   134 {

   135   return a.len<b.len;

   136 }

   138 std::vector<UEdge> edges;

   140 void triangle()

   141 {

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

   143   std::vector<Pedge> pedges;

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

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

   146       {

   147 	Pedge p;

   148 	p.a=n;

   149 	p.b=m;

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

   151 	pedges.push_back(p);

   152       }

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

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

   155     {

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

   157       UEdgeIt e(g);

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

   159       UEdge ne;

   160       if(e==INVALID) {

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

   162 	edges.push_back(ne);

   163 	cnt++;

   164       }

   165     }

   166 }

   168 void sparse(int d)

   169 {

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

   171   Bfs<ListUGraph> bfs(g);

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

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

   174     {

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

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

   177       g.erase(*ei);

   178       bfs.run(a,b);

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

   180 	g.addEdge(a,b);

   181       else cnt++;

   182     }

   183 }

   185 void sparse2(int d)

   186 {

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

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

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

   190     {

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

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

   193       g.erase(*ei);

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

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

   196       int k=sur.run(2);

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

   198 	g.addEdge(a,b);

   199       else cnt++;

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

   201     }

   202 }

   204 void sparseTriangle(int d)

   205 {

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

   207   std::vector<Pedge> pedges;

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

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

   210       {

   211 	Pedge p;

   212 	p.a=n;

   213 	p.b=m;

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

   215 	pedges.push_back(p);

   216       }

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

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

   219     {

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

   221       UEdgeIt e(g);

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

   223       UEdge ne;

   224       if(e==INVALID) {

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

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

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

   228 	int k=sur.run(2);

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

   230 	  {

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

   232 	    edges.push_back(ne);

   233 	    cnt++;

   234 	  }

   235       }

   236     }

   237 }

   239 void minTree() {

   240   int en=0;

   241   int pr=0;

   242   std::vector<Pedge> pedges;

   243   Timer T;

   244   std::cout << T.realTime() << "s: Setting up the edges...\n";

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

   246     {

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

   248 	{

   249 	  Pedge p;

   250 	  p.a=n;

   251 	  p.b=m;

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

   253 	  pedges.push_back(p);

   254 	}

   255       en++;

   256       if(progress && en>=pr*double(N)/100)

   257 	{

   258 	  std::cout << pr << "%  \r" << std::flush;

   259 	  pr++;

   260 	}

   261     }

   262   std::cout << T.realTime() << "s: Sorting the edges...\n";

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

   264   std::cout << T.realTime() << "s: Creating the tree...\n";

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

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

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

   268     uf.insert(it);

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

   270     {

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

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

   273 	if(en>=N-1)

   274 	  break;

   275       }

   276     }

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

   278 }

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

   283 {

   284   ArgParser ap(argc,argv);

   286 //   bool eps;

   287   bool disc_d, square_d, gauss_d;

   288 //   bool tsp_a,two_a,tree_a;

   289   int num_of_cities=1;

   290   double area=1;

   291   N=100;

   292 //   girth=10;

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

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

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

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

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

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

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

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

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

   302     .refOption("gauss",

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

   304 	    gauss_d)

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

   306 //     .mandatoryGroup("dist")

   307     .onlyOneGroup("dist")

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

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

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

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

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

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

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

   315     .onlyOneGroup("alg")

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

   317     .run();

   319   std::string prefix;

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

   321     {

   322     case 0:

   323       prefix="lgf-gen-out";

   324       break;

   325     case 1:

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

   327       break;

   328     default:

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

   330       exit(1);

   331     }

   333   double sum_sizes=0;

   334   std::vector<double> sizes;

   335   std::vector<double> cum_sizes;

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

   337     {

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

   339       double d=rnd();

   340       sum_sizes+=d;

   341       sizes.push_back(d);

   342       cum_sizes.push_back(sum_sizes);

   343     }

   344   int i=0;

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

   346     {

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

   348       if(gauss_d)

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

   350 	  Node n=g.addNode();

   351 	  nodes.push_back(n);

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

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

   354 	}

   355       else if(square_d)

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

   357 	  Node n=g.addNode();

   358 	  nodes.push_back(n);

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

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

   361 	}

   362       else if(disc_d || true)

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

   364 	  Node n=g.addNode();

   365 	  nodes.push_back(n);

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

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

   368 	}

   369     }

   371   if(ap["tsp"]) {

   372     tsp();

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

   374   }

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

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

   377     //   triangle();

   378     sparseTriangle(ap["g"]);

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

   380     sparse2(ap["g"]);

   381   }

   382   else if(ap["tree"]) {

   383     minTree();

   384   }

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

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

   389   double tlen=0;

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

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

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

   393   if(ap["eps"])

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

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

   396       coords(coords).run();

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

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

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

   401     run();

   402 }

author	alpar
	Tue, 13 Mar 2007 16:32:35 +0000
changeset 2409	fe0a8fe16271
parent 2391	14a343be7a5a
child 2410	fe46b61da4e3
permissions	-rw-r--r--