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