Changeset 2447:260ce674cc65 in lemon-0.x
- Timestamp:
- 06/05/07 13:49:19 (17 years ago)
- Branch:
- default
- Phase:
- public
- Convert:
- svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3284
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
tools/lgf-gen.cc
r2446 r2447 30 30 #include <cmath> 31 31 #include <algorithm> 32 #include <lemon/ unionfind.h>32 #include <lemon/kruskal.h> 33 33 #include <lemon/time_measure.h> 34 34 … … 147 147 std::vector<UEdge> edges; 148 148 149 void triangle() 150 { 149 namespace _delaunay_bits { 150 151 struct Part { 152 int prev, curr, next; 153 154 Part(int p, int c, int n) : prev(p), curr(c), next(n) {} 155 }; 156 157 inline std::ostream& operator<<(std::ostream& os, const Part& part) { 158 os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')'; 159 return os; 160 } 161 162 inline double circle_point(const Point& p, const Point& q, const Point& r) { 163 double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y); 164 if (a == 0) return std::numeric_limits<double>::quiet_NaN(); 165 166 double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) + 167 (q.x * q.x + q.y * q.y) * (r.y - p.y) + 168 (r.x * r.x + r.y * r.y) * (p.y - q.y); 169 170 double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) + 171 (q.x * q.x + q.y * q.y) * (r.x - p.x) + 172 (r.x * r.x + r.y * r.y) * (p.x - q.x); 173 174 double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) + 175 (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) + 176 (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y); 177 178 return d / (2 * a) + sqrt((d * d + e * e) / (4 * a * a) + f / a); 179 } 180 181 inline bool circle_form(const Point& p, const Point& q, const Point& r) { 182 return rot90(q - p) * (r - q) < 0.0; 183 } 184 185 inline double intersection(const Point& p, const Point& q, double sx) { 186 const double epsilon = 1e-8; 187 188 if (p.x == q.x) return (p.y + q.y) / 2.0; 189 190 if (sx < p.x + epsilon) return p.y; 191 if (sx < q.x + epsilon) return q.y; 192 193 double a = q.x - p.x; 194 double b = (q.x - sx) * p.y - (p.x - sx) * q.y; 195 double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare(); 196 return (b - sqrt(d)) / a; 197 } 198 199 struct YLess { 200 201 202 YLess(const std::vector<Point>& points, double& sweep) 203 : _points(points), _sweep(sweep) {} 204 205 bool operator()(const Part& l, const Part& r) const { 206 const double epsilon = 1e-8; 207 208 // std::cerr << l << " vs " << r << std::endl; 209 double lbx = l.prev != -1 ? 210 intersection(_points[l.prev], _points[l.curr], _sweep) : 211 - std::numeric_limits<double>::infinity(); 212 double rbx = r.prev != -1 ? 213 intersection(_points[r.prev], _points[r.curr], _sweep) : 214 - std::numeric_limits<double>::infinity(); 215 double lex = l.next != -1 ? 216 intersection(_points[l.curr], _points[l.next], _sweep) : 217 std::numeric_limits<double>::infinity(); 218 double rex = r.next != -1 ? 219 intersection(_points[r.curr], _points[r.next], _sweep) : 220 std::numeric_limits<double>::infinity(); 221 222 if (lbx > lex) std::swap(lbx, lex); 223 if (rbx > rex) std::swap(rbx, rex); 224 225 if (lex < epsilon + rex && lbx + epsilon < rex) return true; 226 if (rex < epsilon + lex && rbx + epsilon < lex) return false; 227 return lex < rex; 228 } 229 230 const std::vector<Point>& _points; 231 double& _sweep; 232 }; 233 234 struct BeachIt; 235 236 typedef std::multimap<double, BeachIt> SpikeHeap; 237 238 typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach; 239 240 struct BeachIt { 241 Beach::iterator it; 242 243 BeachIt(Beach::iterator iter) : it(iter) {} 244 }; 245 246 } 247 248 inline void delaunay() { 151 249 Counter cnt("Number of edges added: "); 152 std::vector<Pedge> pedges; 153 for(NodeIt n(g);n!=INVALID;++n) 154 for(NodeIt m=++(NodeIt(n));m!=INVALID;++m) 250 251 using namespace _delaunay_bits; 252 253 typedef _delaunay_bits::Part Part; 254 typedef std::vector<std::pair<double, int> > SiteHeap; 255 256 257 std::vector<Point> points; 258 std::vector<Node> nodes; 259 260 for (NodeIt it(g); it != INVALID; ++it) { 261 nodes.push_back(it); 262 points.push_back(coords[it]); 263 } 264 265 SiteHeap siteheap(points.size()); 266 267 double sweep; 268 269 270 for (int i = 0; i < int(siteheap.size()); ++i) { 271 siteheap[i] = std::make_pair(points[i].x, i); 272 } 273 274 std::sort(siteheap.begin(), siteheap.end()); 275 sweep = siteheap.front().first; 276 277 YLess yless(points, sweep); 278 Beach beach(yless); 279 280 SpikeHeap spikeheap; 281 282 std::set<std::pair<int, int> > edges; 283 284 beach.insert(std::make_pair(Part(-1, siteheap[0].second, -1), 285 spikeheap.end())); 286 int siteindex = 1; 287 288 while (siteindex < int(points.size()) || !spikeheap.empty()) { 289 290 SpikeHeap::iterator spit = spikeheap.begin(); 291 292 if (siteindex < int(points.size()) && 293 (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) { 294 int site = siteheap[siteindex].second; 295 sweep = siteheap[siteindex].first; 296 297 Beach::iterator bit = beach.upper_bound(Part(site, site, site)); 298 299 if (bit->second != spikeheap.end()) { 300 spikeheap.erase(bit->second); 301 } 302 303 int prev = bit->first.prev; 304 int curr = bit->first.curr; 305 int next = bit->first.next; 306 307 beach.erase(bit); 308 309 SpikeHeap::iterator pit = spikeheap.end(); 310 if (prev != -1 && 311 circle_form(points[prev], points[curr], points[site])) { 312 double x = circle_point(points[prev], points[curr], points[site]); 313 pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); 314 pit->second.it = 315 beach.insert(std::make_pair(Part(prev, curr, site), pit)); 316 } else { 317 beach.insert(std::make_pair(Part(prev, curr, site), pit)); 318 } 319 320 beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end())); 321 322 SpikeHeap::iterator nit = spikeheap.end(); 323 if (next != -1 && 324 circle_form(points[site], points[curr],points[next])) { 325 double x = circle_point(points[site], points[curr], points[next]); 326 nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); 327 nit->second.it = 328 beach.insert(std::make_pair(Part(site, curr, next), nit)); 329 } else { 330 beach.insert(std::make_pair(Part(site, curr, next), nit)); 331 } 332 333 ++siteindex; 334 } else { 335 sweep = spit->first; 336 337 Beach::iterator bit = spit->second.it; 338 339 int prev = bit->first.prev; 340 int curr = bit->first.curr; 341 int next = bit->first.next; 342 155 343 { 156 Pedge p; 157 p.a=n; 158 p.b=m; 159 p.len=(coords[m]-coords[n]).normSquare(); 160 pedges.push_back(p); 161 } 162 std::sort(pedges.begin(),pedges.end(),pedgeLess); 163 for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi) 164 { 165 Line li(pi->a,pi->b); 166 UEdgeIt e(g); 167 for(;e!=INVALID && !cross(e,li);++e) ; 168 UEdge ne; 169 if(e==INVALID) { 170 ne=g.addEdge(pi->a,pi->b); 171 edges.push_back(ne); 172 cnt++; 173 } 174 } 344 std::pair<int, int> edge; 345 346 edge = prev < curr ? 347 std::make_pair(prev, curr) : std::make_pair(curr, prev); 348 349 if (edges.find(edge) == edges.end()) { 350 edges.insert(edge); 351 g.addEdge(nodes[prev], nodes[curr]); 352 ++cnt; 353 } 354 355 edge = curr < next ? 356 std::make_pair(curr, next) : std::make_pair(next, curr); 357 358 if (edges.find(edge) == edges.end()) { 359 edges.insert(edge); 360 g.addEdge(nodes[curr], nodes[next]); 361 ++cnt; 362 } 363 } 364 365 Beach::iterator pbit = bit; --pbit; 366 int ppv = pbit->first.prev; 367 Beach::iterator nbit = bit; ++nbit; 368 int nnt = nbit->first.next; 369 370 if (bit->second != spikeheap.end()) spikeheap.erase(bit->second); 371 if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second); 372 if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second); 373 374 beach.erase(nbit); 375 beach.erase(bit); 376 beach.erase(pbit); 377 378 SpikeHeap::iterator pit = spikeheap.end(); 379 if (ppv != -1 && ppv != next && 380 circle_form(points[ppv], points[prev], points[next])) { 381 double x = circle_point(points[ppv], points[prev], points[next]); 382 if (x < sweep) x = sweep; 383 pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); 384 pit->second.it = 385 beach.insert(std::make_pair(Part(ppv, prev, next), pit)); 386 } else { 387 beach.insert(std::make_pair(Part(ppv, prev, next), pit)); 388 } 389 390 SpikeHeap::iterator nit = spikeheap.end(); 391 if (nnt != -1 && prev != nnt && 392 circle_form(points[prev], points[next], points[nnt])) { 393 double x = circle_point(points[prev], points[next], points[nnt]); 394 if (x < sweep) x = sweep; 395 nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); 396 nit->second.it = 397 beach.insert(std::make_pair(Part(prev, next, nnt), nit)); 398 } else { 399 beach.insert(std::make_pair(Part(prev, next, nnt), nit)); 400 } 401 402 } 403 } 404 405 for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) { 406 int curr = it->first.curr; 407 int next = it->first.next; 408 409 if (next == -1) continue; 410 411 std::pair<int, int> edge; 412 413 edge = curr < next ? 414 std::make_pair(curr, next) : std::make_pair(next, curr); 415 416 if (edges.find(edge) == edges.end()) { 417 edges.insert(edge); 418 g.addEdge(nodes[curr], nodes[next]); 419 ++cnt; 420 } 421 } 175 422 } 176 423 … … 246 493 } 247 494 495 template <typename UGraph, typename CoordMap> 496 class LengthSquareMap { 497 public: 498 typedef typename UGraph::UEdge Key; 499 typedef typename CoordMap::Value::Value Value; 500 501 LengthSquareMap(const UGraph& ugraph, const CoordMap& coords) 502 : _ugraph(ugraph), _coords(coords) {} 503 504 Value operator[](const Key& key) const { 505 return (_coords[_ugraph.target(key)] - 506 _coords[_ugraph.source(key)]).normSquare(); 507 } 508 509 private: 510 511 const UGraph& _ugraph; 512 const CoordMap& _coords; 513 }; 514 248 515 void minTree() { 249 int en=0;250 int pr=0;251 516 std::vector<Pedge> pedges; 252 517 Timer T; 253 std::cout << T.realTime() << "s: Setting up the edges...\n"; 254 for(NodeIt n(g);n!=INVALID;++n) 255 { 256 for(NodeIt m=++(NodeIt(n));m!=INVALID;++m) 257 { 258 Pedge p; 259 p.a=n; 260 p.b=m; 261 p.len=(coords[m]-coords[n]).normSquare(); 262 pedges.push_back(p); 263 } 264 if(progress && en>=pr*double(N)/100) 265 { 266 std::cout << pr << "% \r" << std::flush; 267 pr++; 268 } 269 } 270 std::cout << T.realTime() << "s: Sorting the edges...\n"; 271 std::sort(pedges.begin(),pedges.end(),pedgeLess); 272 std::cout << T.realTime() << "s: Creating the tree...\n"; 273 ListUGraph::NodeMap<int> comp(g); 274 UnionFind<ListUGraph::NodeMap<int> > uf(comp); 275 for (NodeIt it(g); it != INVALID; ++it) 276 uf.insert(it); 277 for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi) 278 { 279 if ( uf.join(pi->a,pi->b) ) { 280 g.addEdge(pi->a,pi->b); 281 en++; 282 if(en>=N-1) 283 break; 284 } 285 } 518 std::cout << T.realTime() << "s: Creating delaunay triangulation...\n"; 519 delaunay(); 520 std::cout << T.realTime() << "s: Calculating spanning tree...\n"; 521 LengthSquareMap<ListUGraph, ListUGraph::NodeMap<Point> > ls(g, coords); 522 ListUGraph::UEdgeMap<bool> tree(g); 523 kruskal(g, ls, tree); 524 std::cout << T.realTime() << "s: Removing non tree edges...\n"; 525 std::vector<UEdge> remove; 526 for (UEdgeIt e(g); e != INVALID; ++e) { 527 if (!tree[e]) remove.push_back(e); 528 } 529 for(int i = 0; i < int(remove.size()); ++i) { 530 g.erase(remove[i]); 531 } 286 532 std::cout << T.realTime() << "s: Done\n"; 287 533 } … … 409 655 .boolOption("tsp2", "Create a TSP tour (tree based)") 410 656 .optionGroup("alg","tsp2") 657 .boolOption("dela", "Delaunay triangulation graph") 658 .optionGroup("alg","dela") 411 659 .onlyOneGroup("alg") 660 .boolOption("rand", "Use time seed for random number generator") 661 .optionGroup("rand", "rand") 662 .intOption("seed", "Random seed", -1) 663 .optionGroup("rand", "seed") 664 .onlyOneGroup("rand") 412 665 .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'") 413 666 .run(); 667 668 if (ap["rand"]) { 669 int seed = time(0); 670 std::cout << "Random number seed: " << seed << std::endl; 671 rnd = Random(seed); 672 } 673 if (ap.given("seed")) { 674 int seed = ap["seed"]; 675 std::cout << "Random number seed: " << seed << std::endl; 676 rnd = Random(seed); 677 } 414 678 415 679 std::string prefix; … … 464 728 } 465 729 } 730 731 // for (ListUGraph::NodeIt n(g); n != INVALID; ++n) { 732 // std::cerr << coords[n] << std::endl; 733 // } 466 734 467 735 if(ap["tsp"]) { … … 482 750 else if(ap["tree"]) { 483 751 minTree(); 752 } 753 else if(ap["dela"]) { 754 delaunay(); 484 755 } 485 756
Note: See TracChangeset
for help on using the changeset viewer.