Experiment on pulling NVTX3 via CPM

cpp/benchmarks/CMakeLists.txt

@@ -83,7 +83,8 @@ ConfigureBench(HAUSDORFF_BENCH
     hausdorff_benchmark.cpp)
 
 ConfigureNVBench(DISTANCES_BENCH
-    pairwise_linestring_distance.cu)
+    pairwise_linestring_distance.cu
+    pairwise_point_polygon_distance.cu)
 
 ConfigureNVBench(POINTS_IN_RANGE_BENCH
     points_in_range.cu)

cpp/benchmarks/pairwise_point_polygon_distance.cu

@@ -0,0 +1,134 @@
+/*
+ * Copyright (c) 2023, 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 <nvbench/nvbench.cuh>
+
+#include <cuspatial_test/geometry_generator.cuh>
+
+#include <cuspatial/experimental/point_polygon_distance.cuh>
+#include <cuspatial/vec_2d.hpp>
+
+#include <rmm/cuda_stream_view.hpp>
+#include <rmm/exec_policy.hpp>
+
+using namespace cuspatial;
+using namespace cuspatial::test;
+
+template <typename T>
+void pairwise_point_polygon_distance_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;
+  rmm::cuda_stream_view stream{rmm::cuda_stream_default};
+
+  auto const num_pairs{static_cast<std::size_t>(state.get_int64("num_pairs"))};
+
+  auto const num_polygons_per_multipolygon{
+    static_cast<std::size_t>(state.get_int64("num_polygons_per_multipolygon"))};
+  auto const num_holes_per_polygon{
+    static_cast<std::size_t>(state.get_int64("num_holes_per_polygon"))};
+  auto const num_edges_per_ring{static_cast<std::size_t>(state.get_int64("num_edges_per_ring"))};
+
+  auto const num_points_per_multipoint{
+    static_cast<std::size_t>(state.get_int64("num_points_per_multipoint"))};
+
+  auto mpoly_generator_param = multipolygon_generator_parameter<T>{
+    num_pairs, num_polygons_per_multipolygon, num_holes_per_polygon, num_edges_per_ring};
+
+  auto mpoint_generator_param = multipoint_generator_parameter<T>{
+    num_pairs, num_points_per_multipoint, vec_2d<T>{-1, -1}, vec_2d<T>{0, 0}};
+
+  auto multipolygons = generate_multipolygon_array<T>(mpoly_generator_param, stream);
+  auto multipoints   = generate_multipoint_array<T>(mpoint_generator_param, stream);
+
+  auto distances = rmm::device_vector<T>(num_pairs);
+  auto out_it    = distances.begin();
+
+  auto mpoly_view  = multipolygons.range();
+  auto mpoint_view = multipoints.range();
+
+  state.add_element_count(num_pairs, "NumPairs");
+  state.add_element_count(mpoly_generator_param.num_polygons(), "NumPolygons");
+  state.add_element_count(mpoly_generator_param.num_rings(), "NumRings");
+  state.add_element_count(mpoly_generator_param.num_coords(), "NumPoints (in mpoly)");
+  state.add_element_count(static_cast<std::size_t>(mpoly_generator_param.num_coords() *
+                                                   mpoly_generator_param.num_rings() *
+                                                   mpoly_generator_param.num_polygons()),
+                          "Multipolygon Complexity");
+  state.add_element_count(mpoint_generator_param.num_points(), "NumPoints (in multipoints)");
+
+  state.add_global_memory_reads<T>(
+    mpoly_generator_param.num_coords() + mpoint_generator_param.num_points(),
+    "CoordinatesReadSize");
+  state.add_global_memory_reads<std::size_t>(
+    (mpoly_generator_param.num_rings() + 1) + (mpoly_generator_param.num_polygons() + 1) +
+      (mpoly_generator_param.num_multipolygons + 1) + (mpoint_generator_param.num_multipoints + 1),
+    "OffsetsDataSize");
+
+  state.add_global_memory_writes<T>(num_pairs);
+
+  state.exec(nvbench::exec_tag::sync,
+             [&mpoly_view, &mpoint_view, &out_it, &stream](nvbench::launch& launch) {
+               pairwise_point_polygon_distance(mpoint_view, mpoly_view, out_it, stream);
+             });
+}
+
+using floating_point_types = nvbench::type_list<float, double>;
+
+// Benchmark scalability with simple multipolygon (3 sides, 0 hole, 1 poly)
+NVBENCH_BENCH_TYPES(pairwise_point_polygon_distance_benchmark,
+                    NVBENCH_TYPE_AXES(floating_point_types))
+  .set_type_axes_names({"CoordsType"})
+  .add_int64_axis("num_pairs", {1, 1'00, 10'000, 1'000'000, 100'000'000})
+  .add_int64_axis("num_polygons_per_multipolygon", {1})
+  .add_int64_axis("num_holes_per_polygon", {0})
+  .add_int64_axis("num_edges_per_ring", {3})
+  .add_int64_axis("num_points_per_multipoint", {1})
+  .set_name("point_polygon_distance_benchmark_simple_polygon");
+
+// Benchmark scalability with complex multipolygon (100 sides, 10 holes, 3 polys)
+NVBENCH_BENCH_TYPES(pairwise_point_polygon_distance_benchmark,
+                    NVBENCH_TYPE_AXES(floating_point_types))
+  .set_type_axes_names({"CoordsType"})
+  .add_int64_axis("num_pairs", {1'000, 10'000, 100'000, 1'000'000})
+  .add_int64_axis("num_polygons_per_multipolygon", {2})
+  .add_int64_axis("num_holes_per_polygon", {3})
+  .add_int64_axis("num_edges_per_ring", {50})
+  .add_int64_axis("num_points_per_multipoint", {1})
+  .set_name("point_polygon_distance_benchmark_complex_polygon");
+
+// // Benchmark impact of rings (100K pairs, 1 polygon, 3 sides)
+NVBENCH_BENCH_TYPES(pairwise_point_polygon_distance_benchmark,
+                    NVBENCH_TYPE_AXES(floating_point_types))
+  .set_type_axes_names({"CoordsType"})
+  .add_int64_axis("num_pairs", {10'000})
+  .add_int64_axis("num_polygons_per_multipolygon", {1})
+  .add_int64_axis("num_holes_per_polygon", {0, 10, 100, 1000})
+  .add_int64_axis("num_edges_per_ring", {3})
+  .add_int64_axis("num_points_per_multipoint", {1})
+  .set_name("point_polygon_distance_benchmark_ring_numbers");
+
+// Benchmark impact of rings (1M pairs, 1 polygon, 0 holes, 3 sides)
+NVBENCH_BENCH_TYPES(pairwise_point_polygon_distance_benchmark,
+                    NVBENCH_TYPE_AXES(floating_point_types))
+  .set_type_axes_names({"CoordsType"})
+  .add_int64_axis("num_pairs", {100})
+  .add_int64_axis("num_polygons_per_multipolygon", {1})
+  .add_int64_axis("num_holes_per_polygon", {0})
+  .add_int64_axis("num_edges_per_ring", {3})
+  .add_int64_axis("num_points_per_multipoint", {50, 5'00, 5'000, 50'000, 500'000})
+  .set_name("point_polygon_distance_benchmark_points_in_multipoint");

cpp/include/cuspatial/detail/nvtx/ranges.hpp

@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2020, NVIDIA CORPORATION.
+ * Copyright (c) 2020-2023, NVIDIA CORPORATION.
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
  * you may not use this file except in compliance with the License.
@@ -20,16 +20,16 @@
 
 namespace cuspatial {
 /**
- * @brief Tag type for libcudf's NVTX domain.
+ * @brief Tag type for libcuspatial's NVTX domain.
  */
 struct libcuspatial_domain {
-  static constexpr char const* name{"libcuspatial"};  ///< Name of the libcudf domain
+  static constexpr char const* name{"libcuspatial"};  ///< Name of the libcuspatial domain
 };
 
 /**
- * @brief Alias for an NVTX range in the libcudf domain.
+ * @brief Alias for an NVTX range in the libcuspatial domain.
  */
