Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <climits>
- #include <cmath>
- #include <algorithm>
- #include <chrono>
- #include <iostream>
- #include <numeric>
- #include <random>
- #include <utility>
- template<typename CT>
- struct Point {
- typedef CT CoordinateType;
- CoordinateType x;
- CoordinateType y;
- };
- template<typename CT>
- std::ostream &operator<<(std::ostream &os, const Point<CT> &p)
- {
- return os << "x:" << p.x << " y:" << p.y;
- }
- template<typename T1, typename T2>
- std::ostream &operator<<(std::ostream &os, const std::pair<T1, T2> &pair)
- {
- return os << " (" << pair.first << "," << pair.second << ") ";
- }
- template<typename CT>
- CT squareRoot(CT n) {
- return static_cast<CT>(std::pow(n, 0.5));
- }
- template<typename CT>
- CT distanceSquared(const Point<CT> &first, const Point<CT> &second) {
- const CT dx = (first.x - second.x), dy = (first.y - second.y);
- return (dx*dx) + (dy*dy);
- }
- template<typename CT>
- CT distance(const Point<CT> &first, const Point<CT> &second) {
- return static_cast<CT>(std::pow(distanceSquared(first, second), 0.5));
- }
- template<typename PointType, typename PairType = std::pair<PointType, PointType>, typename CT = typename PointType::CoordinateType >
- std::pair<PairType, CT> exhaustiveSearch(const PointType p[], int begin, int end)
- {
- PairType pair;
- CT minSquared = std::numeric_limits<CT>::max();
- for (int i = begin; i < (end - 1); ++i) {
- const PointType &first(p[i]);
- for (int j = i + 1; j < end; ++j) {
- const PointType &second(p[j]);
- const CT distSquared = distanceSquared(second, first);
- if (distSquared < minSquared) {
- minSquared = distSquared;
- pair = std::make_pair(first, second);
- }
- }
- }
- return std::make_pair(pair, squareRoot(minSquared));
- }
- template<typename PointType, typename PairType = std::pair<PointType, PointType>, typename CT = typename PointType::CoordinateType >
- std::pair<PairType, CT> recursiveSearch(const PointType p[], int begin, int end)
- {
- const int span = end - begin;
- if (span < 64) {
- return exhaustiveSearch(p, begin, end);
- }
- const int mid = begin + (span + 1)/2;
- std::pair<PairType, CT> lresult = recursiveSearch(p, begin, mid);
- std::pair<PairType, CT> rresult = recursiveSearch(p, mid, end);
- PairType &lpair = lresult.first;
- CT &lmin = lresult.second;
- PairType &rpair = rresult.first;
- CT &rmin = rresult.second;
- {
- CT min = std::min(lmin, rmin);
- CT minSquared = min * min;
- CT upperBound = p[mid - 1].x + min;
- for (int i = mid; i < end; ++i) {
- const PointType &rp = p[i];
- if (rp.x >= upperBound) {
- break;
- } else {
- CT lowerBound = rp.x - min;
- for (int j = mid - 1; j >= begin; --j) {
- const PointType &lp = p[j];
- if (lp.x <= lowerBound) {
- break;
- } else {
- const CT distSquared = distanceSquared(lp, rp);
- if (distSquared < minSquared) {
- minSquared = distSquared;
- min = squareRoot(distSquared);
- lmin = min;
- lpair = std::make_pair(lp, rp);
- upperBound = p[mid - 1].x + min;
- lowerBound = rp.x - min;
- }
- }
- }
- }
- }
- }
- {
- CT min = std::min(lmin, rmin);
- CT minSquared = min * min;
- CT lowerBound = p[mid].x - min;
- for (int i = mid - 1; i >= begin; --i) {
- const PointType &lp = p[i];
- if (lp.x <= lowerBound) {
- break;
- } else {
- CT upperBound = lp.x + min;
- for (int j = mid; j < end; ++j) {
- const PointType &rp = p[j];
- if (rp.x >= upperBound) {
- break;
- } else {
- const CT distSquared = distanceSquared(lp, rp);
- if (distSquared < minSquared) {
- minSquared = distSquared;
- min = squareRoot(distSquared);
- rmin = min;
- rpair = std::make_pair(lp, rp);
- lowerBound = p[mid].x - min;
- upperBound = lp.x + min;
- }
- }
- }
- }
- }
- }
- return std::make_pair((rmin < lmin ? rpair : lpair), (rmin < lmin ? rmin : lmin));
- }
- template<typename PointType, typename PairType = std::pair<PointType, PointType>, typename CT = typename PointType::CoordinateType >
- std::pair<PairType, CT> neighbourSearch(const PointType p[], int begin, int end)
- {
- PairType pair;
- CT minSquared = distanceSquared(p[begin], p[begin + 1]);
- CT min = squareRoot(minSquared);
- for (int i = begin; i < (end - 1); ++i) {
- const PointType &first = p[i];
- CT upperBound = first.x + min;
- for (int j = i + 1; j < end; ++j) {
- const PointType &second = p[j];
- if (second.x >= upperBound) {
- break;
- } else {
- const CT distSquared = distanceSquared(first, second);
- if (distSquared < minSquared) {
- minSquared = distSquared;
- min = squareRoot(distSquared);
- pair = std::make_pair(first, second);
- upperBound = first.x + min;
- }
- }
- }
- }
- return std::make_pair(pair, min);
- }
- template<typename CT>
- void timeOperations()
- {
- #ifdef USE_EXHAUSTIVE
- const int N = 40 * 1000;
- #else
- const int N = 400 * 1000;
- #endif
- std::vector<Point<CT> > points(N);
- std::vector<Point<CT> > sorted(N);
- std::default_random_engine rng(std::chrono::system_clock::to_time_t(std::chrono::system_clock::now()));
- std::uniform_real_distribution<CT> dst(0.0, 1.0);
- for (int i = 0; i < N; ++i) { points[i] = Point<CT>{ dst(rng), dst(rng) }; }
- std::pair<std::pair<Point<CT>, Point<CT> >, CT> min;
- std::chrono::system_clock::time_point tp = std::chrono::system_clock::now();
- std::chrono::system_clock::duration elapsed;
- #ifdef USE_EXHAUSTIVE
- min = exhaustiveSearch(points.data(), 0, N);
- elapsed = std::chrono::system_clock::now() - tp;
- std::cout << "exhaustive min:" << min.second << " distance:" << min.first << " elapsed:" << elapsed.count() / (1000 * 1000) << std::endl;
- tp = std::chrono::system_clock::now();
- #endif
- for (int i = 0; i < N; ++i) { sorted[i] = points[i]; }
- std::sort(std::begin(sorted), std::end(sorted), [](const Point<CT> &l, const Point<CT> &r) { return l.x < r.x; });
- min = recursiveSearch(sorted.data(), 0, N);
- elapsed = std::chrono::system_clock::now() - tp;
- std::cout << "recursive min:" << min.second << " distance:" << min.first << " elapsed:" << elapsed.count() / (1000 * 1000) << std::endl;
- tp = std::chrono::system_clock::now();
- for (int i = 0; i < N; ++i) { sorted[i] = points[i]; }
- std::sort(std::begin(sorted), std::end(sorted), [](const Point<CT> &l, const Point<CT> &r) { return l.x < r.x; });
- min = neighbourSearch(sorted.data(), 0, N);
- elapsed = std::chrono::system_clock::now() - tp;
- std::cout << "neighbour min:" << min.second << " distance:" << min.first << " elapsed:" << elapsed.count() / (1000 * 1000) << std::endl;
- }
- int main(int,char**)
- {
- timeOperations<float>();
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement