[523] | 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 |
---|
[584] | 26 | /// \code |
---|
| 27 | /// lgf-gen --help |
---|
| 28 | /// \endcode |
---|
[523] | 29 | /// for more info on the usage. |
---|
| 30 | |
---|
| 31 | #include <algorithm> |
---|
| 32 | #include <set> |
---|
[570] | 33 | #include <ctime> |
---|
[523] | 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 | |
---|
| 65 | double totalLen(){ |
---|
| 66 | double tlen=0; |
---|
| 67 | for(EdgeIt e(g);e!=INVALID;++e) |
---|
[612] | 68 | tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare()); |
---|
[523] | 69 | return tlen; |
---|
| 70 | } |
---|
| 71 | |
---|
| 72 | int tsp_impr_num=0; |
---|
| 73 | |
---|
| 74 | const double EPSILON=1e-8; |
---|
| 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 | |
---|
| 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 | |
---|
| 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 | |
---|
| 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; |
---|
| 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)]) {} |
---|
| 132 | }; |
---|
| 133 | |
---|
| 134 | inline std::ostream& operator<<(std::ostream &os, const Line &l) |
---|
| 135 | { |
---|
| 136 | os << l.a << "->" << l.b; |
---|
| 137 | return os; |
---|
| 138 | } |
---|
| 139 | |
---|
| 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 | |
---|
| 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 | |
---|
| 167 | Part(int p, int c, int n) : prev(p), curr(c), next(n) {} |
---|
| 168 | }; |
---|
| 169 | |
---|
| 170 | inline std::ostream& operator<<(std::ostream& os, const Part& part) { |
---|
| 171 | os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')'; |
---|
| 172 | return os; |
---|
| 173 | } |
---|
| 174 | |
---|
| 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 | |
---|
[612] | 191 | return d / (2 * a) + std::sqrt((d * d + e * e) / (4 * a * a) + f / a); |
---|
[523] | 192 | } |
---|
| 193 | |
---|
| 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 | |
---|
| 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(); |
---|
[612] | 209 | return (b - std::sqrt(d)) / a; |
---|
[523] | 210 | } |
---|
| 211 | |
---|
| 212 | struct YLess { |
---|
| 213 | |
---|
| 214 | |
---|
| 215 | YLess(const std::vector<Point>& points, double& sweep) |
---|
| 216 | : _points(points), _sweep(sweep) {} |
---|
| 217 | |
---|
| 218 | 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 | |
---|
| 256 | BeachIt(Beach::iterator iter) : it(iter) {} |
---|
| 257 | }; |
---|
| 258 | |
---|
| 259 | } |
---|
| 260 | |
---|
| 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 | spikeheap.erase(bit->second); |
---|
| 333 | } |
---|
| 334 | |
---|
| 335 | int prev = bit->first.prev; |
---|
| 336 | int curr = bit->first.curr; |
---|
| 337 | int next = bit->first.next; |
---|
| 338 | |
---|
| 339 | beach.erase(bit); |
---|
| 340 | |
---|
| 341 | SpikeHeap::iterator pit = spikeheap.end(); |
---|
| 342 | if (prev != -1 && |
---|
| 343 | circle_form(points[prev], points[curr], points[site])) { |
---|
| 344 | double x = circle_point(points[prev], points[curr], points[site]); |
---|
| 345 | pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); |
---|
| 346 | pit->second.it = |
---|
| 347 | beach.insert(std::make_pair(Part(prev, curr, site), pit)); |
---|
| 348 | } else { |
---|
| 349 | beach.insert(std::make_pair(Part(prev, curr, site), pit)); |
---|
| 350 | } |
---|
| 351 | |
---|
| 352 | beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end())); |
---|
| 353 | |
---|
| 354 | SpikeHeap::iterator nit = spikeheap.end(); |
---|
| 355 | if (next != -1 && |
---|
| 356 | circle_form(points[site], points[curr],points[next])) { |
---|
| 357 | double x = circle_point(points[site], points[curr], points[next]); |
---|
| 358 | nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); |
---|
| 359 | nit->second.it = |
---|
| 360 | beach.insert(std::make_pair(Part(site, curr, next), nit)); |
---|
| 361 | } else { |
---|
| 362 | beach.insert(std::make_pair(Part(site, curr, next), nit)); |
---|
| 363 | } |
---|
| 364 | |
---|
| 365 | ++siteindex; |
---|
| 366 | } else { |
---|
| 367 | sweep = spit->first; |
---|
| 368 | |
---|
| 369 | Beach::iterator bit = spit->second.it; |
---|
| 370 | |
---|
| 371 | int prev = bit->first.prev; |
---|
| 372 | int curr = bit->first.curr; |
---|
| 373 | int next = bit->first.next; |
---|
| 374 | |
---|
| 375 | { |
---|
| 376 | std::pair<int, int> arc; |
---|
| 377 | |
---|
| 378 | arc = prev < curr ? |
---|
| 379 | std::make_pair(prev, curr) : std::make_pair(curr, prev); |
---|
| 380 | |
---|
| 381 | if (arcs.find(arc) == arcs.end()) { |
---|
| 382 | arcs.insert(arc); |
---|
| 383 | g.addEdge(nodes[prev], nodes[curr]); |
---|
| 384 | ++cnt; |
---|
| 385 | } |
---|
| 386 | |
---|
| 387 | arc = curr < next ? |
---|
| 388 | std::make_pair(curr, next) : std::make_pair(next, curr); |
---|
| 389 | |
---|
| 390 | if (arcs.find(arc) == arcs.end()) { |
---|
| 391 | arcs.insert(arc); |
---|
| 392 | g.addEdge(nodes[curr], nodes[next]); |
---|
| 393 | ++cnt; |
---|
| 394 | } |
---|
| 395 | } |
---|
| 396 | |
---|
| 397 | Beach::iterator pbit = bit; --pbit; |
---|
| 398 | int ppv = pbit->first.prev; |
---|
| 399 | Beach::iterator nbit = bit; ++nbit; |
---|
| 400 | int nnt = nbit->first.next; |
---|
| 401 | |
---|
| 402 | if (bit->second != spikeheap.end()) spikeheap.erase(bit->second); |
---|
| 403 | if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second); |
---|
| 404 | if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second); |
---|
| 405 | |
---|
| 406 | beach.erase(nbit); |
---|
| 407 | beach.erase(bit); |
---|
| 408 | beach.erase(pbit); |
---|
| 409 | |
---|
| 410 | SpikeHeap::iterator pit = spikeheap.end(); |
---|
| 411 | if (ppv != -1 && ppv != next && |
---|
| 412 | circle_form(points[ppv], points[prev], points[next])) { |
---|
| 413 | double x = circle_point(points[ppv], points[prev], points[next]); |
---|
| 414 | if (x < sweep) x = sweep; |
---|
| 415 | pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); |
---|
| 416 | pit->second.it = |
---|
| 417 | beach.insert(std::make_pair(Part(ppv, prev, next), pit)); |
---|
| 418 | } else { |
---|
| 419 | beach.insert(std::make_pair(Part(ppv, prev, next), pit)); |
---|
| 420 | } |
---|
| 421 | |
---|
| 422 | SpikeHeap::iterator nit = spikeheap.end(); |
---|
| 423 | if (nnt != -1 && prev != nnt && |
---|
| 424 | circle_form(points[prev], points[next], points[nnt])) { |
---|
| 425 | double x = circle_point(points[prev], points[next], points[nnt]); |
---|
| 426 | if (x < sweep) x = sweep; |
---|
| 427 | nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); |
---|
| 428 | nit->second.it = |
---|
| 429 | beach.insert(std::make_pair(Part(prev, next, nnt), nit)); |
---|
| 430 | } else { |
---|
| 431 | beach.insert(std::make_pair(Part(prev, next, nnt), nit)); |
---|
| 432 | } |
---|
| 433 | |
---|
| 434 | } |
---|
| 435 | } |
---|
| 436 | |
---|
| 437 | for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) { |
---|
| 438 | int curr = it->first.curr; |
---|
| 439 | int next = it->first.next; |
---|
| 440 | |
---|
| 441 | if (next == -1) continue; |
---|
| 442 | |
---|
| 443 | std::pair<int, int> arc; |
---|
| 444 | |
---|
| 445 | arc = curr < next ? |
---|
| 446 | std::make_pair(curr, next) : std::make_pair(next, curr); |
---|
| 447 | |
---|
| 448 | if (arcs.find(arc) == arcs.end()) { |
---|
| 449 | arcs.insert(arc); |
---|
| 450 | g.addEdge(nodes[curr], nodes[next]); |
---|
| 451 | ++cnt; |
---|
| 452 | } |
---|
| 453 | } |
---|
| 454 | } |
---|
| 455 | |
---|
| 456 | void sparse(int d) |
---|
| 457 | { |
---|
| 458 | Counter cnt("Number of arcs removed: "); |
---|
| 459 | Bfs<ListGraph> bfs(g); |
---|
| 460 | for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin(); |
---|
| 461 | ei!=arcs.rend();++ei) |
---|
| 462 | { |
---|
| 463 | Node a=g.u(*ei); |
---|
| 464 | Node b=g.v(*ei); |
---|
| 465 | g.erase(*ei); |
---|
| 466 | bfs.run(a,b); |
---|
| 467 | if(bfs.predArc(b)==INVALID || bfs.dist(b)>d) |
---|
| 468 | g.addEdge(a,b); |
---|
| 469 | else cnt++; |
---|
| 470 | } |
---|
| 471 | } |
---|
| 472 | |
---|
| 473 | void sparse2(int d) |
---|
| 474 | { |
---|
| 475 | Counter cnt("Number of arcs removed: "); |
---|
| 476 | for(std::vector<Edge>::reverse_iterator ei=arcs.rbegin(); |
---|
| 477 | ei!=arcs.rend();++ei) |
---|
| 478 | { |
---|
| 479 | Node a=g.u(*ei); |
---|
| 480 | Node b=g.v(*ei); |
---|
| 481 | g.erase(*ei); |
---|
| 482 | ConstMap<Arc,int> cegy(1); |
---|
[623] | 483 | Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy); |
---|
| 484 | int k=sur.run(a,b,2); |
---|
[523] | 485 | if(k<2 || sur.totalLength()>d) |
---|
| 486 | g.addEdge(a,b); |
---|
| 487 | else cnt++; |
---|
| 488 | // else std::cout << "Remove arc " << g.id(a) << "-" << g.id(b) << '\n'; |
---|
| 489 | } |
---|
| 490 | } |
---|
| 491 | |
---|
| 492 | void sparseTriangle(int d) |
---|
| 493 | { |
---|
| 494 | Counter cnt("Number of arcs added: "); |
---|
| 495 | std::vector<Parc> pedges; |
---|
| 496 | for(NodeIt n(g);n!=INVALID;++n) |
---|
| 497 | for(NodeIt m=++(NodeIt(n));m!=INVALID;++m) |
---|
| 498 | { |
---|
| 499 | Parc p; |
---|
| 500 | p.a=n; |
---|
| 501 | p.b=m; |
---|
| 502 | p.len=(coords[m]-coords[n]).normSquare(); |
---|
| 503 | pedges.push_back(p); |
---|
| 504 | } |
---|
| 505 | std::sort(pedges.begin(),pedges.end(),pedgeLess); |
---|
| 506 | for(std::vector<Parc>::iterator pi=pedges.begin();pi!=pedges.end();++pi) |
---|
| 507 | { |
---|
| 508 | Line li(pi->a,pi->b); |
---|
| 509 | EdgeIt e(g); |
---|
| 510 | for(;e!=INVALID && !cross(e,li);++e) ; |
---|
| 511 | Edge ne; |
---|
| 512 | if(e==INVALID) { |
---|
| 513 | ConstMap<Arc,int> cegy(1); |
---|
[623] | 514 | Suurballe<ListGraph,ConstMap<Arc,int> > sur(g,cegy); |
---|
| 515 | int k=sur.run(pi->a,pi->b,2); |
---|
[523] | 516 | if(k<2 || sur.totalLength()>d) |
---|
| 517 | { |
---|
| 518 | ne=g.addEdge(pi->a,pi->b); |
---|
| 519 | arcs.push_back(ne); |
---|
| 520 | cnt++; |
---|
| 521 | } |
---|
| 522 | } |
---|
| 523 | } |
---|
| 524 | } |
---|
| 525 | |
---|
| 526 | template <typename Graph, typename CoordMap> |
---|
| 527 | class LengthSquareMap { |
---|
| 528 | public: |
---|
| 529 | typedef typename Graph::Edge Key; |
---|
| 530 | typedef typename CoordMap::Value::Value Value; |
---|
| 531 | |
---|
| 532 | LengthSquareMap(const Graph& graph, const CoordMap& coords) |
---|
| 533 | : _graph(graph), _coords(coords) {} |
---|
| 534 | |
---|
| 535 | Value operator[](const Key& key) const { |
---|
| 536 | return (_coords[_graph.v(key)] - |
---|
| 537 | _coords[_graph.u(key)]).normSquare(); |
---|
| 538 | } |
---|
| 539 | |
---|
| 540 | private: |
---|
| 541 | |
---|
| 542 | const Graph& _graph; |
---|
| 543 | const CoordMap& _coords; |
---|
| 544 | }; |
---|
| 545 | |
---|
| 546 | void minTree() { |
---|
| 547 | std::vector<Parc> pedges; |
---|
| 548 | Timer T; |
---|
| 549 | std::cout << T.realTime() << "s: Creating delaunay triangulation...\n"; |
---|
| 550 | delaunay(); |
---|
| 551 | std::cout << T.realTime() << "s: Calculating spanning tree...\n"; |
---|
| 552 | LengthSquareMap<ListGraph, ListGraph::NodeMap<Point> > ls(g, coords); |
---|
| 553 | ListGraph::EdgeMap<bool> tree(g); |
---|
| 554 | kruskal(g, ls, tree); |
---|
| 555 | std::cout << T.realTime() << "s: Removing non tree arcs...\n"; |
---|
| 556 | std::vector<Edge> remove; |
---|
| 557 | for (EdgeIt e(g); e != INVALID; ++e) { |
---|
| 558 | if (!tree[e]) remove.push_back(e); |
---|
| 559 | } |
---|
| 560 | for(int i = 0; i < int(remove.size()); ++i) { |
---|
| 561 | g.erase(remove[i]); |
---|
| 562 | } |
---|
| 563 | std::cout << T.realTime() << "s: Done\n"; |
---|
| 564 | } |
---|
| 565 | |
---|
| 566 | void tsp2() |
---|
| 567 | { |
---|
| 568 | std::cout << "Find a tree..." << std::endl; |
---|
| 569 | |
---|
| 570 | minTree(); |
---|
| 571 | |
---|
| 572 | std::cout << "Total arc length (tree) : " << totalLen() << std::endl; |
---|
| 573 | |
---|
| 574 | std::cout << "Make it Euler..." << std::endl; |
---|
| 575 | |
---|
| 576 | { |
---|
| 577 | std::vector<Node> leafs; |
---|
| 578 | for(NodeIt n(g);n!=INVALID;++n) |
---|
| 579 | if(countIncEdges(g,n)%2==1) leafs.push_back(n); |
---|
| 580 | |
---|
| 581 | // for(unsigned int i=0;i<leafs.size();i+=2) |
---|
| 582 | // g.addArc(leafs[i],leafs[i+1]); |
---|
| 583 | |
---|
| 584 | std::vector<Parc> pedges; |
---|
| 585 | for(unsigned int i=0;i<leafs.size()-1;i++) |
---|
| 586 | for(unsigned int j=i+1;j<leafs.size();j++) |
---|
| 587 | { |
---|
| 588 | Node n=leafs[i]; |
---|
| 589 | Node m=leafs[j]; |
---|
| 590 | Parc p; |
---|
| 591 | p.a=n; |
---|
| 592 | p.b=m; |
---|
| 593 | p.len=(coords[m]-coords[n]).normSquare(); |
---|
| 594 | pedges.push_back(p); |
---|
| 595 | } |
---|
| 596 | std::sort(pedges.begin(),pedges.end(),pedgeLess); |
---|
| 597 | for(unsigned int i=0;i<pedges.size();i++) |
---|
| 598 | if(countIncEdges(g,pedges[i].a)%2 && |
---|
| 599 | countIncEdges(g,pedges[i].b)%2) |
---|
| 600 | g.addEdge(pedges[i].a,pedges[i].b); |
---|
| 601 | } |
---|
| 602 | |
---|
| 603 | for(NodeIt n(g);n!=INVALID;++n) |
---|
| 604 | if(countIncEdges(g,n)%2 || countIncEdges(g,n)==0 ) |
---|
| 605 | std::cout << "GEBASZ!!!" << std::endl; |
---|
| 606 | |
---|
| 607 | for(EdgeIt e(g);e!=INVALID;++e) |
---|
| 608 | if(g.u(e)==g.v(e)) |
---|
| 609 | std::cout << "LOOP GEBASZ!!!" << std::endl; |
---|
| 610 | |
---|
| 611 | std::cout << "Number of arcs : " << countEdges(g) << std::endl; |
---|
| 612 | |
---|
| 613 | std::cout << "Total arc length (euler) : " << totalLen() << std::endl; |
---|
| 614 | |
---|
| 615 | ListGraph::EdgeMap<Arc> enext(g); |
---|
| 616 | { |
---|
| 617 | EulerIt<ListGraph> e(g); |
---|
| 618 | Arc eo=e; |
---|
| 619 | Arc ef=e; |
---|
| 620 | // std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl; |
---|
| 621 | for(++e;e!=INVALID;++e) |
---|
| 622 | { |
---|
| 623 | // std::cout << "Tour arc: " << g.id(Edge(e)) << std::endl; |
---|
| 624 | enext[eo]=e; |
---|
| 625 | eo=e; |
---|
| 626 | } |
---|
| 627 | enext[eo]=ef; |
---|
| 628 | } |
---|
| 629 | |
---|
| 630 | std::cout << "Creating a tour from that..." << std::endl; |
---|
| 631 | |
---|
| 632 | int nnum = countNodes(g); |
---|
| 633 | int ednum = countEdges(g); |
---|
| 634 | |
---|
| 635 | for(Arc p=enext[EdgeIt(g)];ednum>nnum;p=enext[p]) |
---|
| 636 | { |
---|
| 637 | // std::cout << "Checking arc " << g.id(p) << std::endl; |
---|
| 638 | Arc e=enext[p]; |
---|
| 639 | Arc f=enext[e]; |
---|
| 640 | Node n2=g.source(f); |
---|
| 641 | Node n1=g.oppositeNode(n2,e); |
---|
| 642 | Node n3=g.oppositeNode(n2,f); |
---|
| 643 | if(countIncEdges(g,n2)>2) |
---|
| 644 | { |
---|
| 645 | // std::cout << "Remove an Arc" << std::endl; |
---|
| 646 | Arc ff=enext[f]; |
---|
| 647 | g.erase(e); |
---|
| 648 | g.erase(f); |
---|
| 649 | if(n1!=n3) |
---|
| 650 | { |
---|
| 651 | Arc ne=g.direct(g.addEdge(n1,n3),n1); |
---|
| 652 | enext[p]=ne; |
---|
| 653 | enext[ne]=ff; |
---|
| 654 | ednum--; |
---|
| 655 | } |
---|
| 656 | else { |
---|
| 657 | enext[p]=ff; |
---|
| 658 | ednum-=2; |
---|
| 659 | } |
---|
| 660 | } |
---|
| 661 | } |
---|
| 662 | |
---|
| 663 | std::cout << "Total arc length (tour) : " << totalLen() << std::endl; |
---|
| 664 | |
---|
| 665 | std::cout << "2-opt the tour..." << std::endl; |
---|
| 666 | |
---|
| 667 | tsp_improve(); |
---|
| 668 | |
---|
| 669 | std::cout << "Total arc length (2-opt tour) : " << totalLen() << std::endl; |
---|
| 670 | } |
---|
| 671 | |
---|
| 672 | |
---|
| 673 | int main(int argc,const char **argv) |
---|
| 674 | { |
---|
| 675 | ArgParser ap(argc,argv); |
---|
| 676 | |
---|
| 677 | // bool eps; |
---|
| 678 | bool disc_d, square_d, gauss_d; |
---|
| 679 | // bool tsp_a,two_a,tree_a; |
---|
| 680 | int num_of_cities=1; |
---|
| 681 | double area=1; |
---|
| 682 | N=100; |
---|
| 683 | // girth=10; |
---|
| 684 | std::string ndist("disc"); |
---|
| 685 | ap.refOption("n", "Number of nodes (default is 100)", N) |
---|
| 686 | .intOption("g", "Girth parameter (default is 10)", 10) |
---|
| 687 | .refOption("cities", "Number of cities (default is 1)", num_of_cities) |
---|
| 688 | .refOption("area", "Full relative area of the cities (default is 1)", area) |
---|
| 689 | .refOption("disc", "Nodes are evenly distributed on a unit disc (default)",disc_d) |
---|
| 690 | .optionGroup("dist", "disc") |
---|
| 691 | .refOption("square", "Nodes are evenly distributed on a unit square", square_d) |
---|
| 692 | .optionGroup("dist", "square") |
---|
| 693 | .refOption("gauss", |
---|
| 694 | "Nodes are located according to a two-dim gauss distribution", |
---|
| 695 | gauss_d) |
---|
| 696 | .optionGroup("dist", "gauss") |
---|
| 697 | // .mandatoryGroup("dist") |
---|
| 698 | .onlyOneGroup("dist") |
---|
| 699 | .boolOption("eps", "Also generate .eps output (prefix.eps)") |
---|
[524] | 700 | .boolOption("nonodes", "Draw the edges only in the generated .eps") |
---|
[523] | 701 | .boolOption("dir", "Directed digraph is generated (each arcs are replaced by two directed ones)") |
---|
| 702 | .boolOption("2con", "Create a two connected planar digraph") |
---|
| 703 | .optionGroup("alg","2con") |
---|
| 704 | .boolOption("tree", "Create a min. cost spanning tree") |
---|
| 705 | .optionGroup("alg","tree") |
---|
| 706 | .boolOption("tsp", "Create a TSP tour") |
---|
| 707 | .optionGroup("alg","tsp") |
---|
| 708 | .boolOption("tsp2", "Create a TSP tour (tree based)") |
---|
| 709 | .optionGroup("alg","tsp2") |
---|
| 710 | .boolOption("dela", "Delaunay triangulation digraph") |
---|
| 711 | .optionGroup("alg","dela") |
---|
| 712 | .onlyOneGroup("alg") |
---|
| 713 | .boolOption("rand", "Use time seed for random number generator") |
---|
| 714 | .optionGroup("rand", "rand") |
---|
| 715 | .intOption("seed", "Random seed", -1) |
---|
| 716 | .optionGroup("rand", "seed") |
---|
| 717 | .onlyOneGroup("rand") |
---|
| 718 | .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'") |
---|
| 719 | .run(); |
---|
| 720 | |
---|
| 721 | if (ap["rand"]) { |
---|
[616] | 722 | int seed = int(time(0)); |
---|
[523] | 723 | std::cout << "Random number seed: " << seed << std::endl; |
---|
| 724 | rnd = Random(seed); |
---|
| 725 | } |
---|
| 726 | if (ap.given("seed")) { |
---|
| 727 | int seed = ap["seed"]; |
---|
| 728 | std::cout << "Random number seed: " << seed << std::endl; |
---|
| 729 | rnd = Random(seed); |
---|
| 730 | } |
---|
| 731 | |
---|
| 732 | std::string prefix; |
---|
| 733 | switch(ap.files().size()) |
---|
| 734 | { |
---|
| 735 | case 0: |
---|
| 736 | prefix="lgf-gen-out"; |
---|
| 737 | break; |
---|
| 738 | case 1: |
---|
| 739 | prefix=ap.files()[0]; |
---|
| 740 | break; |
---|
| 741 | default: |
---|
| 742 | std::cerr << "\nAt most one prefix can be given\n\n"; |
---|
| 743 | exit(1); |
---|
| 744 | } |
---|
| 745 | |
---|
| 746 | double sum_sizes=0; |
---|
| 747 | std::vector<double> sizes; |
---|
| 748 | std::vector<double> cum_sizes; |
---|
| 749 | for(int s=0;s<num_of_cities;s++) |
---|
| 750 | { |
---|
| 751 | // sum_sizes+=rnd.exponential(); |
---|
| 752 | double d=rnd(); |
---|
| 753 | sum_sizes+=d; |
---|
| 754 | sizes.push_back(d); |
---|
| 755 | cum_sizes.push_back(sum_sizes); |
---|
| 756 | } |
---|
| 757 | int i=0; |
---|
| 758 | for(int s=0;s<num_of_cities;s++) |
---|
| 759 | { |
---|
| 760 | Point center=(num_of_cities==1?Point(0,0):rnd.disc()); |
---|
| 761 | if(gauss_d) |
---|
| 762 | for(;i<N*(cum_sizes[s]/sum_sizes);i++) { |
---|
| 763 | Node n=g.addNode(); |
---|
| 764 | nodes.push_back(n); |
---|
| 765 | coords[n]=center+rnd.gauss2()*area* |
---|
| 766 | std::sqrt(sizes[s]/sum_sizes); |
---|
| 767 | } |
---|
| 768 | else if(square_d) |
---|
| 769 | for(;i<N*(cum_sizes[s]/sum_sizes);i++) { |
---|
| 770 | Node n=g.addNode(); |
---|
| 771 | nodes.push_back(n); |
---|
| 772 | coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area* |
---|
| 773 | std::sqrt(sizes[s]/sum_sizes); |
---|
| 774 | } |
---|
| 775 | else if(disc_d || true) |
---|
| 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.disc()*area* |
---|
| 780 | std::sqrt(sizes[s]/sum_sizes); |
---|
| 781 | } |
---|
| 782 | } |
---|
| 783 | |
---|
| 784 | // for (ListGraph::NodeIt n(g); n != INVALID; ++n) { |
---|
| 785 | // std::cerr << coords[n] << std::endl; |
---|
| 786 | // } |
---|
| 787 | |
---|
| 788 | if(ap["tsp"]) { |
---|
| 789 | tsp(); |
---|
| 790 | std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl; |
---|
| 791 | } |
---|
| 792 | if(ap["tsp2"]) { |
---|
| 793 | tsp2(); |
---|
| 794 | std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl; |
---|
| 795 | } |
---|
| 796 | else if(ap["2con"]) { |
---|
| 797 | std::cout << "Make triangles\n"; |
---|
| 798 | // triangle(); |
---|
| 799 | sparseTriangle(ap["g"]); |
---|
| 800 | std::cout << "Make it sparser\n"; |
---|
| 801 | sparse2(ap["g"]); |
---|
| 802 | } |
---|
| 803 | else if(ap["tree"]) { |
---|
| 804 | minTree(); |
---|
| 805 | } |
---|
| 806 | else if(ap["dela"]) { |
---|
| 807 | delaunay(); |
---|
| 808 | } |
---|
| 809 | |
---|
| 810 | |
---|
| 811 | std::cout << "Number of nodes : " << countNodes(g) << std::endl; |
---|
| 812 | std::cout << "Number of arcs : " << countEdges(g) << std::endl; |
---|
| 813 | double tlen=0; |
---|
| 814 | for(EdgeIt e(g);e!=INVALID;++e) |
---|
[612] | 815 | tlen+=std::sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare()); |
---|
[523] | 816 | std::cout << "Total arc length : " << tlen << std::endl; |
---|
| 817 | |
---|
| 818 | if(ap["eps"]) |
---|
| 819 | graphToEps(g,prefix+".eps").scaleToA4(). |
---|
| 820 | scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false). |
---|
[524] | 821 | coords(coords).hideNodes(ap.given("nonodes")).run(); |
---|
[523] | 822 | |
---|
| 823 | if(ap["dir"]) |
---|
| 824 | DigraphWriter<ListGraph>(g,prefix+".lgf"). |
---|
| 825 | nodeMap("coordinates_x",scaleMap(xMap(coords),600)). |
---|
| 826 | nodeMap("coordinates_y",scaleMap(yMap(coords),600)). |
---|
| 827 | run(); |
---|
| 828 | else GraphWriter<ListGraph>(g,prefix+".lgf"). |
---|
| 829 | nodeMap("coordinates_x",scaleMap(xMap(coords),600)). |
---|
| 830 | nodeMap("coordinates_y",scaleMap(yMap(coords),600)). |
---|
| 831 | run(); |
---|
| 832 | } |
---|
| 833 | |
---|