1 /* 2 Copyright (c) 2005-2021 Intel Corporation 3 4 Licensed under the Apache License, Version 2.0 (the "License"); 5 you may not use this file except in compliance with the License. 6 You may obtain a copy of the License at 7 8 http://www.apache.org/licenses/LICENSE-2.0 9 10 Unless required by applicable law or agreed to in writing, software 11 distributed under the License is distributed on an "AS IS" BASIS, 12 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 See the License for the specific language governing permissions and 14 limitations under the License. 15 */ 16 17 /* 18 This file contains a few implementations, so it may look overly complicated. 19 The most efficient implementation is also separated into convex_hull_sample.cpp 20 */ 21 #include "convex_hull.hpp" 22 23 typedef util::point<double> point_t; 24 25 #ifndef USETBB 26 #define USETBB 1 27 #endif 28 #ifndef USECONCVEC 29 #define USECONCVEC 1 30 #endif 31 32 #if !USETBB // Serial implementation of Quick Hull algorithm 33 34 typedef std::vector<point_t> pointVec_t; 35 36 void serial_initialize(pointVec_t &points); 37 38 class FindXExtremum : public std::unary_function<const point_t &, void> { 39 public: 40 typedef enum { minX, maxX } extremumType; 41 42 FindXExtremum(const point_t &frstPoint, extremumType exType_) 43 : extrXPoint(frstPoint), 44 exType(exType_) {} 45 46 void operator()(const point_t &p) { 47 if (closerToExtremum(p)) 48 extrXPoint = p; 49 } 50 51 operator point_t() { 52 return extrXPoint; 53 } 54 55 private: 56 const extremumType exType; 57 point_t extrXPoint; 58 59 bool closerToExtremum(const point_t &p) const { 60 switch (exType) { 61 case minX: return p.x < extrXPoint.x; break; 62 case maxX: return p.x > extrXPoint.x; break; 63 } 64 return false; // avoid warning 65 } 66 }; 67 68 template <FindXExtremum::extremumType type> 69 point_t extremum(const pointVec_t &points) { 70 assert(!points.empty()); 71 return std::for_each(points.begin(), points.end(), FindXExtremum(points[0], type)); 72 } 73 74 class SplitByCP : public std::unary_function<const point_t &, void> { 75 pointVec_t &reducedSet; 76 point_t p1, p2; 77 point_t farPoint; 78 double howFar; 79 80 public: 81 SplitByCP(point_t _p1, point_t _p2, pointVec_t &_reducedSet) 82 : p1(_p1), 83 p2(_p2), 84 reducedSet(_reducedSet), 85 howFar(0), 86 farPoint(p1) {} 87 88 void operator()(const point_t &p) { 89 double cp; 90 if ((p != p1) && (p != p2)) { 91 cp = util::cross_product(p1, p2, p); 92 if (cp > 0) { 93 reducedSet.push_back(p); 94 if (cp > howFar) { 95 farPoint = p; 96 howFar = cp; 97 } 98 } 99 } 100 } 101 102 operator point_t() { 103 return farPoint; 104 } 105 }; 106 107 point_t divide(const pointVec_t &P, pointVec_t &P_reduced, const point_t &p1, const point_t &p2) { 108 SplitByCP splitByCP(p1, p2, P_reduced); 109 point_t farPoint = std::for_each(P.begin(), P.end(), splitByCP); 110 111 if (util::verbose) { 112 std::stringstream ss; 113 ss << P.size() << " nodes in bucket" 114 << ", " 115 << "dividing by: [ " << p1 << ", " << p2 << " ], " 116 << "farthest node: " << farPoint; 117 util::OUTPUT.push_back(ss.str()); 118 } 119 120 return farPoint; 121 } 122 123 void divide_and_conquer(const pointVec_t &P, pointVec_t &H, point_t p1, point_t p2) { 124 assert(P.size() >= 2); 125 pointVec_t P_reduced; 126 pointVec_t H1, H2; 127 point_t p_far = divide(P, P_reduced, p1, p2); 128 if (P_reduced.size() < 2) { 129 H.push_back(p1); 130 H.insert(H.end(), P_reduced.begin(), P_reduced.end()); 131 } 132 else { 133 divide_and_conquer(P_reduced, H1, p1, p_far); 134 divide_and_conquer(P_reduced, H2, p_far, p2); 135 136 H.insert(H.end(), H1.begin(), H1.end()); 137 H.insert(H.end(), H2.begin(), H2.end()); 138 } 139 } 140 141 void quickhull(const pointVec_t &points, pointVec_t &hull) { 142 if (points.size() < 2) { 143 hull.insert(hull.end(), points.begin(), points.end()); 144 return; 145 } 146 point_t p_maxx = extremum<FindXExtremum::maxX>(points); 147 point_t p_minx = extremum<FindXExtremum::minX>(points); 148 149 pointVec_t H; 150 151 divide_and_conquer(points, hull, p_maxx, p_minx); 152 divide_and_conquer(points, H, p_minx, p_maxx); 153 hull.insert(hull.end(), H.begin(), H.end()); 154 } 155 156 int main(int argc, char *argv[]) { 157 utility::thread_number_range single_thread([] { 158 return -1; 159 }); 160 util::ParseInputArgs(argc, argv, single_thread); 161 162 pointVec_t points; 163 pointVec_t hull; 164 util::my_time_t tm_init, tm_start, tm_end; 165 166 std::cout << "Starting serial version of QUICK HULL algorithm" 167 << "\n"; 168 169 tm_init = util::gettime(); 170 serial_initialize(points); 171 tm_start = util::gettime(); 172 std::cout << "Init time: " << util::time_diff(tm_init, tm_start) 173 << " Points in input: " << points.size() << "\n"; 174 tm_start = util::gettime(); 175 quickhull(points, hull); 176 tm_end = util::gettime(); 177 std::cout << "Serial time: " << util::time_diff(tm_start, tm_end) 178 << " Points in hull: " << hull.size() << "\n"; 179 180 return 0; 181 } 182 183 #else // USETBB - parallel version of Quick Hull algorithm 184 185 #include "oneapi/tbb/parallel_for.h" 186 #include "oneapi/tbb/parallel_reduce.h" 187 #include "oneapi/tbb/blocked_range.h" 188 189 typedef oneapi::tbb::blocked_range<std::size_t> range_t; 190 191 #if USECONCVEC 192 #include "oneapi/tbb/concurrent_vector.h" 193 194 typedef oneapi::tbb::concurrent_vector<point_t> pointVec_t; 195 196 void appendVector(const point_t *src, std::size_t srcSize, pointVec_t &dest) { 197 std::copy(src, src + srcSize, dest.grow_by(srcSize)); 198 } 199 200 void appendVector(const pointVec_t &src, pointVec_t &dest) { 201 std::copy(src.begin(), src.end(), dest.grow_by(src.size())); 202 } 203 204 void grow_vector_to_at_least(pointVec_t &vect, std::size_t size) { 205 vect.grow_to_at_least(size); 206 } 207 #else // USE STD::VECTOR - include spin_mutex.h and lock vector operations 208 #include "oneapi/tbb/spin_mutex.h" 209 210 typedef oneapi::tbb::spin_mutex mutex_t; 211 typedef std::vector<point_t> pointVec_t; 212 213 void appendVector(mutex_t &insertMutex, const pointVec_t &src, pointVec_t &dest) { 214 mutex_t::scoped_lock lock(insertMutex); 215 dest.insert(dest.end(), src.begin(), src.end()); 216 } 217 218 void appendVector(mutex_t &insertMutex, const point_t *src, std::size_t srcSize, pointVec_t &dest) { 219 mutex_t::scoped_lock lock(insertMutex); 220 dest.insert(dest.end(), src, src + srcSize); 221 } 222 223 void grow_vector_to_at_least(mutex_t &mutex, pointVec_t &vect, std::size_t size) { 224 mutex_t::scoped_lock lock(mutex); 225 if (vect.size() < size) { 226 vect.resize(size); 227 } 228 } 229 #endif // USECONCVEC 230 231 class FillRNDPointsVector { 232 pointVec_t &points; 233 234 public: 235 static const std::size_t grainSize = cfg::generateGrainSize; 236 #if !USECONCVEC 237 static mutex_t pushBackMutex; 238 #endif // USECONCVEC 239 240 explicit FillRNDPointsVector(pointVec_t &_points) : points(_points) {} 241 242 void operator()(const range_t &range) const { 243 util::rng the_rng(range.begin()); 244 const std::size_t i_end = range.end(); 245 std::size_t count = 0; 246 #if USECONCVEC 247 points.grow_to_at_least(i_end); 248 #else // Locked enlarge to a not thread-safe STD::VECTOR 249 grow_vector_to_at_least(pushBackMutex, points, i_end); 250 #endif // USECONCVEC 251 252 for (std::size_t i = range.begin(); i != i_end; ++i) { 253 points[i] = util::GenerateRNDPoint<double>(count, the_rng, util::rng::max_rand); 254 } 255 } 256 }; 257 258 class FillRNDPointsVector_buf { 259 pointVec_t &points; 260 261 public: 262 static const std::size_t grainSize = cfg::generateGrainSize; 263 #if !USECONCVEC 264 static mutex_t insertMutex; 265 #endif // USECONCVEC 266 267 explicit FillRNDPointsVector_buf(pointVec_t &_points) : points(_points) {} 268 269 void operator()(const range_t &range) const { 270 util::rng the_rng(range.begin()); 271 const std::size_t i_end = range.end(); 272 std::size_t count = 0, j = 0; 273 point_t tmp_vec[grainSize]; 274 275 for (std::size_t i = range.begin(); i != i_end; ++i) { 276 tmp_vec[j++] = util::GenerateRNDPoint<double>(count, the_rng, util::rng::max_rand); 277 } 278 #if USECONCVEC 279 grow_vector_to_at_least(points, range.end()); 280 #else // USE STD::VECTOR 281 grow_vector_to_at_least(insertMutex, points, range.end()); 282 #endif // USECONCVEC 283 std::copy(tmp_vec, tmp_vec + j, points.begin() + range.begin()); 284 } 285 }; 286 287 #if !USECONCVEC 288 mutex_t FillRNDPointsVector::pushBackMutex{}; 289 mutex_t FillRNDPointsVector_buf::insertMutex{}; 290 #endif 291 292 template <typename BodyType> 293 void initialize(pointVec_t &points) { 294 //This function generate the same series of point on every call. 295 //Reproducibility is needed for benchmarking to produce reliable results. 296 //It is achieved through the following points: 297 // - FillRNDPointsVector_buf instance has its own local instance 298 // of random number generator, which in turn does not use any global data 299 // - oneapi::tbb::simple_partitioner produce the same set of ranges on every call to 300 // oneapi::tbb::parallel_for 301 // - local RNG instances are seeded by the starting indexes of corresponding ranges 302 // - grow_to_at_least() enables putting points into the resulting vector in deterministic order 303 // (unlike concurrent push_back or grow_by). 304 305 // In the buffered version, a temporary storage for as much as grainSize elements 306 // is allocated inside the body. Since auto_partitioner may increase effective 307 // range size which would cause a crash, simple partitioner has to be used. 308 309 oneapi::tbb::parallel_for(range_t(0, cfg::numberOfPoints, BodyType::grainSize), 310 BodyType(points), 311 oneapi::tbb::simple_partitioner()); 312 } 313 314 class FindXExtremum { 315 public: 316 typedef enum { minX, maxX } extremumType; 317 318 static const std::size_t grainSize = cfg::findExtremumGrainSize; 319 320 FindXExtremum(const pointVec_t &points_, extremumType exType_) 321 : points(points_), 322 exType(exType_), 323 extrXPoint(points[0]) {} 324 325 FindXExtremum(const FindXExtremum &fxex, oneapi::tbb::split) 326 : points(fxex.points), 327 exType(fxex.exType), 328 extrXPoint(fxex.extrXPoint) {} 329 330 void operator()(const range_t &range) { 331 const std::size_t i_end = range.end(); 332 if (!range.empty()) { 333 for (std::size_t i = range.begin(); i != i_end; ++i) { 334 if (closerToExtremum(points[i])) { 335 extrXPoint = points[i]; 336 } 337 } 338 } 339 } 340 341 void join(const FindXExtremum &rhs) { 342 if (closerToExtremum(rhs.extrXPoint)) { 343 extrXPoint = rhs.extrXPoint; 344 } 345 } 346 347 point_t extremeXPoint() { 348 return extrXPoint; 349 } 350 351 private: 352 const pointVec_t &points; 353 const extremumType exType; 354 point_t extrXPoint; 355 bool closerToExtremum(const point_t &p) const { 356 switch (exType) { 357 case minX: return p.x < extrXPoint.x; break; 358 case maxX: return p.x > extrXPoint.x; break; 359 } 360 return false; // avoid warning 361 } 362 }; 363 364 template <FindXExtremum::extremumType type> 365 point_t extremum(const pointVec_t &P) { 366 FindXExtremum fxBody(P, type); 367 oneapi::tbb::parallel_reduce(range_t(0, P.size(), FindXExtremum::grainSize), fxBody); 368 return fxBody.extremeXPoint(); 369 } 370 371 class SplitByCP { 372 const pointVec_t &initialSet; 373 pointVec_t &reducedSet; 374 point_t p1, p2; 375 point_t farPoint; 376 double howFar; 377 378 public: 379 static const std::size_t grainSize = cfg::divideGrainSize; 380 #if !USECONCVEC 381 static mutex_t pushBackMutex; 382 #endif // USECONCVEC 383 384 SplitByCP(point_t _p1, point_t _p2, const pointVec_t &_initialSet, pointVec_t &_reducedSet) 385 : p1(_p1), 386 p2(_p2), 387 initialSet(_initialSet), 388 reducedSet(_reducedSet), 389 howFar(0), 390 farPoint(p1) {} 391 392 SplitByCP(SplitByCP &sbcp, oneapi::tbb::split) 393 : p1(sbcp.p1), 394 p2(sbcp.p2), 395 initialSet(sbcp.initialSet), 396 reducedSet(sbcp.reducedSet), 397 howFar(0), 398 farPoint(p1) {} 399 400 void operator()(const range_t &range) { 401 const std::size_t i_end = range.end(); 402 double cp; 403 for (std::size_t i = range.begin(); i != i_end; ++i) { 404 if ((initialSet[i] != p1) && (initialSet[i] != p2)) { 405 cp = util::cross_product(p1, p2, initialSet[i]); 406 if (cp > 0) { 407 #if USECONCVEC 408 reducedSet.push_back(initialSet[i]); 409 #else // Locked push_back to a not thread-safe STD::VECTOR 410 { 411 mutex_t::scoped_lock lock(pushBackMutex); 412 reducedSet.push_back(initialSet[i]); 413 } 414 #endif // USECONCVEC 415 if (cp > howFar) { 416 farPoint = initialSet[i]; 417 howFar = cp; 418 } 419 } 420 } 421 } 422 } 423 424 void join(const SplitByCP &rhs) { 425 if (rhs.howFar > howFar) { 426 howFar = rhs.howFar; 427 farPoint = rhs.farPoint; 428 } 429 } 430 431 point_t farthestPoint() const { 432 return farPoint; 433 } 434 }; 435 436 class SplitByCP_buf { 437 const pointVec_t &initialSet; 438 pointVec_t &reducedSet; 439 point_t p1, p2; 440 point_t farPoint; 441 double howFar; 442 443 public: 444 static const std::size_t grainSize = cfg::divideGrainSize; 445 #if !USECONCVEC 446 static mutex_t insertMutex; 447 #endif // USECONCVEC 448 449 SplitByCP_buf(point_t _p1, point_t _p2, const pointVec_t &_initialSet, pointVec_t &_reducedSet) 450 : p1(_p1), 451 p2(_p2), 452 initialSet(_initialSet), 453 reducedSet(_reducedSet), 454 howFar(0), 455 farPoint(p1) {} 456 457 SplitByCP_buf(SplitByCP_buf &sbcp, oneapi::tbb::split) 458 : p1(sbcp.p1), 459 p2(sbcp.p2), 460 initialSet(sbcp.initialSet), 461 reducedSet(sbcp.reducedSet), 462 howFar(0), 463 farPoint(p1) {} 464 465 void operator()(const range_t &range) { 466 const std::size_t i_end = range.end(); 467 std::size_t j = 0; 468 double cp; 469 point_t tmp_vec[grainSize]; 470 for (std::size_t i = range.begin(); i != i_end; ++i) { 471 if ((initialSet[i] != p1) && (initialSet[i] != p2)) { 472 cp = util::cross_product(p1, p2, initialSet[i]); 473 if (cp > 0) { 474 tmp_vec[j++] = initialSet[i]; 475 if (cp > howFar) { 476 farPoint = initialSet[i]; 477 howFar = cp; 478 } 479 } 480 } 481 } 482 483 #if USECONCVEC 484 appendVector(tmp_vec, j, reducedSet); 485 #else // USE STD::VECTOR 486 appendVector(insertMutex, tmp_vec, j, reducedSet); 487 #endif // USECONCVEC 488 } 489 490 void join(const SplitByCP_buf &rhs) { 491 if (rhs.howFar > howFar) { 492 howFar = rhs.howFar; 493 farPoint = rhs.farPoint; 494 } 495 } 496 497 point_t farthestPoint() const { 498 return farPoint; 499 } 500 }; 501 502 #if !USECONCVEC 503 mutex_t SplitByCP::pushBackMutex{}; 504 mutex_t SplitByCP_buf::insertMutex{}; 505 #endif 506 507 template <typename BodyType> 508 point_t divide(const pointVec_t &P, pointVec_t &P_reduced, const point_t &p1, const point_t &p2) { 509 BodyType body(p1, p2, P, P_reduced); 510 // Must use simple_partitioner (see the comment in initialize() above) 511 oneapi::tbb::parallel_reduce( 512 range_t(0, P.size(), BodyType::grainSize), body, oneapi::tbb::simple_partitioner()); 513 514 if (util::verbose) { 515 std::stringstream ss; 516 ss << P.size() << " nodes in bucket" 517 << ", " 518 << "dividing by: [ " << p1 << ", " << p2 << " ], " 519 << "farthest node: " << body.farthestPoint(); 520 util::OUTPUT.push_back(ss.str()); 521 } 522 523 return body.farthestPoint(); 524 } 525 526 void divide_and_conquer(const pointVec_t &P, pointVec_t &H, point_t p1, point_t p2, bool buffered) { 527 assert(P.size() >= 2); 528 pointVec_t P_reduced; 529 pointVec_t H1, H2; 530 point_t p_far; 531 532 if (buffered) { 533 p_far = divide<SplitByCP_buf>(P, P_reduced, p1, p2); 534 } 535 else { 536 p_far = divide<SplitByCP>(P, P_reduced, p1, p2); 537 } 538 539 if (P_reduced.size() < 2) { 540 H.push_back(p1); 541 #if USECONCVEC 542 appendVector(P_reduced, H); 543 #else // insert into STD::VECTOR 544 H.insert(H.end(), P_reduced.begin(), P_reduced.end()); 545 #endif 546 } 547 else { 548 divide_and_conquer(P_reduced, H1, p1, p_far, buffered); 549 divide_and_conquer(P_reduced, H2, p_far, p2, buffered); 550 551 #if USECONCVEC 552 appendVector(H1, H); 553 appendVector(H2, H); 554 #else // insert into STD::VECTOR 555 H.insert(H.end(), H1.begin(), H1.end()); 556 H.insert(H.end(), H2.begin(), H2.end()); 557 #endif 558 } 559 } 560 561 void quickhull(const pointVec_t &points, pointVec_t &hull, bool buffered) { 562 if (points.size() < 2) { 563 #if USECONCVEC 564 appendVector(points, hull); 565 #else // STD::VECTOR 566 hull.insert(hull.end(), points.begin(), points.end()); 567 #endif // USECONCVEC 568 return; 569 } 570 571 point_t p_maxx = extremum<FindXExtremum::maxX>(points); 572 point_t p_minx = extremum<FindXExtremum::minX>(points); 573 574 pointVec_t H; 575 576 divide_and_conquer(points, hull, p_maxx, p_minx, buffered); 577 divide_and_conquer(points, H, p_minx, p_maxx, buffered); 578 #if USECONCVEC 579 appendVector(H, hull); 580 #else // STD::VECTOR 581 hull.insert(hull.end(), H.begin(), H.end()); 582 #endif // USECONCVEC 583 } 584 585 int main(int argc, char *argv[]) { 586 utility::thread_number_range threads(utility::get_default_num_threads); 587 util::ParseInputArgs(argc, argv, threads); 588 589 int nthreads; 590 util::my_time_t tm_init, tm_start, tm_end; 591 592 #if USECONCVEC 593 std::cout << "Starting TBB unbuffered push_back version of QUICK HULL algorithm" 594 << "\n"; 595 #else 596 std::cout << "Starting STL locked unbuffered push_back version of QUICK HULL algorithm" 597 << "\n"; 598 #endif // USECONCVEC 599 600 for (nthreads = threads.first; nthreads <= threads.last; nthreads = threads.step(nthreads)) { 601 pointVec_t points; 602 pointVec_t hull; 603 604 oneapi::tbb::global_control c(oneapi::tbb::global_control::max_allowed_parallelism, 605 nthreads); 606 tm_init = util::gettime(); 607 initialize<FillRNDPointsVector>(points); 608 tm_start = util::gettime(); 609 std::cout << "Parallel init time on " << nthreads 610 << " threads: " << util::time_diff(tm_init, tm_start) 611 << " Points in input: " << points.size() << "\n"; 612 613 tm_start = util::gettime(); 614 quickhull(points, hull, false); 615 tm_end = util::gettime(); 616 std::cout << "Time on " << nthreads << " threads: " << util::time_diff(tm_start, tm_end) 617 << " Points in hull: " << hull.size() << "\n"; 618 } 619 620 #if USECONCVEC 621 std::cout << "Starting TBB buffered version of QUICK HULL algorithm" 622 << "\n"; 623 #else 624 std::cout << "Starting STL locked buffered version of QUICK HULL algorithm" 625 << "\n"; 626 #endif 627 628 for (nthreads = threads.first; nthreads <= threads.last; nthreads = threads.step(nthreads)) { 629 pointVec_t points; 630 pointVec_t hull; 631 632 oneapi::tbb::global_control c(oneapi::tbb::global_control::max_allowed_parallelism, 633 nthreads); 634 635 tm_init = util::gettime(); 636 initialize<FillRNDPointsVector_buf>(points); 637 tm_start = util::gettime(); 638 std::cout << "Init time on " << nthreads 639 << " threads: " << util::time_diff(tm_init, tm_start) 640 << " Points in input: " << points.size() << "\n"; 641 642 tm_start = util::gettime(); 643 quickhull(points, hull, true); 644 tm_end = util::gettime(); 645 std::cout << "Time on " << nthreads << " threads: " << util::time_diff(tm_start, tm_end) 646 << " Points in hull: " << hull.size() << "\n"; 647 } 648 649 return 0; 650 } 651 652 #endif // USETBB 653 654 void serial_initialize(pointVec_t &points) { 655 points.reserve(cfg::numberOfPoints); 656 657 unsigned int rseed = 1; 658 for (std::size_t i = 0, count = 0; long(i) < cfg::numberOfPoints; ++i) { 659 points.push_back(util::GenerateRNDPoint<double>(count, &std::rand, RAND_MAX)); 660 } 661 } 662