Add benchmark for points_in_spatial_window

cpp/benchmarks/CMakeLists.txt

@@ -84,3 +84,6 @@ ConfigureBench(HAUSDORFF_BENCH
 
 ConfigureNVBench(DISTANCES_BENCH
     pairwise_linestring_distance.cu)
+
+ConfigureNVBench(SPATIAL_WINDOW_BENCH
+    spatial_window.cu)

cpp/benchmarks/pairwise_linestring_distance.cu

@@ -20,6 +20,7 @@
 #include <cuspatial/detail/iterator.hpp>
 #include <cuspatial/experimental/linestring_distance.cuh>
 #include <cuspatial/experimental/type_utils.hpp>
+#include <cuspatial/vec_2d.hpp>
 
 #include <rmm/device_vector.hpp>
 #include <rmm/exec_policy.hpp>
@@ -93,8 +94,8 @@ void pairwise_linestring_distance_benchmark(nvbench::state& state, nvbench::type
   // TODO: to be replaced by nvbench fixture once it's ready
   cuspatial::rmm_pool_raii rmm_pool;
 
-  auto const num_string_pairs{state.get_int64("NumStrings")},
-    num_segments_per_string{state.get_int64("NumSegmentsPerString")};
+  auto const num_string_pairs{state.get_int64("NumStrings")};
+  auto const num_segments_per_string{state.get_int64("NumSegmentsPerString")};
 
   auto [ls1, ls1_offset] =
     generate_linestring<T>(num_string_pairs, num_segments_per_string, 1, {0, 0});

cpp/benchmarks/spatial_window.cu

@@ -0,0 +1,114 @@
+/*
+ * Copyright (c) 2022, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include <benchmarks/fixture/rmm_pool_raii.hpp>
+#include <benchmarks/utility/random.cuh>
+
+#include <cuspatial/detail/iterator.hpp>
+#include <cuspatial/spatial_window.hpp>
+#include <cuspatial/vec_2d.hpp>
+
+#include <rmm/device_uvector.hpp>
+#include <rmm/exec_policy.hpp>
+
+#include <nvbench/nvbench.cuh>
+
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/random/linear_congruential_engine.h>
+#include <thrust/random/normal_distribution.h>
+#include <thrust/random/uniform_int_distribution.h>
+
+#include <memory>
+
+using namespace cuspatial;
+
+/**
+ * @brief Helper to generate random points within a rectangular window
+ *
+ * @p begin and @p end must be iterators to device-accessible memory
+ *
+ * @tparam PointsIter The type of the iterator to the output points container
+ * @tparam T The floating point type for the coordinates
+ * @param begin The start of the range of points to generate
+ * @param end The end of the range of points to generate
+ *
+ * @param window_min the lower left window corner
+ * @param window_max the upper right window corner
+ *
+ */
+template <class PointsIter, typename T>
+void generate_points(PointsIter begin, PointsIter end, vec_2d<T> window_min, vec_2d<T> window_max)
+{
+  auto engine_x = deterministic_engine(std::distance(begin, end));
+  auto engine_y = deterministic_engine(2 * std::distance(begin, end));
+
+  auto x_dist = make_uniform_dist(window_min.x, window_max.x);
+  auto y_dist = make_uniform_dist(window_min.y, window_max.y);
+
+  auto x_gen = value_generator{window_min.x, window_max.x, engine_x, x_dist};
+  auto y_gen = value_generator{window_min.y, window_max.y, engine_y, y_dist};
+
+  thrust::tabulate(rmm::exec_policy(), begin, end, [x_gen, y_gen] __device__(size_t n) mutable {
+    return vec_2d<T>{x_gen(n), y_gen(n)};
+  });
+}
+
+template <typename T>
+void points_in_spatial_window_benchmark(nvbench::state& state, nvbench::type_list<T>)
+{
+  // TODO: to be replaced by nvbench fixture once it's ready
+  cuspatial::rmm_pool_raii rmm_pool;
+
+  auto const num_points{state.get_int64("NumPoints")};
+
+  auto window_min = vec_2d<T>{-100, -100};
+  auto window_max = vec_2d<T>{100, 100};
+
+  auto range_min = vec_2d<T>{-200, -200};
+  auto range_max = vec_2d<T>{200, 200};
+
+  auto d_points = rmm::device_uvector<vec_2d<T>>(num_points, rmm::cuda_stream_default);
+  generate_points(d_points.begin(), d_points.end(), range_min, range_max);
+
+  auto d_x = rmm::device_uvector<T>(num_points, rmm::cuda_stream_default);
+  auto d_y = rmm::device_uvector<T>(num_points, rmm::cuda_stream_default);
+
+  thrust::transform(
+    rmm::exec_policy(), d_points.begin(), d_points.end(), d_x.begin(), [] __device__(auto point) {
+      return point.x;
+    });
+  thrust::transform(
+    rmm::exec_policy(), d_points.begin(), d_points.end(), d_y.begin(), [] __device__(auto point) {
+      return point.y;
+    });
+
+  auto xs = cudf::column(cudf::data_type{cudf::type_to_id<T>()}, num_points, d_x.release());
+  auto ys = cudf::column(cudf::data_type{cudf::type_to_id<T>()}, num_points, d_y.release());
+
+  state.add_element_count(num_points);
+
+  points_in_spatial_window(window_min.x, window_max.x, window_min.y, window_max.y, xs, ys);
+
+  state.exec(nvbench::exec_tag::sync, [&](nvbench::launch& launch) {
+    auto points_in =
+      points_in_spatial_window(window_min.x, window_max.x, window_min.y, window_max.y, xs, ys);
+  });
+}
+
+using floating_point_types = nvbench::type_list<float, double>;
+NVBENCH_BENCH_TYPES(points_in_spatial_window_benchmark, NVBENCH_TYPE_AXES(floating_point_types))
+  .set_type_axes_names({"CoordsType"})
+  .add_int64_axis("NumPoints", {100'000, 1'000'000, 10'000'000, 100'000'000});

cpp/benchmarks/utility/random.cuh

@@ -0,0 +1,152 @@
+/*
+ * Copyright (c) 2020-2022, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+
+#include <cuspatial/cuda_utils.hpp>
+#include <cuspatial/error.hpp>
+
+#include <rmm/device_uvector.hpp>
+
+#include <thrust/execution_policy.h>
+#include <thrust/random.h>
+#include <thrust/random/normal_distribution.h>
+#include <thrust/random/uniform_int_distribution.h>
+#include <thrust/tabulate.h>
+
+#include <cuda/std/type_traits>
+
+#include <algorithm>
+#include <memory>
+
+/**
+ * @brief Identifies a probability distribution type.
+ */
+enum class distribution_id : int8_t {
+  UNIFORM,    ///< Uniform sampling between the given bounds. Provides the best coverage of the
+              ///< overall value range. Real data rarely has this distribution.
+  NORMAL,     ///< Gaussian sampling - most samples are close to the middle of the range. Good for
+              ///< simulating real-world numeric data.
+  GEOMETRIC,  ///< Geometric sampling - highest chance to sample close to the lower bound. Good for
+              ///< simulating real data with asymmetric distribution (unsigned values, timestamps).
+};
+
+/**
+ * @brief Real Type that has atleast number of bits of integral type in its mantissa.
+ *  number of bits of integrals < 23 bits of mantissa in float
+ * to allow full range of integer bits to be generated.
+ * @tparam T integral type
+ */
+template <typename T>
+using integral_to_realType =
+  std::conditional_t<std::is_floating_point_v<T>,
+                     T,
+                     std::conditional_t<sizeof(T) * 8 <= 23, float, double>>;
+
+/**
+ * @brief Generates a normal distribution between zero and upper_bound.
+ */
+template <typename T>
+auto make_normal_dist(T lower_bound, T upper_bound)
+{
+  using realT    = integral_to_realType<T>;
+  T const mean   = lower_bound + (upper_bound - lower_bound) / 2;
+  T const stddev = (upper_bound - lower_bound) / 6;
+  return thrust::random::normal_distribution<realT>(mean, stddev);
+}
+
+template <typename T, std::enable_if_t<std::is_integral_v<T>, T>* = nullptr>
+auto make_uniform_dist(T range_start, T range_end)
+{
+  return thrust::uniform_int_distribution<T>(range_start, range_end);
+}
+
+template <typename T, std::enable_if_t<std::is_floating_point_v<T>>* = nullptr>
+auto make_uniform_dist(T range_start, T range_end)
+{
+  return thrust::uniform_real_distribution<T>(range_start, range_end);
+}
+
+template <typename T>
+double geometric_dist_p(T range_size)
+{
+  constexpr double percentage_in_range = 0.99;
+  double const p                       = 1 - exp(log(1 - percentage_in_range) / range_size);
+  return p ? p : std::numeric_limits<double>::epsilon();
+}
+
+/**
+ * @brief Generates a geometric distribution between lower_bound and upper_bound.
+ * This distribution is an approximation generated using normal distribution.
+ *
+ * @tparam T Result type of the number to produce.
+ */
+template <typename T>
+class geometric_distribution : public thrust::random::normal_distribution<integral_to_realType<T>> {
+  using realType = integral_to_realType<T>;
+  using super_t  = thrust::random::normal_distribution<realType>;
+  T _lower_bound;
+  T _upper_bound;
+
+ public:
+  using result_type = T;
+  __host__ __device__ explicit geometric_distribution(T lower_bound, T upper_bound)
+    : super_t(0, std::labs(upper_bound - lower_bound) / 4.0),
+      _lower_bound(lower_bound),
+      _upper_bound(upper_bound)
+  {
+  }
+
+  template <typename UniformRandomNumberGenerator>
+  __host__ __device__ result_type operator()(UniformRandomNumberGenerator& urng)
+  {
+    return _lower_bound < _upper_bound ? std::abs(super_t::operator()(urng)) + _lower_bound
+                                       : _lower_bound - std::abs(super_t::operator()(urng));
+  }
+};
+
+template <typename T, typename Generator>
+struct value_generator {
+  using result_type = T;
+
+  value_generator(T lower_bound, T upper_bound, thrust::minstd_rand& engine, Generator gen)
+    : lower_bound(std::min(lower_bound, upper_bound)),
+      upper_bound(std::max(lower_bound, upper_bound)),
+      engine(engine),
+      dist(gen)
+  {
+  }
+
+  __device__ T operator()(size_t n)
+  {
+    engine.discard(n);
+    if constexpr (std::is_integral_v<T> && std::is_floating_point_v<decltype(dist(engine))>) {
+      return std::clamp(static_cast<T>(std::round(dist(engine))), lower_bound, upper_bound);
+    } else {
+      return std::clamp(dist(engine), lower_bound, upper_bound);
+    }
+  }
+
+  T lower_bound;
+  T upper_bound;
+  thrust::minstd_rand engine;
+  Generator dist;
+};
+
+/**
+ * @brief LCG pseudo-random engine.
+ */
+auto deterministic_engine(unsigned seed) { return thrust::minstd_rand{seed}; }

cpp/include/cuspatial/spatial_window.hpp

@@ -16,6 +16,8 @@
 
 #pragma once
 
+#include <cudf/column/column_view.hpp>
+#include <cudf/table/table.hpp>
 #include <cudf/types.hpp>
 
 #include <rmm/mr/device/per_device_resource.hpp>