3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2007
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
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.
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
20 ///\brief Special plane graph generator.
22 ///Graph generator application for various types of plane graphs.
25 /// This program converts various DIMACS formats to the LEMON Graph Format
29 /// ./tools/lgf-gen [-2con|-tree|-tsp|-tsp2|-dela] [-disc|-square|-gauss]
30 /// [-rand|-seed int] [--help|-h|-help] [-area num] [-cities int] [-dir]
31 /// [-eps] [-g int] [-n int] [prefix]
34 /// Prefix of the output files. Default is 'lgf-gen-out'
36 /// Print a short help message
38 /// Create a two connected planar graph
40 /// Full relative area of the cities (default is 1)
42 /// Number of cities (default is 1)
44 /// Delaunay triangulation graph
46 /// Directed graph is generated (each edges are replaced by two directed ones)
48 /// Nodes are evenly distributed on a unit disc (default)
50 /// Also generate .eps output (prefix.eps)
52 /// Girth parameter (default is 10)
54 /// Nodes are located according to a two-dim gauss distribution
56 /// Number of nodes (default is 100)
58 /// Use time seed for random number generator
62 /// Nodes are evenly distributed on a unit square
64 /// Create a min. cost spanning tree
68 /// Create a TSP tour (tree based)
70 /// \image html plane_tree.png
71 /// \image latex plane_tree.eps "Eucledian spanning tree" width=\textwidth
75 #include <lemon/list_graph.h>
76 #include <lemon/graph_utils.h>
77 #include <lemon/random.h>
78 #include <lemon/dim2.h>
79 #include <lemon/bfs.h>
80 #include <lemon/counter.h>
81 #include <lemon/suurballe.h>
82 #include <lemon/graph_to_eps.h>
83 #include <lemon/graph_writer.h>
84 #include <lemon/arg_parser.h>
85 #include <lemon/euler.h>
88 #include <lemon/kruskal.h>
89 #include <lemon/time_measure.h>
91 using namespace lemon;
93 typedef dim2::Point<double> Point;
95 UGRAPH_TYPEDEFS(ListUGraph);
104 std::vector<Node> nodes;
105 ListUGraph::NodeMap<Point> coords(g);
110 for(UEdgeIt e(g);e!=INVALID;++e)
111 tlen+=sqrt((coords[g.source(e)]-coords[g.target(e)]).normSquare());
117 const double EPSILON=1e-8;
118 bool tsp_improve(Node u, Node v)
120 double luv=std::sqrt((coords[v]-coords[u]).normSquare());
125 for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e);
128 if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON>
129 std::sqrt((coords[u]-coords[u2]).normSquare())+
130 std::sqrt((coords[v]-coords[v2]).normSquare()))
132 g.erase(findUEdge(g,u,v));
133 g.erase(findUEdge(g,u2,v2));
143 bool tsp_improve(Node u)
145 for(IncEdgeIt e(g,u);e!=INVALID;++e)
146 if(tsp_improve(u,g.runningNode(e))) return true;
155 for(NodeIt n(g);n!=INVALID;++n)
156 if(tsp_improve(n)) b=true;
162 for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);
171 Line(Point _a,Point _b) :a(_a),b(_b) {}
172 Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}
173 Line(const Edge &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
174 Line(const UEdge &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
177 inline std::ostream& operator<<(std::ostream &os, const Line &l)
179 os << l.a << "->" << l.b;
183 bool cross(Line a, Line b)
185 Point ao=rot90(a.b-a.a);
186 Point bo=rot90(b.b-b.a);
187 return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 &&
188 (bo*(a.a-b.a))*(bo*(a.b-b.a))<0;
198 bool pedgeLess(Pedge a,Pedge b)
203 std::vector<UEdge> edges;
205 namespace _delaunay_bits {
208 int prev, curr, next;
210 Part(int p, int c, int n) : prev(p), curr(c), next(n) {}
213 inline std::ostream& operator<<(std::ostream& os, const Part& part) {
214 os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';
218 inline double circle_point(const Point& p, const Point& q, const Point& r) {
219 double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);
220 if (a == 0) return std::numeric_limits<double>::quiet_NaN();
222 double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +
223 (q.x * q.x + q.y * q.y) * (r.y - p.y) +
224 (r.x * r.x + r.y * r.y) * (p.y - q.y);
226 double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +
227 (q.x * q.x + q.y * q.y) * (r.x - p.x) +
228 (r.x * r.x + r.y * r.y) * (p.x - q.x);
230 double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +
231 (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +
232 (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);
234 return d / (2 * a) + sqrt((d * d + e * e) / (4 * a * a) + f / a);
237 inline bool circle_form(const Point& p, const Point& q, const Point& r) {
238 return rot90(q - p) * (r - q) < 0.0;
241 inline double intersection(const Point& p, const Point& q, double sx) {
242 const double epsilon = 1e-8;
244 if (p.x == q.x) return (p.y + q.y) / 2.0;
246 if (sx < p.x + epsilon) return p.y;
247 if (sx < q.x + epsilon) return q.y;
249 double a = q.x - p.x;
250 double b = (q.x - sx) * p.y - (p.x - sx) * q.y;
251 double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();
252 return (b - sqrt(d)) / a;
258 YLess(const std::vector<Point>& points, double& sweep)
259 : _points(points), _sweep(sweep) {}
261 bool operator()(const Part& l, const Part& r) const {
262 const double epsilon = 1e-8;
264 // std::cerr << l << " vs " << r << std::endl;
265 double lbx = l.prev != -1 ?
266 intersection(_points[l.prev], _points[l.curr], _sweep) :
267 - std::numeric_limits<double>::infinity();
268 double rbx = r.prev != -1 ?
269 intersection(_points[r.prev], _points[r.curr], _sweep) :
270 - std::numeric_limits<double>::infinity();
271 double lex = l.next != -1 ?
272 intersection(_points[l.curr], _points[l.next], _sweep) :
273 std::numeric_limits<double>::infinity();
274 double rex = r.next != -1 ?
275 intersection(_points[r.curr], _points[r.next], _sweep) :
276 std::numeric_limits<double>::infinity();
278 if (lbx > lex) std::swap(lbx, lex);
279 if (rbx > rex) std::swap(rbx, rex);
281 if (lex < epsilon + rex && lbx + epsilon < rex) return true;
282 if (rex < epsilon + lex && rbx + epsilon < lex) return false;
286 const std::vector<Point>& _points;
292 typedef std::multimap<double, BeachIt> SpikeHeap;
294 typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;
299 BeachIt(Beach::iterator iter) : it(iter) {}
304 inline void delaunay() {
305 Counter cnt("Number of edges added: ");
307 using namespace _delaunay_bits;
309 typedef _delaunay_bits::Part Part;
310 typedef std::vector<std::pair<double, int> > SiteHeap;
313 std::vector<Point> points;
314 std::vector<Node> nodes;
316 for (NodeIt it(g); it != INVALID; ++it) {
318 points.push_back(coords[it]);
321 SiteHeap siteheap(points.size());
326 for (int i = 0; i < int(siteheap.size()); ++i) {
327 siteheap[i] = std::make_pair(points[i].x, i);
330 std::sort(siteheap.begin(), siteheap.end());
331 sweep = siteheap.front().first;
333 YLess yless(points, sweep);
338 std::set<std::pair<int, int> > edges;
344 while (siteindex < int(siteheap.size()) &&
345 siteheap[0].first == siteheap[siteindex].first) {
346 front.push_back(std::make_pair(points[siteheap[siteindex].second].y,
347 siteheap[siteindex].second));
351 std::sort(front.begin(), front.end());
353 for (int i = 0; i < int(front.size()); ++i) {
354 int prev = (i == 0 ? -1 : front[i - 1].second);
355 int curr = front[i].second;
356 int next = (i + 1 == int(front.size()) ? -1 : front[i + 1].second);
358 beach.insert(std::make_pair(Part(prev, curr, next),
363 while (siteindex < int(points.size()) || !spikeheap.empty()) {
365 SpikeHeap::iterator spit = spikeheap.begin();
367 if (siteindex < int(points.size()) &&
368 (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {
369 int site = siteheap[siteindex].second;
370 sweep = siteheap[siteindex].first;
372 Beach::iterator bit = beach.upper_bound(Part(site, site, site));
374 if (bit->second != spikeheap.end()) {
375 spikeheap.erase(bit->second);
378 int prev = bit->first.prev;
379 int curr = bit->first.curr;
380 int next = bit->first.next;
384 SpikeHeap::iterator pit = spikeheap.end();
386 circle_form(points[prev], points[curr], points[site])) {
387 double x = circle_point(points[prev], points[curr], points[site]);
388 pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
390 beach.insert(std::make_pair(Part(prev, curr, site), pit));
392 beach.insert(std::make_pair(Part(prev, curr, site), pit));
395 beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
397 SpikeHeap::iterator nit = spikeheap.end();
399 circle_form(points[site], points[curr],points[next])) {
400 double x = circle_point(points[site], points[curr], points[next]);
401 nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
403 beach.insert(std::make_pair(Part(site, curr, next), nit));
405 beach.insert(std::make_pair(Part(site, curr, next), nit));
412 Beach::iterator bit = spit->second.it;
414 int prev = bit->first.prev;
415 int curr = bit->first.curr;
416 int next = bit->first.next;
419 std::pair<int, int> edge;
422 std::make_pair(prev, curr) : std::make_pair(curr, prev);
424 if (edges.find(edge) == edges.end()) {
426 g.addEdge(nodes[prev], nodes[curr]);
431 std::make_pair(curr, next) : std::make_pair(next, curr);
433 if (edges.find(edge) == edges.end()) {
435 g.addEdge(nodes[curr], nodes[next]);
440 Beach::iterator pbit = bit; --pbit;
441 int ppv = pbit->first.prev;
442 Beach::iterator nbit = bit; ++nbit;
443 int nnt = nbit->first.next;
445 if (bit->second != spikeheap.end()) spikeheap.erase(bit->second);
446 if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second);
447 if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second);
453 SpikeHeap::iterator pit = spikeheap.end();
454 if (ppv != -1 && ppv != next &&
455 circle_form(points[ppv], points[prev], points[next])) {
456 double x = circle_point(points[ppv], points[prev], points[next]);
457 if (x < sweep) x = sweep;
458 pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
460 beach.insert(std::make_pair(Part(ppv, prev, next), pit));
462 beach.insert(std::make_pair(Part(ppv, prev, next), pit));
465 SpikeHeap::iterator nit = spikeheap.end();
466 if (nnt != -1 && prev != nnt &&
467 circle_form(points[prev], points[next], points[nnt])) {
468 double x = circle_point(points[prev], points[next], points[nnt]);
469 if (x < sweep) x = sweep;
470 nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
472 beach.insert(std::make_pair(Part(prev, next, nnt), nit));
474 beach.insert(std::make_pair(Part(prev, next, nnt), nit));
480 for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
481 int curr = it->first.curr;
482 int next = it->first.next;
484 if (next == -1) continue;
486 std::pair<int, int> edge;
489 std::make_pair(curr, next) : std::make_pair(next, curr);
491 if (edges.find(edge) == edges.end()) {
493 g.addEdge(nodes[curr], nodes[next]);
501 Counter cnt("Number of edges removed: ");
502 Bfs<ListUGraph> bfs(g);
503 for(std::vector<UEdge>::reverse_iterator ei=edges.rbegin();
504 ei!=edges.rend();++ei)
506 Node a=g.source(*ei);
507 Node b=g.target(*ei);
510 if(bfs.predEdge(b)==INVALID || bfs.dist(b)>d)
518 Counter cnt("Number of edges removed: ");
519 for(std::vector<UEdge>::reverse_iterator ei=edges.rbegin();
520 ei!=edges.rend();++ei)
522 Node a=g.source(*ei);
523 Node b=g.target(*ei);
525 ConstMap<Edge,int> cegy(1);
526 Suurballe<ListUGraph,ConstMap<Edge,int> > sur(g,cegy,a,b);
528 if(k<2 || sur.totalLength()>d)
531 // else std::cout << "Remove edge " << g.id(a) << "-" << g.id(b) << '\n';
535 void sparseTriangle(int d)
537 Counter cnt("Number of edges added: ");
538 std::vector<Pedge> pedges;
539 for(NodeIt n(g);n!=INVALID;++n)
540 for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
545 p.len=(coords[m]-coords[n]).normSquare();
548 std::sort(pedges.begin(),pedges.end(),pedgeLess);
549 for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
551 Line li(pi->a,pi->b);
553 for(;e!=INVALID && !cross(e,li);++e) ;
556 ConstMap<Edge,int> cegy(1);
557 Suurballe<ListUGraph,ConstMap<Edge,int> >
558 sur(g,cegy,pi->a,pi->b);
560 if(k<2 || sur.totalLength()>d)
562 ne=g.addEdge(pi->a,pi->b);
570 template <typename UGraph, typename CoordMap>
571 class LengthSquareMap {
573 typedef typename UGraph::UEdge Key;
574 typedef typename CoordMap::Value::Value Value;
576 LengthSquareMap(const UGraph& ugraph, const CoordMap& coords)
577 : _ugraph(ugraph), _coords(coords) {}
579 Value operator[](const Key& key) const {
580 return (_coords[_ugraph.target(key)] -
581 _coords[_ugraph.source(key)]).normSquare();
586 const UGraph& _ugraph;
587 const CoordMap& _coords;
591 std::vector<Pedge> pedges;
593 std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
595 std::cout << T.realTime() << "s: Calculating spanning tree...\n";
596 LengthSquareMap<ListUGraph, ListUGraph::NodeMap<Point> > ls(g, coords);
597 ListUGraph::UEdgeMap<bool> tree(g);
598 kruskal(g, ls, tree);
599 std::cout << T.realTime() << "s: Removing non tree edges...\n";
600 std::vector<UEdge> remove;
601 for (UEdgeIt e(g); e != INVALID; ++e) {
602 if (!tree[e]) remove.push_back(e);
604 for(int i = 0; i < int(remove.size()); ++i) {
607 std::cout << T.realTime() << "s: Done\n";
612 std::cout << "Find a tree..." << std::endl;
616 std::cout << "Total edge length (tree) : " << totalLen() << std::endl;
618 std::cout << "Make it Euler..." << std::endl;
621 std::vector<Node> leafs;
622 for(NodeIt n(g);n!=INVALID;++n)
623 if(countIncEdges(g,n)%2==1) leafs.push_back(n);
625 // for(unsigned int i=0;i<leafs.size();i+=2)
626 // g.addEdge(leafs[i],leafs[i+1]);
628 std::vector<Pedge> pedges;
629 for(unsigned int i=0;i<leafs.size()-1;i++)
630 for(unsigned int j=i+1;j<leafs.size();j++)
637 p.len=(coords[m]-coords[n]).normSquare();
640 std::sort(pedges.begin(),pedges.end(),pedgeLess);
641 for(unsigned int i=0;i<pedges.size();i++)
642 if(countIncEdges(g,pedges[i].a)%2 &&
643 countIncEdges(g,pedges[i].b)%2)
644 g.addEdge(pedges[i].a,pedges[i].b);
647 for(NodeIt n(g);n!=INVALID;++n)
648 if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
649 std::cout << "GEBASZ!!!" << std::endl;
651 for(UEdgeIt e(g);e!=INVALID;++e)
652 if(g.source(e)==g.target(e))
653 std::cout << "LOOP GEBASZ!!!" << std::endl;
655 std::cout << "Number of edges : " << countUEdges(g) << std::endl;
657 std::cout << "Total edge length (euler) : " << totalLen() << std::endl;
659 ListUGraph::UEdgeMap<Edge> enext(g);
661 UEulerIt<ListUGraph> e(g);
664 // std::cout << "Tour edge: " << g.id(UEdge(e)) << std::endl;
665 for(++e;e!=INVALID;++e)
667 // std::cout << "Tour edge: " << g.id(UEdge(e)) << std::endl;
674 std::cout << "Creating a tour from that..." << std::endl;
676 int nnum = countNodes(g);
677 int ednum = countUEdges(g);
679 for(Edge p=enext[UEdgeIt(g)];ednum>nnum;p=enext[p])
681 // std::cout << "Checking edge " << g.id(p) << std::endl;
685 Node n1=g.oppositeNode(n2,e);
686 Node n3=g.oppositeNode(n2,f);
687 if(countIncEdges(g,n2)>2)
689 // std::cout << "Remove an Edge" << std::endl;
695 Edge ne=g.direct(g.addEdge(n1,n3),n1);
707 std::cout << "Total edge length (tour) : " << totalLen() << std::endl;
709 std::cout << "2-opt the tour..." << std::endl;
713 std::cout << "Total edge length (2-opt tour) : " << totalLen() << std::endl;
717 int main(int argc,const char **argv)
719 ArgParser ap(argc,argv);
722 bool disc_d, square_d, gauss_d;
723 // bool tsp_a,two_a,tree_a;
728 std::string ndist("disc");
729 ap.refOption("n", "Number of nodes (default is 100)", N)
730 .intOption("g", "Girth parameter (default is 10)", 10)
731 .refOption("cities", "Number of cities (default is 1)", num_of_cities)
732 .refOption("area", "Full relative area of the cities (default is 1)", area)
733 .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d)
734 .optionGroup("dist", "disc")
735 .refOption("square", "Nodes are evenly distributed on a unit square", square_d)
736 .optionGroup("dist", "square")
738 "Nodes are located according to a two-dim gauss distribution",
740 .optionGroup("dist", "gauss")
741 // .mandatoryGroup("dist")
742 .onlyOneGroup("dist")
743 .boolOption("eps", "Also generate .eps output (prefix.eps)")
744 .boolOption("dir", "Directed graph is generated (each edges are replaced by two directed ones)")
745 .boolOption("2con", "Create a two connected planar graph")
746 .optionGroup("alg","2con")
747 .boolOption("tree", "Create a min. cost spanning tree")
748 .optionGroup("alg","tree")
749 .boolOption("tsp", "Create a TSP tour")
750 .optionGroup("alg","tsp")
751 .boolOption("tsp2", "Create a TSP tour (tree based)")
752 .optionGroup("alg","tsp2")
753 .boolOption("dela", "Delaunay triangulation graph")
754 .optionGroup("alg","dela")
756 .boolOption("rand", "Use time seed for random number generator")
757 .optionGroup("rand", "rand")
758 .intOption("seed", "Random seed", -1)
759 .optionGroup("rand", "seed")
760 .onlyOneGroup("rand")
761 .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
766 std::cout << "Random number seed: " << seed << std::endl;
769 if (ap.given("seed")) {
770 int seed = ap["seed"];
771 std::cout << "Random number seed: " << seed << std::endl;
776 switch(ap.files().size())
779 prefix="lgf-gen-out";
782 prefix=ap.files()[0];
785 std::cerr << "\nAt most one prefix can be given\n\n";
790 std::vector<double> sizes;
791 std::vector<double> cum_sizes;
792 for(int s=0;s<num_of_cities;s++)
794 // sum_sizes+=rnd.exponential();
798 cum_sizes.push_back(sum_sizes);
801 for(int s=0;s<num_of_cities;s++)
803 Point center=(num_of_cities==1?Point(0,0):rnd.disc());
805 for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
808 coords[n]=center+rnd.gauss2()*area*
809 std::sqrt(sizes[s]/sum_sizes);
812 for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
815 coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
816 std::sqrt(sizes[s]/sum_sizes);
818 else if(disc_d || true)
819 for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
822 coords[n]=center+rnd.disc()*area*
823 std::sqrt(sizes[s]/sum_sizes);
827 // for (ListUGraph::NodeIt n(g); n != INVALID; ++n) {
828 // std::cerr << coords[n] << std::endl;
833 std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
837 std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
839 else if(ap["2con"]) {
840 std::cout << "Make triangles\n";
842 sparseTriangle(ap["g"]);
843 std::cout << "Make it sparser\n";
846 else if(ap["tree"]) {
849 else if(ap["dela"]) {
854 std::cout << "Number of nodes : " << countNodes(g) << std::endl;
855 std::cout << "Number of edges : " << countUEdges(g) << std::endl;
857 for(UEdgeIt e(g);e!=INVALID;++e)
858 tlen+=sqrt((coords[g.source(e)]-coords[g.target(e)]).normSquare());
859 std::cout << "Total edge length : " << tlen << std::endl;
862 graphToEps(g,prefix+".eps").scaleToA4().
863 scale(600).nodeScale(.2).edgeWidthScale(.001).preScale(false).
864 coords(coords).run();
867 GraphWriter<ListUGraph>(prefix+".lgf",g).
868 writeNodeMap("coordinates_x",scaleMap(xMap(coords),600)).
869 writeNodeMap("coordinates_y",scaleMap(yMap(coords),600)).
871 else UGraphWriter<ListUGraph>(prefix+".lgf",g).
872 writeNodeMap("coordinates_x",scaleMap(xMap(coords),600)).
873 writeNodeMap("coordinates_y",scaleMap(yMap(coords),600)).