[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 |
---|
| 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)") |
---|
[524] | 702 | .boolOption("nonodes", "Draw the edges only in the generated .eps") |
---|
[523] | 703 | .boolOption("dir", "Directed digraph is generated (each arcs are replaced by two directed ones)") |
---|
| 704 | .boolOption("2con", "Create a two connected planar digraph") |
---|
| 705 | .optionGroup("alg","2con") |
---|
| 706 | .boolOption("tree", "Create a min. cost spanning tree") |
---|
| 707 | .optionGroup("alg","tree") |
---|
| 708 | .boolOption("tsp", "Create a TSP tour") |
---|
| 709 | .optionGroup("alg","tsp") |
---|
| 710 | .boolOption("tsp2", "Create a TSP tour (tree based)") |
---|
| 711 | .optionGroup("alg","tsp2") |
---|
| 712 | .boolOption("dela", "Delaunay triangulation digraph") |
---|
| 713 | .optionGroup("alg","dela") |
---|
| 714 | .onlyOneGroup("alg") |
---|
| 715 | .boolOption("rand", "Use time seed for random number generator") |
---|
| 716 | .optionGroup("rand", "rand") |
---|
| 717 | .intOption("seed", "Random seed", -1) |
---|
| 718 | .optionGroup("rand", "seed") |
---|
| 719 | .onlyOneGroup("rand") |
---|
| 720 | .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'") |
---|
| 721 | .run(); |
---|
| 722 | |
---|
| 723 | if (ap["rand"]) { |
---|
| 724 | int seed = time(0); |
---|
| 725 | std::cout << "Random number seed: " << seed << std::endl; |
---|
| 726 | rnd = Random(seed); |
---|
| 727 | } |
---|
| 728 | if (ap.given("seed")) { |
---|
| 729 | int seed = ap["seed"]; |
---|
| 730 | std::cout << "Random number seed: " << seed << std::endl; |
---|
| 731 | rnd = Random(seed); |
---|
| 732 | } |
---|
| 733 | |
---|
| 734 | std::string prefix; |
---|
| 735 | switch(ap.files().size()) |
---|
| 736 | { |
---|
| 737 | case 0: |
---|
| 738 | prefix="lgf-gen-out"; |
---|
| 739 | break; |
---|
| 740 | case 1: |
---|
| 741 | prefix=ap.files()[0]; |
---|
| 742 | break; |
---|
| 743 | default: |
---|
| 744 | std::cerr << "\nAt most one prefix can be given\n\n"; |
---|
| 745 | exit(1); |
---|
| 746 | } |
---|
| 747 | |
---|
| 748 | double sum_sizes=0; |
---|
| 749 | std::vector<double> sizes; |
---|
| 750 | std::vector<double> cum_sizes; |
---|
| 751 | for(int s=0;s<num_of_cities;s++) |
---|
| 752 | { |
---|
| 753 | // sum_sizes+=rnd.exponential(); |
---|
| 754 | double d=rnd(); |
---|
| 755 | sum_sizes+=d; |
---|
| 756 | sizes.push_back(d); |
---|
| 757 | cum_sizes.push_back(sum_sizes); |
---|
| 758 | } |
---|
| 759 | int i=0; |
---|
| 760 | for(int s=0;s<num_of_cities;s++) |
---|
| 761 | { |
---|
| 762 | Point center=(num_of_cities==1?Point(0,0):rnd.disc()); |
---|
| 763 | if(gauss_d) |
---|
| 764 | for(;i<N*(cum_sizes[s]/sum_sizes);i++) { |
---|
| 765 | Node n=g.addNode(); |
---|
| 766 | nodes.push_back(n); |
---|
| 767 | coords[n]=center+rnd.gauss2()*area* |
---|
| 768 | std::sqrt(sizes[s]/sum_sizes); |
---|
| 769 | } |
---|
| 770 | else if(square_d) |
---|
| 771 | for(;i<N*(cum_sizes[s]/sum_sizes);i++) { |
---|
| 772 | Node n=g.addNode(); |
---|
| 773 | nodes.push_back(n); |
---|
| 774 | coords[n]=center+Point(rnd()*2-1,rnd()*2-1)*area* |
---|
| 775 | std::sqrt(sizes[s]/sum_sizes); |
---|
| 776 | } |
---|
| 777 | else if(disc_d || true) |
---|
| 778 | for(;i<N*(cum_sizes[s]/sum_sizes);i++) { |
---|
| 779 | Node n=g.addNode(); |
---|
| 780 | nodes.push_back(n); |
---|
| 781 | coords[n]=center+rnd.disc()*area* |
---|
| 782 | std::sqrt(sizes[s]/sum_sizes); |
---|
| 783 | } |
---|
| 784 | } |
---|
| 785 | |
---|
| 786 | // for (ListGraph::NodeIt n(g); n != INVALID; ++n) { |
---|
| 787 | // std::cerr << coords[n] << std::endl; |
---|
| 788 | // } |
---|
| 789 | |
---|
| 790 | if(ap["tsp"]) { |
---|
| 791 | tsp(); |
---|
| 792 | std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl; |
---|
| 793 | } |
---|
| 794 | if(ap["tsp2"]) { |
---|
| 795 | tsp2(); |
---|
| 796 | std::cout << "#2-opt improvements: " << tsp_impr_num << std::endl; |
---|
| 797 | } |
---|
| 798 | else if(ap["2con"]) { |
---|
| 799 | std::cout << "Make triangles\n"; |
---|
| 800 | // triangle(); |
---|
| 801 | sparseTriangle(ap["g"]); |
---|
| 802 | std::cout << "Make it sparser\n"; |
---|
| 803 | sparse2(ap["g"]); |
---|
| 804 | } |
---|
| 805 | else if(ap["tree"]) { |
---|
| 806 | minTree(); |
---|
| 807 | } |
---|
| 808 | else if(ap["dela"]) { |
---|
| 809 | delaunay(); |
---|
| 810 | } |
---|
| 811 | |
---|
| 812 | |
---|
| 813 | std::cout << "Number of nodes : " << countNodes(g) << std::endl; |
---|
| 814 | std::cout << "Number of arcs : " << countEdges(g) << std::endl; |
---|
| 815 | double tlen=0; |
---|
| 816 | for(EdgeIt e(g);e!=INVALID;++e) |
---|
| 817 | tlen+=sqrt((coords[g.v(e)]-coords[g.u(e)]).normSquare()); |
---|
| 818 | std::cout << "Total arc length : " << tlen << std::endl; |
---|
| 819 | |
---|
| 820 | if(ap["eps"]) |
---|
| 821 | graphToEps(g,prefix+".eps").scaleToA4(). |
---|
| 822 | scale(600).nodeScale(.005).arcWidthScale(.001).preScale(false). |
---|
[524] | 823 | coords(coords).hideNodes(ap.given("nonodes")).run(); |
---|
[523] | 824 | |
---|
| 825 | if(ap["dir"]) |
---|
| 826 | DigraphWriter<ListGraph>(g,prefix+".lgf"). |
---|
| 827 | nodeMap("coordinates_x",scaleMap(xMap(coords),600)). |
---|
| 828 | nodeMap("coordinates_y",scaleMap(yMap(coords),600)). |
---|
| 829 | run(); |
---|
| 830 | else GraphWriter<ListGraph>(g,prefix+".lgf"). |
---|
| 831 | nodeMap("coordinates_x",scaleMap(xMap(coords),600)). |
---|
| 832 | nodeMap("coordinates_y",scaleMap(yMap(coords),600)). |
---|
| 833 | run(); |
---|
| 834 | } |
---|
| 835 | |
---|