-using thread_range = ::nvtx3::domain_thread_range<libcudf_domain>;
+using thread_range = ::nvtx3::domain_thread_range<libcuspatial_domain>;
 
 }  // namespace cuspatial
 

cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh

@@ -18,6 +18,7 @@
 
 #include <cuspatial/cuda_utils.hpp>
 #include <cuspatial/detail/iterator.hpp>
+#include <cuspatial/detail/nvtx/ranges.hpp>
 #include <cuspatial/detail/utility/device_atomics.cuh>
 #include <cuspatial/detail/utility/linestring.cuh>
 #include <cuspatial/detail/utility/zero_data.cuh>
@@ -120,6 +121,8 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints,
                                          OutputIt distances_first,
                                          rmm::cuda_stream_view stream)
 {
+  CUSPATIAL_FUNC_RANGE();
+
   using T       = typename MultiPointRange::element_t;
   using index_t = typename MultiPointRange::index_t;
 

cpp/include/cuspatial_test/geometry_generator.cuh

@@ -14,6 +14,7 @@
  * limitations under the License.
  */
 
+#include <cuspatial_test/random.cuh>
 #include <cuspatial_test/vector_factories.cuh>
 
 #include <cuspatial/cuda_utils.hpp>
@@ -249,5 +250,60 @@ auto generate_multipolygon_array(multipolygon_generator_parameter<T> params,
                                                          std::move(coordinates));
 }
 
