|
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 digraph generator. |
|
22 /// |
|
23 /// Graph generator application for various types of plane graphs. |
|
24 /// |
|
25 /// See |
|
26 /// \verbatim |
|
27 /// lgf-gen --help |
|
28 /// \endverbatim |
|
29 /// for more info on the usage. |
|
30 /// |
|
31 |
|
32 |
|
33 #include <algorithm> |
|
34 #include <set> |
|
35 #include <lemon/list_graph.h> |
|
36 #include <lemon/random.h> |
|
37 #include <lemon/dim2.h> |
|
38 #include <lemon/bfs.h> |
|
39 #include <lemon/counter.h> |
|
40 #include <lemon/suurballe.h> |
|
41 #include <lemon/graph_to_eps.h> |
|
42 #include <lemon/lgf_writer.h> |
|
43 #include <lemon/arg_parser.h> |
|
44 #include <lemon/euler.h> |
|
45 #include <lemon/math.h> |
|
46 #include <lemon/kruskal.h> |
|
47 #include <lemon/time_measure.h> |
|
48 |
|
49 using namespace lemon; |
|
50 |
|
51 typedef dim2::Point<double> Point; |
|
52 |
|
53 GRAPH_TYPEDEFS(ListGraph); |
|
54 |
|
55 bool progress=true; |
|
56 |
|
57 int N; |
|
58 // int girth; |
|
59 |
|
60 ListGraph g; |
|
61 |
|
62 std::vector<Node> nodes; |
|
63 ListGraph::NodeMap<Point> coords(g); |
|
64 |
|
65 |
|
66 double totalLen(){ |
|
67 double tlen=0; |
|
68 for(EdgeIt e(g);e!=INVALID;++e) |
|
69 tlen+=sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare()); |
|
70 return tlen; |
|
71 } |
|
72 |
|
73 int tsp_impr_num=0; |
|
74 |
|
75 const double EPSILON=1e-8; |
|
76 bool tsp_improve(Node u, Node v) |
|
77 { |
|
78 double luv=std::sqrt((coords[v]-coords[u]).normSquare()); |
|
79 Node u2=u; |
|
80 Node v2=v; |
|
81 do { |
|
82 Node n; |
|
83 for(IncEdgeIt e(g,v2);(n=g.runningNode(e))==u2;++e) { } |
|
84 u2=v2; |
|
85 v2=n; |
|
86 if(luv+std::sqrt((coords[v2]-coords[u2]).normSquare())-EPSILON> |
|
87 std::sqrt((coords[u]-coords[u2]).normSquare())+ |
|
88 std::sqrt((coords[v]-coords[v2]).normSquare())) |
|
89 { |
|
90 g.erase(findEdge(g,u,v)); |
|
91 g.erase(findEdge(g,u2,v2)); |
|
92 g.addEdge(u2,u); |
|
93 g.addEdge(v,v2); |
|
94 tsp_impr_num++; |
|
95 return true; |
|
96 } |
|
97 } while(v2!=u); |
|
98 return false; |
|
99 } |
|
100 |
|
101 bool tsp_improve(Node u) |
|
102 { |
|
103 for(IncEdgeIt e(g,u);e!=INVALID;++e) |
|
104 if(tsp_improve(u,g.runningNode(e))) return true; |
|
105 return false; |
|
106 } |
|
107 |
|
108 void tsp_improve() |
|
109 { |
|
110 bool b; |
|
111 do { |
|
112 b=false; |
|
113 for(NodeIt n(g);n!=INVALID;++n) |
|
114 if(tsp_improve(n)) b=true; |
|
115 } while(b); |
|
116 } |
|
117 |
|
118 void tsp() |
|
119 { |
|
120 for(int i=0;i<N;i++) g.addEdge(nodes[i],nodes[(i+1)%N]); |
|
121 tsp_improve(); |
|
122 } |
|
123 |
|
124 class Line |
|
125 { |
|
126 public: |
|
127 Point a; |
|
128 Point b; |
|
129 Line(Point _a,Point _b) :a(_a),b(_b) {} |
|
130 Line(Node _a,Node _b) : a(coords[_a]),b(coords[_b]) {} |
|
131 Line(const Arc &e) : a(coords[g.source(e)]),b(coords[g.target(e)]) {} |
|
132 Line(const Edge &e) : a(coords[g.u(e)]),b(coords[g.v(e)]) {} |
|
133 }; |
|
134 |
|
135 inline std::ostream& operator<<(std::ostream &os, const Line &l) |
|
136 { |
|
137 os << l.a << "->" << l.b; |
|
138 return os; |
|
139 } |
|
140 |
|
141 bool cross(Line a, Line b) |
|
142 { |
|
143 Point ao=rot90(a.b-a.a); |
|
144 Point bo=rot90(b.b-b.a); |
|
145 return (ao*(b.a-a.a))*(ao*(b.b-a.a))<0 && |
|
146 (bo*(a.a-b.a))*(bo*(a.b-b.a))<0; |
|
147 } |
|
148 |
|
149 struct Parc |
|
150 { |
|
151 Node a; |
|
152 Node b; |
|
153 double len; |
|
154 }; |
|
155 |
|
156 bool pedgeLess(Parc a,Parc b) |
|
157 { |
|
158 return a.len<b.len; |
|
159 } |
|
160 |
|
161 std::vector<Edge> arcs; |
|
162 |
|
163 namespace _delaunay_bits { |
|
164 |
|
165 struct Part { |
|
166 int prev, curr, next; |
|
167 |
|
168 Part(int p, int c, int n) : prev(p), curr(c), next(n) {} |
|
169 }; |
|
170 |
|
171 inline std::ostream& operator<<(std::ostream& os, const Part& part) { |
|
172 os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')'; |
|
173 return os; |
|
174 } |
|
175 |
|
176 inline double circle_point(const Point& p, const Point& q, const Point& r) { |
|
177 double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y); |
|
178 if (a == 0) return std::numeric_limits<double>::quiet_NaN(); |
|
179 |
|
180 double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) + |
|
181 (q.x * q.x + q.y * q.y) * (r.y - p.y) + |
|
182 (r.x * r.x + r.y * r.y) * (p.y - q.y); |
|
183 |
|
184 double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) + |
|
185 (q.x * q.x + q.y * q.y) * (r.x - p.x) + |
|
186 (r.x * r.x + r.y * r.y) * (p.x - q.x); |
|
187 |
|
188 double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) + |
|
189 (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) + |
|
190 (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y); |
|
191 |
|
192 return d / (2 * a) + sqrt((d * d + e * e) / (4 * a * a) + f / a); |
|
193 } |
|
194 |
|
195 inline bool circle_form(const Point& p, const Point& q, const Point& r) { |
|
196 return rot90(q - p) * (r - q) < 0.0; |
|
197 } |
|
198 |
|
199 inline double intersection(const Point& p, const Point& q, double sx) { |
|
200 const double epsilon = 1e-8; |
|
201 |
|
202 if (p.x == q.x) return (p.y + q.y) / 2.0; |
|
203 |
|
204 if (sx < p.x + epsilon) return p.y; |
|
205 if (sx < q.x + epsilon) return q.y; |
|
206 |
|
207 double a = q.x - p.x; |
|
208 double b = (q.x - sx) * p.y - (p.x - sx) * q.y; |
|
209 double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare(); |
|
210 return (b - sqrt(d)) / a; |
|
211 } |
|
212 |
|
213 struct YLess { |
|
214 |
|
215 |
|
216 YLess(const std::vector<Point>& points, double& sweep) |
|
217 : _points(points), _sweep(sweep) {} |
|
218 |
|
219 bool operator()(const Part& l, const Part& r) const { |
|
220 const double epsilon = 1e-8; |
|
221 |
|
222 // std::cerr << l << " vs " << r << std::endl; |
|
223 double lbx = l.prev != -1 ? |
|
224 intersection(_points[l.prev], _points[l.curr], _sweep) : |
|
225 - std::numeric_limits<double>::infinity(); |
|
226 double rbx = r.prev != -1 ? |
|
227 intersection(_points[r.prev], _points[r.curr], _sweep) : |
|
228 - std::numeric_limits<double>::infinity(); |
|
229 double lex = l.next != -1 ? |
|
230 intersection(_points[l.curr], _points[l.next], _sweep) : |
|
231 std::numeric_limits<double>::infinity(); |
|
232 double rex = r.next != -1 ? |
|
233 intersection(_points[r.curr], _points[r.next], _sweep) : |
|
234 std::numeric_limits<double>::infinity(); |
|
235 |
|
236 if (lbx > lex) std::swap(lbx, lex); |
|
237 if (rbx > rex) std::swap(rbx, rex); |
|
238 |
|
239 if (lex < epsilon + rex && lbx + epsilon < rex) return true; |
|
240 if (rex < epsilon + lex && rbx + epsilon < lex) return false; |
|
241 return lex < rex; |
|
242 } |
|
243 |
|
244 const std::vector<Point>& _points; |
|
245 double& _sweep; |
|
246 }; |
|
247 |
|
248 struct BeachIt; |
|
249 |
|
250 typedef std::multimap<double, BeachIt> SpikeHeap; |
|
251 |
|
252 typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach; |
|
253 |
|
254 struct BeachIt { |
|
255 Beach::iterator it; |
|
256 |
|
257 BeachIt(Beach::iterator iter) : it(iter) {} |
|
258 }; |
|
259 |
|
260 } |
|
261 |
|
262 inline void delaunay() { |
|
263 Counter cnt("Number of arcs added: "); |
|
264 |
|
265 using namespace _delaunay_bits; |
|
266 |
|
267 typedef _delaunay_bits::Part Part; |
|
268 typedef std::vector<std::pair<double, int> > SiteHeap; |
|
269 |
|
270 |
|
271 std::vector<Point> points; |
|
272 std::vector<Node> nodes; |
|
273 |
|
274 for (NodeIt it(g); it != INVALID; ++it) { |
|
275 nodes.push_back(it); |
|
276 points.push_back(coords[it]); |
|
277 } |
|
278 |
|
279 SiteHeap siteheap(points.size()); |
|
280 |
|
281 double sweep; |
|
282 |
|
283 |
|
284 for (int i = 0; i < int(siteheap.size()); ++i) { |
|
285 siteheap[i] = std::make_pair(points[i].x, i); |
|
286 } |
|
287 |
|
288 std::sort(siteheap.begin(), siteheap.end()); |
|
289 sweep = siteheap.front().first; |
|
290 |
|
291 YLess yless(points, sweep); |
|
292 Beach beach(yless); |
|
293 |
|
294 SpikeHeap spikeheap; |
|
295 |
|
296 std::set<std::pair<int, int> > arcs; |
|
297 |
|
298 int siteindex = 0; |
|
299 { |
|
300 SiteHeap front; |
|
301 |
|
302 while (siteindex < int(siteheap.size()) && |
|
303 siteheap[0].first == siteheap[siteindex].first) { |
|
304 front.push_back(std::make_pair(points[siteheap[siteindex].second].y, |
|
305 siteheap[siteindex].second)); |
|
306 ++siteindex; |
|
307 } |
|
308 |
|
309 std::sort(front.begin(), front.end()); |
|
310 |
|
311 for (int i = 0; i < int(front.size()); ++i) { |
|
312 int prev = (i == 0 ? -1 : front[i - 1].second); |
|
313 int curr = front[i].second; |
|
314 int next = (i + 1 == int(front.size()) ? -1 : front[i + 1].second); |
|
315 |
|
316 beach.insert(std::make_pair(Part(prev, curr, next), |
|
317 spikeheap.end())); |
|
318 } |
|
319 } |
|
320 |
|
321 while (siteindex < int(points.size()) || !spikeheap.empty()) { |
|
322 |
|
323 SpikeHeap::iterator spit = spikeheap.begin(); |
|
324 |
|
325 if (siteindex < int(points.size()) && |
|
326 (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) { |
|
327 int site = siteheap[siteindex].second; |
|
328 sweep = siteheap[siteindex].first; |
|
329 |
|
330 Beach::iterator bit = beach.upper_bound(Part(site, site, site)); |
|
331 |
|
332 if (bit->second != spikeheap.end()) { |
|
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, 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, 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()) spikeheap.erase(bit->second); |
|
404 if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second); |
|
405 if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second); |
|
406 |
|
407 beach.erase(nbit); |
|
408 beach.erase(bit); |
|
409 beach.erase(pbit); |
|
410 |
|
411 SpikeHeap::iterator pit = spikeheap.end(); |
|
412 if (ppv != -1 && ppv != next && |
|
413 circle_form(points[ppv], points[prev], points[next])) { |
|
414 double x = circle_point(points[ppv], points[prev], points[next]); |
|
415 if (x < sweep) x = sweep; |
|
416 pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); |
|
417 pit->second.it = |
|
418 beach.insert(std::make_pair(Part(ppv, prev, next), pit)); |
|
419 } else { |
|
420 beach.insert(std::make_pair(Part(ppv, prev, next), pit)); |
|
421 } |
|
422 |
|
423 SpikeHeap::iterator nit = spikeheap.end(); |
|
424 if (nnt != -1 && prev != nnt && |
|
425 circle_form(points[prev], points[next], points[nnt])) { |
|
426 double x = circle_point(points[prev], points[next], points[nnt]); |
|
427 if (x < sweep) x = sweep; |
|
428 nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); |
|
429 nit->second.it = |
|
430 beach.insert(std::make_pair(Part(prev, next, nnt), nit)); |
|
431 } else { |
|
432 beach.insert(std::make_pair(Part(prev, next, nnt), nit)); |
|
433 } |
|
434 |
|
435 } |
|
436 } |
|
437 |
|
438 for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) { |
|
439 int curr = it->first.curr; |
|
440 int next = it->first.next; |
|
441 |
|
442 if (next == -1) continue; |
|
443 |
|
444 std::pair<int, int> arc; |
|
445 |
|
446 arc = curr < next ? |
|
447 std::make_pair(curr, next) : std::make_pair(next, curr); |
|
448 |
|
449 if (arcs.find(arc) == arcs.end()) { |
|
450 arcs.insert(arc); |
|
451 g.addEdge(nodes[curr], nodes[next]); |
|
452 ++cnt; |
|
453 } |
|
454 } |
|
455 } |
|
456 |
|
457 void sparse(int d) |
|
458 { |
|
459 Counter cnt("Number of arcs removed: "); |
|
460 Bfs<ListGraph> bfs(g); |
|
461 for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin(); |
|
462 ei!=arcs.rend();++ei) |
|
463 { |
|
464 Node a=g.u(*ei); |
|
465 Node b=g.v(*ei); |
|
466 g.erase(*ei); |
|
467 bfs.run(a,b); |
|
468 if(bfs.predArc(b)==INVALID || bfs.dist(b)>d) |
|
469 g.addEdge(a,b); |
|
470 else cnt++; |
|
471 } |
|
472 } |
|
473 |
|
474 void sparse2(int d) |
|
475 { |
|
476 Counter cnt("Number of arcs removed: "); |
|
477 for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin(); |
|
478 ei!=arcs.rend();++ei) |
|
479 { |
|
480 Node a=g.u(*ei); |
|
481 Node b=g.v(*ei); |
|
482 g.erase(*ei); |
|
483 ConstMap<Arc,int> cegy(1); |
|
484 Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy,a,b); |
|
485 int k=sur.run(2); |
|
486 if(k<2 || sur.totalLength()>d) |
|
487 g.addEdge(a,b); |
|
488 else cnt++; |
|
489 // else std::cout << "Remove arc " << g.id(a) << "-" << g.id(b) << '\n'; |
|
490 } |
|
491 } |
|
492 |
|
493 void sparseTriangle(int d) |
|
494 { |
|
495 Counter cnt("Number of arcs added: "); |
|
496 std::vector<Parc> pedges; |
|
497 for(NodeIt n(g);n!=INVALID;++n) |
|
498 for(NodeIt m=++(NodeIt(n));m!=INVALID;++m) |
|
499 { |
|
500 Parc p; |
|
501 p.a=n; |
|
502 p.b=m; |
|
503 p.len=(coords[m]-coords[n]).normSquare(); |
|
504 pedges.push_back(p); |
|
505 } |
|
506 std::sort(pedges.begin(),pedges.end(),pedgeLess); |
|
507 for(std::vector<Parc>::iterator pi=pedges.begin();pi!=pedges.end();++pi) |
|
508 { |
|
509 Line li(pi->a,pi->b); |
|
510 EdgeIt e(g); |
|
511 for(;e!=INVALID && !cross(e,li);++e) ; |
|
512 Edge ne; |
|
513 if(e==INVALID) { |
|
514 ConstMap<Arc,int> cegy(1); |
|
515 Suurballe<ListGraph,ConstMap<Arc,int> > |
|
516 sur(g,cegy,pi->a,pi->b); |
|
517 int k=sur.run(2); |
|
518 if(k<2 || sur.totalLength()>d) |
|
519 { |
|
520 ne=g.addEdge(pi->a,pi->b); |
|
521 arcs.push_back(ne); |
|
522 cnt++; |
|
523 } |
|
524 } |
|
525 } |
|
526 } |
|
527 |
|
528 template <typename Graph, typename CoordMap> |
|
529 class LengthSquareMap { |
|
530 public: |
|
531 typedef typename Graph::Edge Key; |
|
532 typedef typename CoordMap::Value::Value Value; |
|
533 |
|
534 LengthSquareMap(const Graph& graph, const CoordMap& coords) |
|
535 : _graph(graph), _coords(coords) {} |
|
536 |
|
537 Value operator[](const Key& key) const { |
|
538 return (_coords[_graph.v(key)] - |
|
539 _coords[_graph.u(key)]).normSquare(); |
|
540 } |
|
541 |
|
542 private: |
|
543 |
|
544 const Graph& _graph; |
|
545 const CoordMap& _coords; |
|
546 }; |
|
547 |
|
548 void minTree() { |
|
549 std::vector<Parc> pedges; |
|
550 Timer T; |
|
551 std::cout << T.realTime() << "s: Creating delaunay triangulation...\n"; |
|
552 delaunay(); |
|
553 std::cout << T.realTime() << "s: Calculating spanning tree...\n"; |
|
554 LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords); |
|
555 ListGraph::EdgeMap<bool> tree(g); |
|
556 kruskal(g, ls, tree); |
|
557 std::cout << T.realTime() << "s: Removing non tree arcs...\n"; |
|
558 std::vector<Edge> remove; |
|
559 for (EdgeIt e(g); e != INVALID; ++e) { |
|
560 if (!tree[e]) remove.push_back(e); |
|
561 } |
|
562 for(int i = 0; i < int(remove.size()); ++i) { |
|
563 g.erase(remove[i]); |
|
564 } |
|
565 std::cout << T.realTime() << "s: Done\n"; |
|
566 } |
|
567 |
|
568 void tsp2() |
|
569 { |
|
570 std::cout << "Find a tree..." << std::endl; |
|
571 |
|
572 minTree(); |
|
573 |
|
574 std::cout << "Total arc length (tree) : " << totalLen() << std::endl; |
|
575 |
|
576 std::cout << "Make it Euler..." << std::endl; |
|
577 |
|
578 { |
|
579 std::vector<Node> leafs; |
|
580 for(NodeIt n(g);n!=INVALID;++n) |
|
581 if(countIncEdges(g,n)%2==1) leafs.push_back(n); |
|
582 |
|
583 // for(unsigned int i=0;i<leafs.size();i+=2) |
|
584 // g.addArc(leafs[i],leafs[i+1]); |
|
585 |
|
586 std::vector<Parc> pedges; |
|
587 for(unsigned int i=0;i<leafs.size()-1;i++) |
|
588 for(unsigned int j=i+1;j<leafs.size();j++) |
|
589 { |
|
590 Node n=leafs[i]; |
|
591 Node m=leafs[j]; |
|
592 Parc p; |
|
593 p.a=n; |
|
594 p.b=m; |
|
595 p.len=(coords[m]-coords[n]).normSquare(); |
|
596 pedges.push_back(p); |
|
597 } |
|
598 std::sort(pedges.begin(),pedges.end(),pedgeLess); |
|
599 for(unsigned int i=0;i<pedges.size();i++) |
|
600 if(countIncEdges(g,pedges[i].a)%2 && |
|
601 countIncEdges(g,pedges[i].b)%2) |
|
602 g.addEdge(pedges[i].a,pedges[i].b); |
|
603 } |
|
604 |
|
605 for(NodeIt n(g);n!=INVALID;++n) |
|
606 if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 ) |
|
607 std::cout << "GEBASZ!!!" << std::endl; |
|
608 |
|
609 for(EdgeIt e(g);e!=INVALID;++e) |
|
610 if(g.u(e)==g.v(e)) |
|
611 std::cout << "LOOP GEBASZ!!!" << std::endl; |
|
612 |
|
613 std::cout << "Number of arcs : " << countEdges(g) << std::endl; |
|
614 |
|
615 std::cout << "Total arc length (euler) : " << totalLen() << std::endl; |
|
616 |
|
617 ListGraph::EdgeMap<Arc> enext(g); |
|
618 { |
|
619 EulerIt<ListGraph> e(g); |
|
620 Arc eo=e; |
|
621 Arc ef=e; |
|
622 // std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl; |
|
623 for(++e;e!=INVALID;++e) |
|
624 { |
|
625 // std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl; |
|
626 enext[eo]=e; |
|
627 eo=e; |
|
628 } |
|
629 enext[eo]=ef; |
|
630 } |
|
631 |
|
632 std::cout << "Creating a tour from that..." << std::endl; |
|
633 |
|
634 int nnum = countNodes(g); |
|
635 int ednum = countEdges(g); |
|
636 |
|
637 for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p]) |
|
638 { |
|
639 // std::cout << "Checking arc " << g.id(p) << std::endl; |
|
640 Arc e=enext[p]; |
|
641 Arc f=enext[e]; |
|
642 Node n2=g.source(f); |
|
643 Node n1=g.oppositeNode(n2,e); |
|
644 Node n3=g.oppositeNode(n2,f); |
|
645 if(countIncEdges(g,n2)>2) |
|
646 { |
|
647 // std::cout << "Remove an Arc" << std::endl; |
|
648 Arc ff=enext[f]; |
|
649 g.erase(e); |
|
650 g.erase(f); |
|
651 if(n1!=n3) |
|
652 { |
|
653 Arc ne=g.direct(g.addEdge(n1,n3),n1); |
|
654 enext[p]=ne; |
|
655 enext[ne]=ff; |
|
656 ednum--; |
|
657 } |
|
658 else { |
|
659 enext[p]=ff; |
|
660 ednum-=2; |
|
661 } |
|
662 } |
|
663 } |
|
664 |
|
665 std::cout << "Total arc length (tour) : " << totalLen() << std::endl; |
|
666 |
|
667 std::cout << "2-opt the tour..." << std::endl; |
|
668 |
|
669 tsp_improve(); |
|
670 |
|
671 std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl; |
|
672 } |
|
673 |
|
674 |
|
675 int main(int argc,const char **argv) |
|
676 { |
|
677 ArgParser ap(argc,argv); |
|
678 |
|
679 // bool eps; |
|
680 bool disc_d, square_d, gauss_d; |
|
681 // bool tsp_a,two_a,tree_a; |
|
682 int num_of_cities=1; |
|
683 double area=1; |
|
684 N=100; |
|
685 // girth=10; |
|
686 std::string ndist("disc"); |
|
687 ap.refOption("n", "Number of nodes (default is 100)", N) |
|
688 .intOption("g", "Girth parameter (default is 10)", 10) |
|
689 .refOption("cities", "Number of cities (default is 1)", num_of_cities) |
|
690 .refOption("area", "Full relative area of the cities (default is 1)", area) |
|
691 .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d) |
|
692 .optionGroup("dist", "disc") |
|
693 .refOption("square", "Nodes are evenly distributed on a unit square", square_d) |
|
694 .optionGroup("dist", "square") |
|
695 .refOption("gauss", |
|
696 "Nodes are located according to a two-dim gauss distribution", |
|
697 gauss_d) |
|
698 .optionGroup("dist", "gauss") |
|
699 // .mandatoryGroup("dist") |
|
700 .onlyOneGroup("dist") |
|
701 .boolOption("eps", "Also generate .eps output (prefix.eps)") |
|
702 .boolOption("dir", "Directed digraph is generated (each arcs are replaced by two directed ones)") |
|
703 .boolOption("2con", "Create a two connected planar digraph") |
|
704 .optionGroup("alg","2con") |
|
705 .boolOption("tree", "Create a min. cost spanning tree") |
|
706 .optionGroup("alg","tree") |
|
707 .boolOption("tsp", "Create a TSP tour") |
|
708 .optionGroup("alg","tsp") |
|
709 .boolOption("tsp2", "Create a TSP tour (tree based)") |
|
710 .optionGroup("alg","tsp2") |
|
711 .boolOption("dela", "Delaunay triangulation digraph") |
|
712 .optionGroup("alg","dela") |
|
713 .onlyOneGroup("alg") |
|
714 .boolOption("rand", "Use time seed for random number generator") |
|
715 .optionGroup("rand", "rand") |
|
716 .intOption("seed", "Random seed", -1) |
|
717 .optionGroup("rand", "seed") |
|
718 .onlyOneGroup("rand") |
|
719 .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'") |
|
720 .run(); |
|
721 |
|
722 if (ap["rand"]) { |
|
723 int seed = time(0); |
|
724 std::cout << "Random number seed: " << seed << std::endl; |
|
725 rnd = Random(seed); |
|
726 } |
|
727 if (ap.given("seed")) { |
|
728 int seed = ap["seed"]; |
|
729 std::cout << "Random number seed: " << seed << std::endl; |
|
730 rnd = Random(seed); |
|
731 } |
|
732 |
|
733 std::string prefix; |
|
734 switch(ap.files().size()) |
|
735 { |
|
736 case 0: |
|
737 prefix="lgf-gen-out"; |
|
738 break; |
|
739 case 1: |
|
740 prefix=ap.files()[0]; |
|
741 break; |
|
742 default: |
|
743 std::cerr << "\nAt most one prefix can be given\n\n"; |
|
744 exit(1); |
|
745 } |
|
746 |
|
747 double sum_sizes=0; |
|
748 std::vector<double> sizes; |
|
749 std::vector<double> cum_sizes; |
|
750 for(int s=0;s<num_of_cities;s++) |
|
751 { |
|
752 // sum_sizes+=rnd.exponential(); |
|
753 double d=rnd(); |
|
754 sum_sizes+=d; |
|
755 sizes.push_back(d); |
|
756 cum_sizes.push_back(sum_sizes); |
|
757 } |
|
758 int i=0; |
|
759 for(int s=0;s<num_of_cities;s++) |
|
760 { |
|
761 Point center=(num_of_cities==1?Point(0,0):rnd.disc()); |
|
762 if(gauss_d) |
|
763 for(;i<N*(cum_sizes[s]/sum_sizes);i++) { |
|
764 Node n=g.addNode(); |
|
765 nodes.push_back(n); |
|
766 coords[n]=center+rnd.gauss2()*area* |
|
767 std::sqrt(sizes[s]/sum_sizes); |
|
768 } |
|
769 else if(square_d) |
|
770 for(;i<N*(cum_sizes[s]/sum_sizes);i++) { |
|
771 Node n=g.addNode(); |
|
772 nodes.push_back(n); |
|
773 coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area* |
|
774 std::sqrt(sizes[s]/sum_sizes); |
|
775 } |
|
776 else if(disc_d || true) |
|
777 for(;i<N*(cum_sizes[s]/sum_sizes);i++) { |
|
778 Node n=g.addNode(); |
|
779 nodes.push_back(n); |
|
780 coords[n]=center+rnd.disc()*area* |
|
781 std::sqrt(sizes[s]/sum_sizes); |
|
782 } |
|
783 } |
|
784 |
|
785 // for (ListGraph::NodeIt n(g); n != INVALID; ++n) { |
|
786 // std::cerr << coords[n] << std::endl; |
|
787 // } |
|
788 |
|
789 if(ap["tsp"]) { |
|
790 tsp(); |
|
791 std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl; |
|
792 } |
|
793 if(ap["tsp2"]) { |
|
794 tsp2(); |
|
795 std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl; |
|
796 } |
|
797 else if(ap["2con"]) { |
|
798 std::cout << "Make triangles\n"; |
|
799 // triangle(); |
|
800 sparseTriangle(ap["g"]); |
|
801 std::cout << "Make it sparser\n"; |
|
802 sparse2(ap["g"]); |
|
803 } |
|
804 else if(ap["tree"]) { |
|
805 minTree(); |
|
806 } |
|
807 else if(ap["dela"]) { |
|
808 delaunay(); |
|
809 } |
|
810 |
|
811 |
|
812 std::cout << "Number of nodes : " << countNodes(g) << std::endl; |
|
813 std::cout << "Number of arcs : " << countEdges(g) << std::endl; |
|
814 double tlen=0; |
|
815 for(EdgeIt e(g);e!=INVALID;++e) |
|
816 tlen+=sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare()); |
|
817 std::cout << "Total arc length : " << tlen << std::endl; |
|
818 |
|
819 if(ap["eps"]) |
|
820 graphToEps(g,prefix+".eps").scaleToA4(). |
|
821 scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false). |
|
822 coords(coords).run(); |
|
823 |
|
824 if(ap["dir"]) |
|
825 DigraphWriter<ListGraph>(g,prefix+".lgf"). |
|
826 nodeMap("coordinates_x",scaleMap(xMap(coords),600)). |
|
827 nodeMap("coordinates_y",scaleMap(yMap(coords),600)). |
|
828 run(); |
|
829 else GraphWriter<ListGraph>(g,prefix+".lgf"). |
|
830 nodeMap("coordinates_x",scaleMap(xMap(coords),600)). |
|
831 nodeMap("coordinates_y",scaleMap(yMap(coords),600)). |
|
832 run(); |
|
833 } |
|
834 |