Benchmarking the speed of side predicates

Readme.txt

The following file benchmarks the side predicates found in the current (as of writing this) develop branch of Boost.Geometry.

At best this provides a rough estimate because it is based on artifical data and run in isolation.

I used the following command to compile the program:

g++ -I../geometry/include -O3 -march=native benchmark.cpp

benchmark.cpp

#include <iostream>
#include <random>
#include <chrono>

#include <boost/geometry.hpp>

#include <boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp>

const int tests = 10'000'000;

namespace bg = boost::geometry;
using p = bg::model::d2::point_xy<double>;
using robust = bg::strategy::side::side_robust<double>;
using adaptive0 = bg::strategy::side::side_robust<double, 0>;
using adaptive1 = bg::strategy::side::side_robust<double, 1>;
using adaptive2 = bg::strategy::side::side_robust<double, 2>;
using approx = bg::strategy::side::side_by_triangle<double>;

void continuous_random(std::vector<p>& problems, std::mt19937& gen) {
  for (int i = 0; i < tests; ++i) {
    std::uniform_real_distribution<> dis(-40.0, 40.0);
    problems.emplace_back(dis(gen), dis(gen));
    problems.emplace_back(dis(gen), dis(gen));
    problems.emplace_back(dis(gen), dis(gen));
  }
}

void continuous_random_exponent(std::vector<p>& problems, std::mt19937& gen) {
  for (int i = 0; i < tests; ++i) {
    std::uniform_real_distribution<> dis(-40.0, 40.0);
    problems.emplace_back(dis(gen) * std::pow(10.0, dis(gen)),
         dis(gen) * std::pow(10.0, dis(gen)));
    problems.emplace_back(dis(gen) * std::pow(10.0, dis(gen)),
         dis(gen) * std::pow(10.0, dis(gen)));
    problems.emplace_back(dis(gen) * std::pow(10.0, dis(gen)),
         dis(gen) * std::pow(10.0, dis(gen)));
  }
}

void perturbed_grid(std::vector<p>& problems, std::mt19937& gen) {
  for (int i = 0; i < tests; ++i) {
    std::uniform_int_distribution<> dis(-40, 40);
    p p1(dis(gen), dis(gen));
    p p2(dis(gen), dis(gen));
    p p3(dis(gen), dis(gen));

    namespace trans = bg::strategy::transform;
    if (dis(gen) > 0) {
      std::uniform_real_distribution<> disa(0.0, 360.0);
      trans::rotate_transformer<bg::degree, double, 2, 2> rotate(disa(gen));
      p rp1, rp2, rp3;
      bg::transform(p1, rp1, rotate);
      bg::transform(p2, rp2, rotate);
      bg::transform(p3, rp3, rotate);
      p1 = rp1;
      p2 = rp2;
      p3 = rp3;
    }
    if (dis(gen) > 0) {
      std::uniform_real_distribution<> diss(-40.0, 40.0);
      trans::scale_transformer<double, 2, 2> scale(diss(gen));
      p sp1, sp2, sp3;
      bg::transform(p1, sp1, scale);
      bg::transform(p2, sp2, scale);
      bg::transform(p3, sp3, scale);
      p1 = sp1;
      p2 = sp2;
      p3 = sp3;
    }
    problems.push_back(p1);
    problems.push_back(p2);
    problems.push_back(p3);
  }
}

void purposefully_degenerate(std::vector<p>& problems) {
  for (int i = 0; i < tests; ++i) {
    auto v = i % 3;
    problems.emplace_back(1, 1);
    problems.emplace_back(1e20, v == 0 ? 1e20 : ( v == 1 ? std::nextafter(1e20, 1e21) : std::nextafter(1e20, 1e19) ) );
    problems.emplace_back(1e40, 1e40);
  }
}

template<typename Strategy>
void compute_and_print(std::vector<p> const& problems) {
  int sum = 0;
  auto start = std::chrono::system_clock::now();
  for(int i = 0; i < tests; ++i)
    sum += Strategy::apply(problems[ 3 * i ], problems[ 3 * i + 1], problems[ 3 * i + 2]);
  auto end = std::chrono::system_clock::now();
  auto elapsed =
    std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
  std::cout << "duration: " << elapsed.count() << "ms, sum: " << sum << "\n";
}

void compute_and_print_all(std::vector<p> const& problems) {
  std::cout << "Approx:\n";
  compute_and_print<approx>(problems);
  std::cout << "Robust:\n";
  compute_and_print<robust>(problems);
  std::cout << "Adaptive 0 (No extra precision):\n";
  compute_and_print<adaptive0>(problems);
  std::cout << "Adaptive 1 (should be correct, if det is around eps):\n";
  compute_and_print<adaptive1>(problems);
  std::cout << "Adaptive 2 (should be correct, if det is around eps^2):\n";
  compute_and_print<adaptive2>(problems);
}

int main() {
  namespace bg = boost::geometry;
  std::random_device rd;
  std::mt19937 gen(rd());
  std::vector<p> problems;
  problems.reserve(3 * tests);

  std::cout << "Continuous random: \n";
  continuous_random(problems, gen);
  compute_and_print_all(problems);
  problems.clear();
  std::cout << "\nRandom exponent: \n";
  continuous_random_exponent(problems, gen);
  compute_and_print_all(problems);
  problems.clear();
  std::cout << "\nPerturbed grid: \n";
  perturbed_grid(problems, gen);
  compute_and_print_all(problems);
  problems.clear();
  purposefully_degenerate(problems);
  std::cout << "\nPurposefully degenerate values: \n";
  compute_and_print_all(problems);

  return 0;
}