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