+/**
+ * @brief Struct to store the parameters of the multipoint aray
+ *
+ * @tparam T Type of the coordinates
+ */
+template <typename T>
+struct multipoint_generator_parameter {
+  using element_t = T;
+
+  std::size_t num_multipoints;
+  std::size_t num_points_per_multipoints;
+  vec_2d<T> lower_left;
+  vec_2d<T> upper_right;
+
+  CUSPATIAL_HOST_DEVICE std::size_t num_points()
+  {
+    return num_multipoints * num_points_per_multipoints;
+  }
+};
+
+/**
+ * @brief Helper to generate random multipoints within a range
+ *
+ * @tparam T The floating point type for the coordinates
+ * @param params Parameters to specify for the multipoints
+ * @param stream The CUDA stream to use for device memory operations and kernel launches
+ * @return a cuspatial::test::multipoint_array object
+ */
+template <typename T>
+auto generate_multipoint_array(multipoint_generator_parameter<T> params,
+                               rmm::cuda_stream_view stream)
+{
+  rmm::device_uvector<vec_2d<T>> coordinates(params.num_points(), stream);
+  rmm::device_uvector<std::size_t> offsets(params.num_multipoints + 1, stream);
+
+  thrust::sequence(rmm::exec_policy(stream),
+                   offsets.begin(),
+                   offsets.end(),
+                   std::size_t{0},
+                   params.num_points_per_multipoints);
+
+  auto engine_x = deterministic_engine(params.num_points());
+  auto engine_y = deterministic_engine(2 * params.num_points());
+
+  auto x_dist = make_uniform_dist(params.lower_left.x, params.upper_right.x);
+  auto y_dist = make_uniform_dist(params.lower_left.y, params.upper_right.y);
+
+  auto point_gen =
+    point_generator(params.lower_left, params.upper_right, engine_x, engine_y, x_dist, y_dist);
+
+  thrust::tabulate(rmm::exec_policy(stream), coordinates.begin(), coordinates.end(), point_gen);
+
+  return make_multipoints_array(std::move(offsets), std::move(coordinates));
+}
+
 }  // namespace test
 }  // namespace cuspatial

cpp/include/cuspatial_test/random.cuh

@@ -33,6 +33,8 @@
 #include <algorithm>
 #include <memory>
 
+using namespace cuspatial;
+
 /**
  * @brief Identifies a probability distribution type.
  */
