1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/tools/lgf-gen.cc Sat Mar 03 16:04:50 2007 +0000
1.3 @@ -0,0 +1,369 @@
1.4 +#include <lemon/list_graph.h>
1.5 +#include <lemon/graph_utils.h>
1.6 +#include <lemon/random.h>
1.7 +#include <lemon/dim2.h>
1.8 +#include <lemon/bfs.h>
1.9 +#include <lemon/counter.h>
1.10 +#include <lemon/suurballe.h>
1.11 +#include <lemon/graph_to_eps.h>
1.12 +#include <lemon/graph_writer.h>
1.13 +#include <lemon/arg_parser.h>
1.14 +#include <cmath>
1.15 +#include <algorithm>
1.16 +#include <lemon/unionfind.h>
1.17 +
1.18 +using namespace lemon;
1.19 +
1.20 +typedef dim2::Point<double> Point;
1.21 +
1.22 +UGRAPH_TYPEDEFS(ListUGraph);
1.23 +
1.24 +int N;
1.25 +int girth;
1.26 +
1.27 +ListUGraph g;
1.28 +
1.29 +std::vector<Node> nodes;
1.30 +ListUGraph::NodeMap<Point> coords(g);
1.31 +
1.32 +int tsp_impr_num=0;
1.33 +
1.34 +const double EPSILON=1e-8;
1.35 +bool tsp_improve(Node u, Node v)
1.36 +{
1.37 + double luv=std::sqrt((coords[v]-coords[u]).normSquare());
1.38 + Node u2=u;
1.39 + Node v2=v;
1.40 + do {
1.41 + Node n;
1.42 + for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e);
1.43 + u2=v2;
1.44 + v2=n;
1.45 + if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON>
1.46 + std::sqrt((coords[u]-coords[u2]).normSquare())+
1.47 + std::sqrt((coords[v]-coords[v2]).normSquare()))
1.48 + {
1.49 + g.erase(findUEdge(g,u,v));
1.50 + g.erase(findUEdge(g,u2,v2));
1.51 + g.addEdge(u2,u);
1.52 + g.addEdge(v,v2);
1.53 + tsp_impr_num++;
1.54 + return true;
1.55 + }
1.56 + } while(v2!=u);
1.57 + return false;
1.58 +}
1.59 +
1.60 +bool tsp_improve(Node u)
1.61 +{
1.62 + for(IncEdgeIt e(g,u);e!=INVALID;++e)
1.63 + if(tsp_improve(u,g.runningNode(e))) return true;
1.64 + return false;
1.65 +}
1.66 +
1.67 +void tsp_improve()
1.68 +{
1.69 + bool b;
1.70 + do {
1.71 + b=false;
1.72 + for(NodeIt n(g);n!=INVALID;++n)
1.73 + if(tsp_improve(n)) b=true;
1.74 + } while(b);
1.75 +}
1.76 +
1.77 +void tsp()
1.78 +{
1.79 + for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);
1.80 + tsp_improve();
1.81 +}
1.82 +
1.83 +class Line
1.84 +{
1.85 +public:
1.86 + Point a;
1.87 + Point b;
1.88 + Line(Point _a,Point _b) :a(_a),b(_b) {}
1.89 + Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}
1.90 + Line(const Edge &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
1.91 + Line(const UEdge &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
1.92 +};
1.93 +
1.94 +inline std::ostream& operator<<(std::ostream &os, const Line &l)
1.95 +{
1.96 + os << l.a << "->" << l.b;
1.97 + return os;
1.98 +}
1.99 +
1.100 +bool cross(Line a, Line b)
1.101 +{
1.102 + Point ao=rot90(a.b-a.a);
1.103 + Point bo=rot90(b.b-b.a);
1.104 + return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 &&
1.105 + (bo*(a.a-b.a))*(bo*(a.b-b.a))<0;
1.106 +}
1.107 +
1.108 +struct Pedge
1.109 +{
1.110 + Node a;
1.111 + Node b;
1.112 + double len;
1.113 +};
1.114 +
1.115 +bool pedgeLess(Pedge a,Pedge b)
1.116 +{
1.117 + return a.len<b.len;
1.118 +}
1.119 +
1.120 +std::vector<UEdge> edges;
1.121 +
1.122 +void triangle()
1.123 +{
1.124 + Counter cnt("Number of edges added: ");
1.125 + std::vector<Pedge> pedges;
1.126 + for(NodeIt n(g);n!=INVALID;++n)
1.127 + for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
1.128 + {
1.129 + Pedge p;
1.130 + p.a=n;
1.131 + p.b=m;
1.132 + p.len=(coords[m]-coords[n]).normSquare();
1.133 + pedges.push_back(p);
1.134 + }
1.135 + std::sort(pedges.begin(),pedges.end(),pedgeLess);
1.136 + for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
1.137 + {
1.138 + Line li(pi->a,pi->b);
1.139 + UEdgeIt e(g);
1.140 + for(;e!=INVALID && !cross(e,li);++e) ;
1.141 + UEdge ne;
1.142 + if(e==INVALID) {
1.143 + ne=g.addEdge(pi->a,pi->b);
1.144 + edges.push_back(ne);
1.145 + cnt++;
1.146 + }
1.147 + }
1.148 +}
1.149 +
1.150 +void sparse(int d)
1.151 +{
1.152 + Counter cnt("Number of edges removed: ");
1.153 + Bfs<ListUGraph> bfs(g);
1.154 + for(std::vector<UEdge>::reverse_iterator ei=edges.rbegin();
1.155 + ei!=edges.rend();++ei)
1.156 + {
1.157 + Node a=g.source(*ei);
1.158 + Node b=g.target(*ei);
1.159 + g.erase(*ei);
1.160 + bfs.run(a,b);
1.161 + if(bfs.predEdge(b)==INVALID || bfs.dist(b)>d)
1.162 + g.addEdge(a,b);
1.163 + else cnt++;
1.164 + }
1.165 +}
1.166 +
1.167 +void sparse2(int d)
1.168 +{
1.169 + Counter cnt("Number of edges removed: ");
1.170 + for(std::vector<UEdge>::reverse_iterator ei=edges.rbegin();
1.171 + ei!=edges.rend();++ei)
1.172 + {
1.173 + Node a=g.source(*ei);
1.174 + Node b=g.target(*ei);
1.175 + g.erase(*ei);
1.176 + ConstMap<Edge,int> cegy(1);
1.177 + Suurballe<ListUGraph,ConstMap<Edge,int> > sur(g,cegy,a,b);
1.178 + int k=sur.run(2);
1.179 + if(k<2 || sur.totalLength()>d)
1.180 + g.addEdge(a,b);
1.181 + else cnt++;
1.182 +// else std::cout << "Remove edge " << g.id(a) << "-" << g.id(b) << '\n';
1.183 + }
1.184 +}
1.185 +
1.186 +void sparseTriangle(int d)
1.187 +{
1.188 + Counter cnt("Number of edges added: ");
1.189 + std::vector<Pedge> pedges;
1.190 + for(NodeIt n(g);n!=INVALID;++n)
1.191 + for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
1.192 + {
1.193 + Pedge p;
1.194 + p.a=n;
1.195 + p.b=m;
1.196 + p.len=(coords[m]-coords[n]).normSquare();
1.197 + pedges.push_back(p);
1.198 + }
1.199 + std::sort(pedges.begin(),pedges.end(),pedgeLess);
1.200 + for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
1.201 + {
1.202 + Line li(pi->a,pi->b);
1.203 + UEdgeIt e(g);
1.204 + for(;e!=INVALID && !cross(e,li);++e) ;
1.205 + UEdge ne;
1.206 + if(e==INVALID) {
1.207 + ConstMap<Edge,int> cegy(1);
1.208 + Suurballe<ListUGraph,ConstMap<Edge,int> >
1.209 + sur(g,cegy,pi->a,pi->b);
1.210 + int k=sur.run(2);
1.211 + if(k<2 || sur.totalLength()>d)
1.212 + {
1.213 + ne=g.addEdge(pi->a,pi->b);
1.214 + edges.push_back(ne);
1.215 + cnt++;
1.216 + }
1.217 + }
1.218 + }
1.219 +}
1.220 +
1.221 +void minTree() {
1.222 + std::vector<Pedge> pedges;
1.223 + for(NodeIt n(g);n!=INVALID;++n)
1.224 + for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
1.225 + {
1.226 + Pedge p;
1.227 + p.a=n;
1.228 + p.b=m;
1.229 + p.len=(coords[m]-coords[n]).normSquare();
1.230 + pedges.push_back(p);
1.231 + }
1.232 + std::sort(pedges.begin(),pedges.end(),pedgeLess);
1.233 + ListUGraph::NodeMap<int> comp(g);
1.234 + UnionFind<ListUGraph::NodeMap<int> > uf(comp);
1.235 + for (NodeIt it(g); it != INVALID; ++it)
1.236 + uf.insert(it);
1.237 +
1.238 + int en=0;
1.239 + for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
1.240 + {
1.241 + if ( uf.join(pi->a,pi->b) ) {
1.242 + g.addEdge(pi->a,pi->b);
1.243 + en++;
1.244 + if(en>=N-1) return;
1.245 + }
1.246 + }
1.247 +}
1.248 +
1.249 +
1.250 +
1.251 +int main(int argc,char **argv)
1.252 +{
1.253 + ArgParser ap(argc,argv);
1.254 +
1.255 + bool eps;
1.256 + bool disc_d, square_d, gauss_d;
1.257 + bool tsp_a,two_a,tree_a;
1.258 + int num_of_cities=1;
1.259 + double area=1;
1.260 + N=100;
1.261 + girth=10;
1.262 + std::string ndist("disc");
1.263 + ap.option("n", "Number of nodes (default is 100)", N)
1.264 + .option("g", "Girth parameter (default is 10)", girth)
1.265 + .option("cities", "Number of cities (default is 1)", num_of_cities)
1.266 + .option("area", "Full relative area of the cities (default is 1)", area)
1.267 + .option("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d)
1.268 + .optionGroup("dist", "disc")
1.269 + .option("square", "Nodes are evenly distributed on a unit square", square_d)
1.270 + .optionGroup("dist", "square")
1.271 + .option("gauss",
1.272 + "Nodes are located according to a two-dim gauss distribution",
1.273 + gauss_d)
1.274 + .optionGroup("dist", "gauss")
1.275 +// .mandatoryGroup("dist")
1.276 + .onlyOneGroup("dist")
1.277 + .option("eps", "Also generate .eps output (prefix.eps)",eps)
1.278 + .option("2con", "Create a two connected planar graph",two_a)
1.279 + .optionGroup("alg","2con")
1.280 + .option("tree", "Create a min. cost spanning tree",tree_a)
1.281 + .optionGroup("alg","tree")
1.282 + .option("tsp", "Create a TSP tour",tsp_a)
1.283 + .optionGroup("alg","tsp")
1.284 + .onlyOneGroup("alg")
1.285 + .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
1.286 + .run();
1.287 +
1.288 + std::string prefix;
1.289 + switch(ap.files().size())
1.290 + {
1.291 + case 0:
1.292 + prefix="lgf-gen-out";
1.293 + break;
1.294 + case 1:
1.295 + prefix=ap.files()[0];
1.296 + break;
1.297 + default:
1.298 + std::cerr << "\nAt most one prefix can be given\n\n";
1.299 + exit(1);
1.300 + }
1.301 +
1.302 + double sum_sizes=0;
1.303 + std::vector<double> sizes;
1.304 + std::vector<double> cum_sizes;
1.305 + for(int s=0;s<num_of_cities;s++)
1.306 + {
1.307 + // sum_sizes+=rnd.exponential();
1.308 + double d=rnd();
1.309 + sum_sizes+=d;
1.310 + sizes.push_back(d);
1.311 + cum_sizes.push_back(sum_sizes);
1.312 + }
1.313 + int i=0;
1.314 + for(int s=0;s<num_of_cities;s++)
1.315 + {
1.316 + Point center=(num_of_cities==1?Point(0,0):rnd.disc());
1.317 + if(gauss_d)
1.318 + for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
1.319 + Node n=g.addNode();
1.320 + nodes.push_back(n);
1.321 + coords[n]=center+rnd.gauss2()*area*
1.322 + std::sqrt(sizes[s]/sum_sizes);
1.323 + }
1.324 + else if(square_d)
1.325 + for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
1.326 + Node n=g.addNode();
1.327 + nodes.push_back(n);
1.328 + coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
1.329 + std::sqrt(sizes[s]/sum_sizes);
1.330 + }
1.331 + else if(disc_d || true)
1.332 + for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
1.333 + Node n=g.addNode();
1.334 + nodes.push_back(n);
1.335 + coords[n]=center+rnd.disc()*area*
1.336 + std::sqrt(sizes[s]/sum_sizes);
1.337 + }
1.338 + }
1.339 +
1.340 + if(tsp_a) {
1.341 + tsp();
1.342 + std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
1.343 + }
1.344 + else if(two_a) {
1.345 + std::cout << "Make triangles\n";
1.346 + // triangle();
1.347 + sparseTriangle(girth);
1.348 + std::cout << "Make it sparser\n";
1.349 + sparse2(girth);
1.350 + }
1.351 + else if(tree_a) {
1.352 + minTree();
1.353 + }
1.354 +
1.355 +
1.356 + std::cout << "Number of nodes : " << countNodes(g) << std::endl;
1.357 + std::cout << "Number of edges : " << countUEdges(g) << std::endl;
1.358 + double tlen=0;
1.359 + for(UEdgeIt e(g);e!=INVALID;++e)
1.360 + tlen+=sqrt((coords[g.source(e)]-coords[g.target(e)]).normSquare());
1.361 + std::cout << "Total edge length : " << tlen << std::endl;
1.362 + if(eps)
1.363 + graphToEps(g,prefix+".eps").
1.364 + scale(600).nodeScale(.2).edgeWidthScale(.001).preScale(false).
1.365 + coords(coords).run();
1.366 +
1.367 + UGraphWriter<ListUGraph>(prefix+".lgf",g).
1.368 + writeNodeMap("coordinates_x",scaleMap(xMap(coords),600)).
1.369 + writeNodeMap("coordinates_y",scaleMap(yMap(coords),600)).
1.370 + run();
1.371 +}
1.372 +