Add missing include header and std:: namespace spec. (#487)
1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
3 * This file is a part of LEMON, a generic C++ optimization library.
5 * Copyright (C) 2003-2009
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
21 /// \brief Special plane graph generator.
23 /// Graph generator application for various types of plane graphs.
29 /// for more information on the usage.
34 #include <lemon/list_graph.h>
35 #include <lemon/random.h>
36 #include <lemon/dim2.h>
37 #include <lemon/bfs.h>
38 #include <lemon/counter.h>
39 #include <lemon/suurballe.h>
40 #include <lemon/graph_to_eps.h>
41 #include <lemon/lgf_writer.h>
42 #include <lemon/arg_parser.h>
43 #include <lemon/euler.h>
44 #include <lemon/math.h>
45 #include <lemon/kruskal.h>
46 #include <lemon/time_measure.h>
48 using namespace lemon;
50 typedef dim2::Point<double> Point;
52 GRAPH_TYPEDEFS(ListGraph);
61 std::vector<Node> nodes;
62 ListGraph::NodeMap<Point> coords(g);
67 for(EdgeIt e(g);e!=INVALID;++e)
68 tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
74 const double EPSILON=1e-8;
75 bool tsp_improve(Node u, Node v)
77 double luv=std::sqrt((coords[v]-coords[u]).normSquare());
82 for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e) { }
85 if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON>
86 std::sqrt((coords[u]-coords[u2]).normSquare())+
87 std::sqrt((coords[v]-coords[v2]).normSquare()))
89 g.erase(findEdge(g,u,v));
90 g.erase(findEdge(g,u2,v2));
100 bool tsp_improve(Node u)
102 for(IncEdgeIt e(g,u);e!=INVALID;++e)
103 if(tsp_improve(u,g.runningNode(e))) return true;
112 for(NodeIt n(g);n!=INVALID;++n)
113 if(tsp_improve(n)) b=true;
119 for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);
128 Line(Point _a,Point _b) :a(_a),b(_b) {}
129 Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}
130 Line(const Arc &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
131 Line(const Edge &e) : a(coords[g.u(e)]),b(coords[g.v(e)]) {}
134 inline std::ostream& operator<<(std::ostream &os, const Line &l)
136 os << l.a << "->" << l.b;
140 bool cross(Line a, Line b)
142 Point ao=rot90(a.b-a.a);
143 Point bo=rot90(b.b-b.a);
144 return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 &&
145 (bo*(a.a-b.a))*(bo*(a.b-b.a))<0;
155 bool pedgeLess(Parc a,Parc b)
160 std::vector<Edge> arcs;
162 namespace _delaunay_bits {
165 int prev, curr, next;
167 Part(int p, int c, int n) : prev(p), curr(c), next(n) {}
170 inline std::ostream& operator<<(std::ostream& os, const Part& part) {
171 os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';
175 inline double circle_point(const Point& p, const Point& q, const Point& r) {
176 double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);
177 if (a == 0) return std::numeric_limits<double>::quiet_NaN();
179 double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +
180 (q.x * q.x + q.y * q.y) * (r.y - p.y) +
181 (r.x * r.x + r.y * r.y) * (p.y - q.y);
183 double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +
184 (q.x * q.x + q.y * q.y) * (r.x - p.x) +
185 (r.x * r.x + r.y * r.y) * (p.x - q.x);
187 double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +
188 (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +
189 (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);
191 return d / (2 * a) + std::sqrt((d * d + e * e) / (4 * a * a) + f / a);
194 inline bool circle_form(const Point& p, const Point& q, const Point& r) {
195 return rot90(q - p) * (r - q) < 0.0;
198 inline double intersection(const Point& p, const Point& q, double sx) {
199 const double epsilon = 1e-8;
201 if (p.x == q.x) return (p.y + q.y) / 2.0;
203 if (sx < p.x + epsilon) return p.y;
204 if (sx < q.x + epsilon) return q.y;
206 double a = q.x - p.x;
207 double b = (q.x - sx) * p.y - (p.x - sx) * q.y;
208 double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();
209 return (b - std::sqrt(d)) / a;
215 YLess(const std::vector<Point>& points, double& sweep)
216 : _points(points), _sweep(sweep) {}
218 bool operator()(const Part& l, const Part& r) const {
219 const double epsilon = 1e-8;
221 // std::cerr << l << " vs " << r << std::endl;
222 double lbx = l.prev != -1 ?
223 intersection(_points[l.prev], _points[l.curr], _sweep) :
224 - std::numeric_limits<double>::infinity();
225 double rbx = r.prev != -1 ?
226 intersection(_points[r.prev], _points[r.curr], _sweep) :
227 - std::numeric_limits<double>::infinity();
228 double lex = l.next != -1 ?
229 intersection(_points[l.curr], _points[l.next], _sweep) :
230 std::numeric_limits<double>::infinity();
231 double rex = r.next != -1 ?
232 intersection(_points[r.curr], _points[r.next], _sweep) :
233 std::numeric_limits<double>::infinity();
235 if (lbx > lex) std::swap(lbx, lex);
236 if (rbx > rex) std::swap(rbx, rex);
238 if (lex < epsilon + rex && lbx + epsilon < rex) return true;
239 if (rex < epsilon + lex && rbx + epsilon < lex) return false;
243 const std::vector<Point>& _points;
249 typedef std::multimap<double, BeachIt*> SpikeHeap;
251 typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;
256 BeachIt(Beach::iterator iter) : it(iter) {}
261 inline void delaunay() {
262 Counter cnt("Number of arcs added: ");
264 using namespace _delaunay_bits;
266 typedef _delaunay_bits::Part Part;
267 typedef std::vector<std::pair<double, int> > SiteHeap;
270 std::vector<Point> points;
271 std::vector<Node> nodes;
273 for (NodeIt it(g); it != INVALID; ++it) {
275 points.push_back(coords[it]);
278 SiteHeap siteheap(points.size());
283 for (int i = 0; i < int(siteheap.size()); ++i) {
284 siteheap[i] = std::make_pair(points[i].x, i);
287 std::sort(siteheap.begin(), siteheap.end());
288 sweep = siteheap.front().first;
290 YLess yless(points, sweep);
295 std::set<std::pair<int, int> > arcs;
301 while (siteindex < int(siteheap.size()) &&
302 siteheap[0].first == siteheap[siteindex].first) {
303 front.push_back(std::make_pair(points[siteheap[siteindex].second].y,
304 siteheap[siteindex].second));
308 std::sort(front.begin(), front.end());
310 for (int i = 0; i < int(front.size()); ++i) {
311 int prev = (i == 0 ? -1 : front[i - 1].second);
312 int curr = front[i].second;
313 int next = (i + 1 == int(front.size()) ? -1 : front[i + 1].second);
315 beach.insert(std::make_pair(Part(prev, curr, next),
320 while (siteindex < int(points.size()) || !spikeheap.empty()) {
322 SpikeHeap::iterator spit = spikeheap.begin();
324 if (siteindex < int(points.size()) &&
325 (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {
326 int site = siteheap[siteindex].second;
327 sweep = siteheap[siteindex].first;
329 Beach::iterator bit = beach.upper_bound(Part(site, site, site));
331 if (bit->second != spikeheap.end()) {
332 delete bit->second->second;
333 spikeheap.erase(bit->second);
336 int prev = bit->first.prev;
337 int curr = bit->first.curr;
338 int next = bit->first.next;
342 SpikeHeap::iterator pit = spikeheap.end();
344 circle_form(points[prev], points[curr], points[site])) {
345 double x = circle_point(points[prev], points[curr], points[site]);
346 pit = spikeheap.insert(std::make_pair(x, new BeachIt(beach.end())));
348 beach.insert(std::make_pair(Part(prev, curr, site), pit));
350 beach.insert(std::make_pair(Part(prev, curr, site), pit));
353 beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
355 SpikeHeap::iterator nit = spikeheap.end();
357 circle_form(points[site], points[curr],points[next])) {
358 double x = circle_point(points[site], points[curr], points[next]);
359 nit = spikeheap.insert(std::make_pair(x, new BeachIt(beach.end())));
361 beach.insert(std::make_pair(Part(site, curr, next), nit));
363 beach.insert(std::make_pair(Part(site, curr, next), nit));
370 Beach::iterator bit = spit->second->it;
372 int prev = bit->first.prev;
373 int curr = bit->first.curr;
374 int next = bit->first.next;
377 std::pair<int, int> arc;
380 std::make_pair(prev, curr) : std::make_pair(curr, prev);
382 if (arcs.find(arc) == arcs.end()) {
384 g.addEdge(nodes[prev], nodes[curr]);
389 std::make_pair(curr, next) : std::make_pair(next, curr);
391 if (arcs.find(arc) == arcs.end()) {
393 g.addEdge(nodes[curr], nodes[next]);
398 Beach::iterator pbit = bit; --pbit;
399 int ppv = pbit->first.prev;
400 Beach::iterator nbit = bit; ++nbit;
401 int nnt = nbit->first.next;
403 if (bit->second != spikeheap.end())
405 delete bit->second->second;
406 spikeheap.erase(bit->second);
408 if (pbit->second != spikeheap.end())
410 delete pbit->second->second;
411 spikeheap.erase(pbit->second);
413 if (nbit->second != spikeheap.end())
415 delete nbit->second->second;
416 spikeheap.erase(nbit->second);
423 SpikeHeap::iterator pit = spikeheap.end();
424 if (ppv != -1 && ppv != next &&
425 circle_form(points[ppv], points[prev], points[next])) {
426 double x = circle_point(points[ppv], points[prev], points[next]);
427 if (x < sweep) x = sweep;
428 pit = spikeheap.insert(std::make_pair(x, new BeachIt(beach.end())));
430 beach.insert(std::make_pair(Part(ppv, prev, next), pit));
432 beach.insert(std::make_pair(Part(ppv, prev, next), pit));
435 SpikeHeap::iterator nit = spikeheap.end();
436 if (nnt != -1 && prev != nnt &&
437 circle_form(points[prev], points[next], points[nnt])) {
438 double x = circle_point(points[prev], points[next], points[nnt]);
439 if (x < sweep) x = sweep;
440 nit = spikeheap.insert(std::make_pair(x, new BeachIt(beach.end())));
442 beach.insert(std::make_pair(Part(prev, next, nnt), nit));
444 beach.insert(std::make_pair(Part(prev, next, nnt), nit));
450 for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
451 int curr = it->first.curr;
452 int next = it->first.next;
454 if (next == -1) continue;
456 std::pair<int, int> arc;
459 std::make_pair(curr, next) : std::make_pair(next, curr);
461 if (arcs.find(arc) == arcs.end()) {
463 g.addEdge(nodes[curr], nodes[next]);
471 Counter cnt("Number of arcs removed: ");
472 Bfs<ListGraph> bfs(g);
473 for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
474 ei!=arcs.rend();++ei)
480 if(bfs.predArc(b)==INVALID || bfs.dist(b)>d)
488 Counter cnt("Number of arcs removed: ");
489 for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
490 ei!=arcs.rend();++ei)
495 ConstMap<Arc,int> cegy(1);
496 Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy);
497 int k=sur.run(a,b,2);
498 if(k<2 || sur.totalLength()>d)
501 // else std::cout << "Remove arc " << g.id(a) << "-" << g.id(b) << '\n';
505 void sparseTriangle(int d)
507 Counter cnt("Number of arcs added: ");
508 std::vector<Parc> pedges;
509 for(NodeIt n(g);n!=INVALID;++n)
510 for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
515 p.len=(coords[m]-coords[n]).normSquare();
518 std::sort(pedges.begin(),pedges.end(),pedgeLess);
519 for(std::vector<Parc>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
521 Line li(pi->a,pi->b);
523 for(;e!=INVALID && !cross(e,li);++e) ;
526 ConstMap<Arc,int> cegy(1);
527 Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy);
528 int k=sur.run(pi->a,pi->b,2);
529 if(k<2 || sur.totalLength()>d)
531 ne=g.addEdge(pi->a,pi->b);
539 template <typename Graph, typename CoordMap>
540 class LengthSquareMap {
542 typedef typename Graph::Edge Key;
543 typedef typename CoordMap::Value::Value Value;
545 LengthSquareMap(const Graph& graph, const CoordMap& coords)
546 : _graph(graph), _coords(coords) {}
548 Value operator[](const Key& key) const {
549 return (_coords[_graph.v(key)] -
550 _coords[_graph.u(key)]).normSquare();
556 const CoordMap& _coords;
560 std::vector<Parc> pedges;
562 std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
564 std::cout << T.realTime() << "s: Calculating spanning tree...\n";
565 LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords);
566 ListGraph::EdgeMap<bool> tree(g);
567 kruskal(g, ls, tree);
568 std::cout << T.realTime() << "s: Removing non tree arcs...\n";
569 std::vector<Edge> remove;
570 for (EdgeIt e(g); e != INVALID; ++e) {
571 if (!tree[e]) remove.push_back(e);
573 for(int i = 0; i < int(remove.size()); ++i) {
576 std::cout << T.realTime() << "s: Done\n";
581 std::cout << "Find a tree..." << std::endl;
585 std::cout << "Total arc length (tree) : " << totalLen() << std::endl;
587 std::cout << "Make it Euler..." << std::endl;
590 std::vector<Node> leafs;
591 for(NodeIt n(g);n!=INVALID;++n)
592 if(countIncEdges(g,n)%2==1) leafs.push_back(n);
594 // for(unsigned int i=0;i<leafs.size();i+=2)
595 // g.addArc(leafs[i],leafs[i+1]);
597 std::vector<Parc> pedges;
598 for(unsigned int i=0;i<leafs.size()-1;i++)
599 for(unsigned int j=i+1;j<leafs.size();j++)
606 p.len=(coords[m]-coords[n]).normSquare();
609 std::sort(pedges.begin(),pedges.end(),pedgeLess);
610 for(unsigned int i=0;i<pedges.size();i++)
611 if(countIncEdges(g,pedges[i].a)%2 &&
612 countIncEdges(g,pedges[i].b)%2)
613 g.addEdge(pedges[i].a,pedges[i].b);
616 for(NodeIt n(g);n!=INVALID;++n)
617 if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
618 std::cout << "GEBASZ!!!" << std::endl;
620 for(EdgeIt e(g);e!=INVALID;++e)
622 std::cout << "LOOP GEBASZ!!!" << std::endl;
624 std::cout << "Number of arcs : " << countEdges(g) << std::endl;
626 std::cout << "Total arc length (euler) : " << totalLen() << std::endl;
628 ListGraph::EdgeMap<Arc> enext(g);
630 EulerIt<ListGraph> e(g);
633 // std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
634 for(++e;e!=INVALID;++e)
636 // std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
643 std::cout << "Creating a tour from that..." << std::endl;
645 int nnum = countNodes(g);
646 int ednum = countEdges(g);
648 for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p])
650 // std::cout << "Checking arc " << g.id(p) << std::endl;
654 Node n1=g.oppositeNode(n2,e);
655 Node n3=g.oppositeNode(n2,f);
656 if(countIncEdges(g,n2)>2)
658 // std::cout << "Remove an Arc" << std::endl;
664 Arc ne=g.direct(g.addEdge(n1,n3),n1);
676 std::cout << "Total arc length (tour) : " << totalLen() << std::endl;
678 std::cout << "2-opt the tour..." << std::endl;
682 std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl;
686 int main(int argc,const char **argv)
688 ArgParser ap(argc,argv);
691 bool disc_d, square_d, gauss_d;
692 // bool tsp_a,two_a,tree_a;
697 std::string ndist("disc");
698 ap.refOption("n", "Number of nodes (default is 100)", N)
699 .intOption("g", "Girth parameter (default is 10)", 10)
700 .refOption("cities", "Number of cities (default is 1)", num_of_cities)
701 .refOption("area", "Full relative area of the cities (default is 1)", area)
702 .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",
704 .optionGroup("dist", "disc")
705 .refOption("square", "Nodes are evenly distributed on a unit square",
707 .optionGroup("dist", "square")
708 .refOption("gauss", "Nodes are located according to a two-dim Gauss "
709 "distribution", gauss_d)
710 .optionGroup("dist", "gauss")
711 .onlyOneGroup("dist")
712 .boolOption("eps", "Also generate .eps output (<prefix>.eps)")
713 .boolOption("nonodes", "Draw only the edges in the generated .eps output")
714 .boolOption("dir", "Directed graph is generated (each edge is replaced by "
715 "two directed arcs)")
716 .boolOption("2con", "Create a two connected planar graph")
717 .optionGroup("alg","2con")
718 .boolOption("tree", "Create a min. cost spanning tree")
719 .optionGroup("alg","tree")
720 .boolOption("tsp", "Create a TSP tour")
721 .optionGroup("alg","tsp")
722 .boolOption("tsp2", "Create a TSP tour (tree based)")
723 .optionGroup("alg","tsp2")
724 .boolOption("dela", "Delaunay triangulation graph")
725 .optionGroup("alg","dela")
727 .boolOption("rand", "Use time seed for random number generator")
728 .optionGroup("rand", "rand")
729 .intOption("seed", "Random seed", -1)
730 .optionGroup("rand", "seed")
731 .onlyOneGroup("rand")
732 .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
736 int seed = int(time(0));
737 std::cout << "Random number seed: " << seed << std::endl;
740 if (ap.given("seed")) {
741 int seed = ap["seed"];
742 std::cout << "Random number seed: " << seed << std::endl;
747 switch(ap.files().size())
750 prefix="lgf-gen-out";
753 prefix=ap.files()[0];
756 std::cerr << "\nAt most one prefix can be given\n\n";
761 std::vector<double> sizes;
762 std::vector<double> cum_sizes;
763 for(int s=0;s<num_of_cities;s++)
765 // sum_sizes+=rnd.exponential();
769 cum_sizes.push_back(sum_sizes);
772 for(int s=0;s<num_of_cities;s++)
774 Point center=(num_of_cities==1?Point(0,0):rnd.disc());
776 for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
779 coords[n]=center+rnd.gauss2()*area*
780 std::sqrt(sizes[s]/sum_sizes);
783 for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
786 coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
787 std::sqrt(sizes[s]/sum_sizes);
789 else if(disc_d || true)
790 for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
793 coords[n]=center+rnd.disc()*area*
794 std::sqrt(sizes[s]/sum_sizes);
798 // for (ListGraph::NodeIt n(g); n != INVALID; ++n) {
799 // std::cerr << coords[n] << std::endl;
804 std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
808 std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
810 else if(ap["2con"]) {
811 std::cout << "Make triangles\n";
813 sparseTriangle(ap["g"]);
814 std::cout << "Make it sparser\n";
817 else if(ap["tree"]) {
820 else if(ap["dela"]) {
825 std::cout << "Number of nodes : " << countNodes(g) << std::endl;
826 std::cout << "Number of arcs : " << countEdges(g) << std::endl;
828 for(EdgeIt e(g);e!=INVALID;++e)
829 tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
830 std::cout << "Total arc length : " << tlen << std::endl;
833 graphToEps(g,prefix+".eps").scaleToA4().
834 scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false).
835 coords(coords).hideNodes(ap.given("nonodes")).run();
838 DigraphWriter<ListGraph>(g,prefix+".lgf").
839 nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
840 nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
842 else GraphWriter<ListGraph>(g,prefix+".lgf").
843 nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
844 nodeMap("coordinates_y",scaleMap(yMap(coords),600)).