1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
5 * Copyright (C) 2003-2009
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 */
18
19 /// \ingroup tools
20 /// \file
21 /// \brief Special plane graph generator.
22 ///
23 /// Graph generator application for various types of plane graphs.
24 ///
25 /// See
26 /// \code
27 /// lgf-gen --help
28 /// \endcode
29 /// for more information on the usage.
30
31 #include <algorithm>
32 #include <set>
33 #include <ctime>
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>
47
48 using namespace lemon;
49
50 typedef dim2::Point<double> Point;
51
52 GRAPH_TYPEDEFS(ListGraph);
53
54 bool progress=true;
55
56 int N;
57 // int girth;
58
59 ListGraph g;
60
61 std::vector<Node> nodes;
62 ListGraph::NodeMap<Point> coords(g);
63
64
totalLen()65 double totalLen(){
66 double tlen=0;
67 for(EdgeIt e(g);e!=INVALID;++e)
68 tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare());
69 return tlen;
70 }
71
72 int tsp_impr_num=0;
73
74 const double EPSILON=1e-8;
tsp_improve(Node u,Node v)75 bool tsp_improve(Node u, Node v)
76 {
77 double luv=std::sqrt((coords[v]-coords[u]).normSquare());
78 Node u2=u;
79 Node v2=v;
80 do {
81 Node n;
82 for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e) { }
83 u2=v2;
84 v2=n;
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()))
88 {
89 g.erase(findEdge(g,u,v));
90 g.erase(findEdge(g,u2,v2));
91 g.addEdge(u2,u);
92 g.addEdge(v,v2);
93 tsp_impr_num++;
94 return true;
95 }
96 } while(v2!=u);
97 return false;
98 }
99
tsp_improve(Node u)100 bool tsp_improve(Node u)
101 {
102 for(IncEdgeIt e(g,u);e!=INVALID;++e)
103 if(tsp_improve(u,g.runningNode(e))) return true;
104 return false;
105 }
106
tsp_improve()107 void tsp_improve()
108 {
109 bool b;
110 do {
111 b=false;
112 for(NodeIt n(g);n!=INVALID;++n)
113 if(tsp_improve(n)) b=true;
114 } while(b);
115 }
116
tsp()117 void tsp()
118 {
119 for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]);
120 tsp_improve();
121 }
122
123 class Line
124 {
125 public:
126 Point a;
127 Point b;
Line(Point _a,Point _b)128 Line(Point _a,Point _b) :a(_a),b(_b) {}
Line(Node _a,Node _b)129 Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {}
Line(const Arc & e)130 Line(const Arc &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {}
Line(const Edge & e)131 Line(const Edge &e) : a(coords[g.u(e)]),b(coords[g.v(e)]) {}
132 };
133
operator <<(std::ostream & os,const Line & l)134 inline std::ostream& operator<<(std::ostream &os, const Line &l)
135 {
136 os << l.a << "->" << l.b;
137 return os;
138 }
139
cross(Line a,Line b)140 bool cross(Line a, Line b)
141 {
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;
146 }
147
148 struct Parc
149 {
150 Node a;
151 Node b;
152 double len;
153 };
154
pedgeLess(Parc a,Parc b)155 bool pedgeLess(Parc a,Parc b)
156 {
157 return a.len<b.len;
158 }
159
160 std::vector<Edge> arcs;
161
162 namespace _delaunay_bits {
163
164 struct Part {
165 int prev, curr, next;
166
Part_delaunay_bits::Part167 Part(int p, int c, int n) : prev(p), curr(c), next(n) {}
168 };
169
operator <<(std::ostream & os,const Part & part)170 inline std::ostream& operator<<(std::ostream& os, const Part& part) {
171 os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';
172 return os;
173 }
174
circle_point(const Point & p,const Point & q,const Point & r)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();
178
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);
182
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);
186
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);
190
191 return d / (2 * a) + std::sqrt((d * d + e * e) / (4 * a * a) + f / a);
192 }
193
circle_form(const Point & p,const Point & q,const Point & r)194 inline bool circle_form(const Point& p, const Point& q, const Point& r) {
195 return rot90(q - p) * (r - q) < 0.0;
196 }
197
intersection(const Point & p,const Point & q,double sx)198 inline double intersection(const Point& p, const Point& q, double sx) {
199 const double epsilon = 1e-8;
200
201 if (p.x == q.x) return (p.y + q.y) / 2.0;
202
203 if (sx < p.x + epsilon) return p.y;
204 if (sx < q.x + epsilon) return q.y;
205
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;
210 }
211
212 struct YLess {
213
214
YLess_delaunay_bits::YLess215 YLess(const std::vector<Point>& points, double& sweep)
216 : _points(points), _sweep(sweep) {}
217
operator ()_delaunay_bits::YLess218 bool operator()(const Part& l, const Part& r) const {
219 const double epsilon = 1e-8;
220
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();
234
235 if (lbx > lex) std::swap(lbx, lex);
236 if (rbx > rex) std::swap(rbx, rex);
237
238 if (lex < epsilon + rex && lbx + epsilon < rex) return true;
239 if (rex < epsilon + lex && rbx + epsilon < lex) return false;
240 return lex < rex;
241 }
242
243 const std::vector<Point>& _points;
244 double& _sweep;
245 };
246
247 struct BeachIt;
248
249 typedef std::multimap<double, BeachIt*> SpikeHeap;
250
251 typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;
252
253 struct BeachIt {
254 Beach::iterator it;
255
BeachIt_delaunay_bits::BeachIt256 BeachIt(Beach::iterator iter) : it(iter) {}
257 };
258
259 }
260
delaunay()261 inline void delaunay() {
262 Counter cnt("Number of arcs added: ");
263
264 using namespace _delaunay_bits;
265
266 typedef _delaunay_bits::Part Part;
267 typedef std::vector<std::pair<double, int> > SiteHeap;
268
269
270 std::vector<Point> points;
271 std::vector<Node> nodes;
272
273 for (NodeIt it(g); it != INVALID; ++it) {
274 nodes.push_back(it);
275 points.push_back(coords[it]);
276 }
277
278 SiteHeap siteheap(points.size());
279
280 double sweep;
281
282
283 for (int i = 0; i < int(siteheap.size()); ++i) {
284 siteheap[i] = std::make_pair(points[i].x, i);
285 }
286
287 std::sort(siteheap.begin(), siteheap.end());
288 sweep = siteheap.front().first;
289
290 YLess yless(points, sweep);
291 Beach beach(yless);
292
293 SpikeHeap spikeheap;
294
295 std::set<std::pair<int, int> > arcs;
296
297 int siteindex = 0;
298 {
299 SiteHeap front;
300
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));
305 ++siteindex;
306 }
307
308 std::sort(front.begin(), front.end());
309
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);
314
315 beach.insert(std::make_pair(Part(prev, curr, next),
316 spikeheap.end()));
317 }
318 }
319
320 while (siteindex < int(points.size()) || !spikeheap.empty()) {
321
322 SpikeHeap::iterator spit = spikeheap.begin();
323
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;
328
329 Beach::iterator bit = beach.upper_bound(Part(site, site, site));
330
331 if (bit->second != spikeheap.end()) {
332 delete bit->second->second;
333 spikeheap.erase(bit->second);
334 }
335
336 int prev = bit->first.prev;
337 int curr = bit->first.curr;
338 int next = bit->first.next;
339
340 beach.erase(bit);
341
342 SpikeHeap::iterator pit = spikeheap.end();
343 if (prev != -1 &&
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())));
347 pit->second->it =
348 beach.insert(std::make_pair(Part(prev, curr, site), pit));
349 } else {
350 beach.insert(std::make_pair(Part(prev, curr, site), pit));
351 }
352
353 beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
354
355 SpikeHeap::iterator nit = spikeheap.end();
356 if (next != -1 &&
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())));
360 nit->second->it =
361 beach.insert(std::make_pair(Part(site, curr, next), nit));
362 } else {
363 beach.insert(std::make_pair(Part(site, curr, next), nit));
364 }
365
366 ++siteindex;
367 } else {
368 sweep = spit->first;
369
370 Beach::iterator bit = spit->second->it;
371
372 int prev = bit->first.prev;
373 int curr = bit->first.curr;
374 int next = bit->first.next;
375
376 {
377 std::pair<int, int> arc;
378
379 arc = prev < curr ?
380 std::make_pair(prev, curr) : std::make_pair(curr, prev);
381
382 if (arcs.find(arc) == arcs.end()) {
383 arcs.insert(arc);
384 g.addEdge(nodes[prev], nodes[curr]);
385 ++cnt;
386 }
387
388 arc = curr < next ?
389 std::make_pair(curr, next) : std::make_pair(next, curr);
390
391 if (arcs.find(arc) == arcs.end()) {
392 arcs.insert(arc);
393 g.addEdge(nodes[curr], nodes[next]);
394 ++cnt;
395 }
396 }
397
398 Beach::iterator pbit = bit; --pbit;
399 int ppv = pbit->first.prev;
400 Beach::iterator nbit = bit; ++nbit;
401 int nnt = nbit->first.next;
402
403 if (bit->second != spikeheap.end())
404 {
405 delete bit->second->second;
406 spikeheap.erase(bit->second);
407 }
408 if (pbit->second != spikeheap.end())
409 {
410 delete pbit->second->second;
411 spikeheap.erase(pbit->second);
412 }
413 if (nbit->second != spikeheap.end())
414 {
415 delete nbit->second->second;
416 spikeheap.erase(nbit->second);
417 }
418
419 beach.erase(nbit);
420 beach.erase(bit);
421 beach.erase(pbit);
422
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())));
429 pit->second->it =
430 beach.insert(std::make_pair(Part(ppv, prev, next), pit));
431 } else {
432 beach.insert(std::make_pair(Part(ppv, prev, next), pit));
433 }
434
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())));
441 nit->second->it =
442 beach.insert(std::make_pair(Part(prev, next, nnt), nit));
443 } else {
444 beach.insert(std::make_pair(Part(prev, next, nnt), nit));
445 }
446
447 }
448 }
449
450 for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
451 int curr = it->first.curr;
452 int next = it->first.next;
453
454 if (next == -1) continue;
455
456 std::pair<int, int> arc;
457
458 arc = curr < next ?
459 std::make_pair(curr, next) : std::make_pair(next, curr);
460
461 if (arcs.find(arc) == arcs.end()) {
462 arcs.insert(arc);
463 g.addEdge(nodes[curr], nodes[next]);
464 ++cnt;
465 }
466 }
467 }
468
sparse(int d)469 void sparse(int d)
470 {
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)
475 {
476 Node a=g.u(*ei);
477 Node b=g.v(*ei);
478 g.erase(*ei);
479 bfs.run(a,b);
480 if(bfs.predArc(b)==INVALID || bfs.dist(b)>d)
481 g.addEdge(a,b);
482 else cnt++;
483 }
484 }
485
sparse2(int d)486 void sparse2(int d)
487 {
488 Counter cnt("Number of arcs removed: ");
489 for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin();
490 ei!=arcs.rend();++ei)
491 {
492 Node a=g.u(*ei);
493 Node b=g.v(*ei);
494 g.erase(*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)
499 g.addEdge(a,b);
500 else cnt++;
501 // else std::cout << "Remove arc " << g.id(a) << "-" << g.id(b) << '\n';
502 }
503 }
504
sparseTriangle(int d)505 void sparseTriangle(int d)
506 {
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)
511 {
512 Parc p;
513 p.a=n;
514 p.b=m;
515 p.len=(coords[m]-coords[n]).normSquare();
516 pedges.push_back(p);
517 }
518 std::sort(pedges.begin(),pedges.end(),pedgeLess);
519 for(std::vector<Parc>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
520 {
521 Line li(pi->a,pi->b);
522 EdgeIt e(g);
523 for(;e!=INVALID && !cross(e,li);++e) ;
524 Edge ne;
525 if(e==INVALID) {
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)
530 {
531 ne=g.addEdge(pi->a,pi->b);
532 arcs.push_back(ne);
533 cnt++;
534 }
535 }
536 }
537 }
538
539 template <typename Graph, typename CoordMap>
540 class LengthSquareMap {
541 public:
542 typedef typename Graph::Edge Key;
543 typedef typename CoordMap::Value::Value Value;
544
LengthSquareMap(const Graph & graph,const CoordMap & coords)545 LengthSquareMap(const Graph& graph, const CoordMap& coords)
546 : _graph(graph), _coords(coords) {}
547
operator [](const Key & key) const548 Value operator[](const Key& key) const {
549 return (_coords[_graph.v(key)] -
550 _coords[_graph.u(key)]).normSquare();
551 }
552
553 private:
554
555 const Graph& _graph;
556 const CoordMap& _coords;
557 };
558
minTree()559 void minTree() {
560 std::vector<Parc> pedges;
561 Timer T;
562 std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
563 delaunay();
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);
572 }
573 for(int i = 0; i < int(remove.size()); ++i) {
574 g.erase(remove[i]);
575 }
576 std::cout << T.realTime() << "s: Done\n";
577 }
578
tsp2()579 void tsp2()
580 {
581 std::cout << "Find a tree..." << std::endl;
582
583 minTree();
584
585 std::cout << "Total arc length (tree) : " << totalLen() << std::endl;
586
587 std::cout << "Make it Euler..." << std::endl;
588
589 {
590 std::vector<Node> leafs;
591 for(NodeIt n(g);n!=INVALID;++n)
592 if(countIncEdges(g,n)%2==1) leafs.push_back(n);
593
594 // for(unsigned int i=0;i<leafs.size();i+=2)
595 // g.addArc(leafs[i],leafs[i+1]);
596
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++)
600 {
601 Node n=leafs[i];
602 Node m=leafs[j];
603 Parc p;
604 p.a=n;
605 p.b=m;
606 p.len=(coords[m]-coords[n]).normSquare();
607 pedges.push_back(p);
608 }
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);
614 }
615
616 for(NodeIt n(g);n!=INVALID;++n)
617 if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 )
618 std::cout << "GEBASZ!!!" << std::endl;
619
620 for(EdgeIt e(g);e!=INVALID;++e)
621 if(g.u(e)==g.v(e))
622 std::cout << "LOOP GEBASZ!!!" << std::endl;
623
624 std::cout << "Number of arcs : " << countEdges(g) << std::endl;
625
626 std::cout << "Total arc length (euler) : " << totalLen() << std::endl;
627
628 ListGraph::EdgeMap<Arc> enext(g);
629 {
630 EulerIt<ListGraph> e(g);
631 Arc eo=e;
632 Arc ef=e;
633 // std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
634 for(++e;e!=INVALID;++e)
635 {
636 // std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl;
637 enext[eo]=e;
638 eo=e;
639 }
640 enext[eo]=ef;
641 }
642
643 std::cout << "Creating a tour from that..." << std::endl;
644
645 int nnum = countNodes(g);
646 int ednum = countEdges(g);
647
648 for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p])
649 {
650 // std::cout << "Checking arc " << g.id(p) << std::endl;
651 Arc e=enext[p];
652 Arc f=enext[e];
653 Node n2=g.source(f);
654 Node n1=g.oppositeNode(n2,e);
655 Node n3=g.oppositeNode(n2,f);
656 if(countIncEdges(g,n2)>2)
657 {
658 // std::cout << "Remove an Arc" << std::endl;
659 Arc ff=enext[f];
660 g.erase(e);
661 g.erase(f);
662 if(n1!=n3)
663 {
664 Arc ne=g.direct(g.addEdge(n1,n3),n1);
665 enext[p]=ne;
666 enext[ne]=ff;
667 ednum--;
668 }
669 else {
670 enext[p]=ff;
671 ednum-=2;
672 }
673 }
674 }
675
676 std::cout << "Total arc length (tour) : " << totalLen() << std::endl;
677
678 std::cout << "2-opt the tour..." << std::endl;
679
680 tsp_improve();
681
682 std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl;
683 }
684
685
main(int argc,const char ** argv)686 int main(int argc,const char **argv)
687 {
688 ArgParser ap(argc,argv);
689
690 // bool eps;
691 bool disc_d, square_d, gauss_d;
692 // bool tsp_a,two_a,tree_a;
693 int num_of_cities=1;
694 double area=1;
695 N=100;
696 // girth=10;
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)",
703 disc_d)
704 .optionGroup("dist", "disc")
705 .refOption("square", "Nodes are evenly distributed on a unit square",
706 square_d)
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")
726 .onlyOneGroup("alg")
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'")
733 .run();
734
735 if (ap["rand"]) {
736 int seed = int(time(0));
737 std::cout << "Random number seed: " << seed << std::endl;
738 rnd = Random(seed);
739 }
740 if (ap.given("seed")) {
741 int seed = ap["seed"];
742 std::cout << "Random number seed: " << seed << std::endl;
743 rnd = Random(seed);
744 }
745
746 std::string prefix;
747 switch(ap.files().size())
748 {
749 case 0:
750 prefix="lgf-gen-out";
751 break;
752 case 1:
753 prefix=ap.files()[0];
754 break;
755 default:
756 std::cerr << "\nAt most one prefix can be given\n\n";
757 exit(1);
758 }
759
760 double sum_sizes=0;
761 std::vector<double> sizes;
762 std::vector<double> cum_sizes;
763 for(int s=0;s<num_of_cities;s++)
764 {
765 // sum_sizes+=rnd.exponential();
766 double d=rnd();
767 sum_sizes+=d;
768 sizes.push_back(d);
769 cum_sizes.push_back(sum_sizes);
770 }
771 int i=0;
772 for(int s=0;s<num_of_cities;s++)
773 {
774 Point center=(num_of_cities==1?Point(0,0):rnd.disc());
775 if(gauss_d)
776 for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
777 Node n=g.addNode();
778 nodes.push_back(n);
779 coords[n]=center+rnd.gauss2()*area*
780 std::sqrt(sizes[s]/sum_sizes);
781 }
782 else if(square_d)
783 for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
784 Node n=g.addNode();
785 nodes.push_back(n);
786 coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area*
787 std::sqrt(sizes[s]/sum_sizes);
788 }
789 else if(disc_d || true)
790 for(;i<N*(cum_sizes[s]/sum_sizes);i++) {
791 Node n=g.addNode();
792 nodes.push_back(n);
793 coords[n]=center+rnd.disc()*area*
794 std::sqrt(sizes[s]/sum_sizes);
795 }
796 }
797
798 // for (ListGraph::NodeIt n(g); n != INVALID; ++n) {
799 // std::cerr << coords[n] << std::endl;
800 // }
801
802 if(ap["tsp"]) {
803 tsp();
804 std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
805 }
806 if(ap["tsp2"]) {
807 tsp2();
808 std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl;
809 }
810 else if(ap["2con"]) {
811 std::cout << "Make triangles\n";
812 // triangle();
813 sparseTriangle(ap["g"]);
814 std::cout << "Make it sparser\n";
815 sparse2(ap["g"]);
816 }
817 else if(ap["tree"]) {
818 minTree();
819 }
820 else if(ap["dela"]) {
821 delaunay();
822 }
823
824
825 std::cout << "Number of nodes : " << countNodes(g) << std::endl;
826 std::cout << "Number of arcs : " << countEdges(g) << std::endl;
827 double tlen=0;
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;
831
832 if(ap["eps"])
833 graphToEps(g,prefix+".eps").scaleToA4().
834 scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false).
835 coords(coords).hideNodes(ap.given("nonodes")).run();
836
837 if(ap["dir"])
838 DigraphWriter<ListGraph>(g,prefix+".lgf").
839 nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
840 nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
841 run();
842 else GraphWriter<ListGraph>(g,prefix+".lgf").
843 nodeMap("coordinates_x",scaleMap(xMap(coords),600)).
844 nodeMap("coordinates_y",scaleMap(yMap(coords),600)).
845 run();
846 }
847
848