@@ -150,14 +152,21 @@ struct value_generator {
 template <typename T, typename Generator>
 struct point_generator {
   using Cart2D = cuspatial::vec_2d<T>;
-  value_generator<T, Generator> vgen;
-
-  point_generator(T lower_bound, T upper_bound, thrust::minstd_rand& engine, Generator gen)
-    : vgen(lower_bound, upper_bound, engine, gen)
+  value_generator<T, Generator> vgenx;
+  value_generator<T, Generator> vgeny;
+
+  point_generator(vec_2d<T> lower_left,
+                  vec_2d<T> upper_right,
+                  thrust::minstd_rand& engine_x,
+                  thrust::minstd_rand& engine_y,
+                  Generator gen_x,
+                  Generator gen_y)
+    : vgenx(lower_left.x, upper_right.x, engine_x, gen_x),
+      vgeny(lower_left.y, upper_right.y, engine_y, gen_y)
   {
   }
 
-  __device__ Cart2D operator()(size_t n) { return {vgen(n), vgen(n)}; }
+  __device__ Cart2D operator()(size_t n) { return {vgenx(n), vgeny(n)}; }
 };
 
 /**

cpp/include/cuspatial_test/vector_factories.cuh

@@ -270,11 +270,22 @@ auto make_multilinestring_array(std::initializer_list<std::size_t> geometry_inl,
 template <typename GeometryArray, typename CoordinateArray>
 class multipoint_array {
  public:
-  multipoint_array(GeometryArray geometry_offsets_array, CoordinateArray coordinate_array)
+  using geometry_t = typename GeometryArray::value_type;
+  using coord_t    = typename CoordinateArray::value_type;
+
+  multipoint_array(thrust::device_vector<geometry_t> geometry_offsets_array,
+                   thrust::device_vector<coord_t> coordinate_array)
     : _geometry_offsets(geometry_offsets_array), _coordinates(coordinate_array)
   {
   }
 
+  multipoint_array(rmm::device_uvector<geometry_t>&& geometry_offsets_array,
+                   rmm::device_uvector<coord_t>&& coordinate_array)
+    : _geometry_offsets(std::move(geometry_offsets_array)),
+      _coordinates(std::move(coordinate_array))
+  {
+  }
+
   /// Return the number of multipoints
   auto size() { return _geometry_offsets.size() - 1; }
 
@@ -322,8 +333,20 @@ auto make_multipoints_array(std::initializer_list<std::initializer_list<vec_2d<T
       return init;
     });
 
-  return multipoint_array{rmm::device_vector<std::size_t>(offsets),
-                          rmm::device_vector<vec_2d<T>>(coordinates)};
+  return multipoint_array<rmm::device_vector<std::size_t>, rmm::device_vector<vec_2d<T>>>{
+    rmm::device_vector<std::size_t>(offsets), rmm::device_vector<vec_2d<T>>(coordinates)};
+}
+
+/**
+ * @brief Factory method to construct multipoint array by moving the offsets and coordinates from
+ * `rmm::device_uvector`.
+ */
+template <typename IndexType, typename T>
+auto make_multipoints_array(rmm::device_uvector<IndexType> geometry_offsets,
+                            rmm::device_uvector<vec_2d<T>> coords)
+{
+  return multipoint_array<rmm::device_uvector<std::size_t>, rmm::device_uvector<vec_2d<T>>>{
+    std::move(geometry_offsets), std::move(coords)};
 }
 
 }  // namespace test

cpp/tests/experimental/spatial/point_distance_test.cu

@@ -57,9 +57,12 @@ struct PairwisePointDistanceTest : public ::testing::Test {
     std::size_t seed,
     rmm::cuda_stream_view stream = rmm::cuda_stream_default)
   {
-    auto engine  = deterministic_engine(0);
-    auto uniform = make_normal_dist<T>(0.0, 1.0);
-    auto pgen    = point_generator(T{0.0}, T{1.0}, engine, uniform);
+    auto engine_x  = deterministic_engine(0);
+    auto engine_y  = deterministic_engine(0);
+    auto uniform_x = make_normal_dist<T>(0.0, 1.0);
+    auto uniform_y = make_normal_dist<T>(0.0, 1.0);
+    auto pgen =
+      point_generator(vec_2d<T>{0, 0}, vec_2d<T>{1, 1}, engine_x, engine_y, uniform_x, uniform_y);
     rmm::device_vector<vec_2d<T>> points(num_points);
     auto counting_iter = thrust::make_counting_iterator(seed);
     thrust::transform(