[skip-ci] Experimental work with Binary Operation Dispatcher (Python)

cpp/CMakeLists.txt

@@ -134,6 +134,8 @@ add_library(cuspatial
     src/spatial/point_distance.cu
     src/spatial/point_linestring_distance.cu
     src/spatial/point_polygon_distance.cu
+    src/spatial/linestring_polygon_distance.cu
+    src/spatial/polygon_distance.cu
     src/spatial/point_linestring_nearest_points.cu
     src/spatial/sinusoidal_projection.cu
     src/trajectory/derive_trajectories.cu

cpp/include/cuspatial/assert.cuh

@@ -0,0 +1,35 @@
+/*
+ * 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.
+ */
+
+#pragma once
+
+#include <cuda_runtime.h>
+
+/**
+ * @brief `assert`-like macro for device code
+ *
+ * This is effectively the same as the standard `assert` macro, except it
+ * relies on the `__PRETTY_FUNCTION__` macro which is specific to GCC and Clang
+ * to produce better assert messages.
+ */
+#if !defined(NDEBUG) && defined(__CUDA_ARCH__) && (defined(__clang__) || defined(__GNUC__))
+#define __ASSERT_STR_HELPER(x) #x
+#define cuspatial_assert(e)   \
+  ((e) ? static_cast<void>(0) \
+       : __assert_fail(__ASSERT_STR_HELPER(e), __FILE__, __LINE__, __PRETTY_FUNCTION__))
+#else
+#define cuspatial_assert(e) (static_cast<void>(0))
+#endif

cpp/include/cuspatial/detail/utility/validation.hpp

@@ -35,10 +35,10 @@
  * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout
  */
 #define CUSPATIAL_EXPECTS_VALID_LINESTRING_SIZES(num_linestring_points, num_linestring_offsets) \
-  CUSPATIAL_EXPECTS(num_linestring_offsets > 0,                                                 \
-                    "Polygon offsets must contain at least one (1) value");                     \
-  CUSPATIAL_EXPECTS(num_linestring_points >= 2 * (num_linestring_offsets - 1),                  \
-                    "Each linestring must have at least two vertices");
+  CUSPATIAL_HOST_DEVICE_EXPECTS(num_linestring_offsets > 0,                                     \
+                                "Polygon offsets must contain at least one (1) value");         \
+  CUSPATIAL_HOST_DEVICE_EXPECTS(num_linestring_points >= 2 * (num_linestring_offsets - 1),      \
+                                "Each linestring must have at least two vertices");
 
 /**
  * @brief Macro for validating the data array sizes for multilinestrings.
@@ -57,10 +57,10 @@
  * [1]: https://github.com/geoarrow/geoarrow/blob/main/format.md
  * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout
  */
-#define CUSPATIAL_EXPECTS_VALID_MULTILINESTRING_SIZES(                              \
-  num_linestring_points, num_multilinestring_offsets, num_linestring_offsets)       \
-  CUSPATIAL_EXPECTS(num_multilinestring_offsets > 0,                                \
-                    "Multilinestring offsets must contain at least one (1) value"); \
+#define CUSPATIAL_EXPECTS_VALID_MULTILINESTRING_SIZES(                                          \
+  num_linestring_points, num_multilinestring_offsets, num_linestring_offsets)                   \
+  CUSPATIAL_HOST_DEVICE_EXPECTS(num_multilinestring_offsets > 0,                                \
+                                "Multilinestring offsets must contain at least one (1) value"); \
   CUSPATIAL_EXPECTS_VALID_LINESTRING_SIZES(num_linestring_points, num_linestring_offsets);
 
 /**
@@ -84,15 +84,16 @@
  * [1]: https://github.com/geoarrow/geoarrow/blob/main/format.md
  * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout
  */
-#define CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES(                                                    \
-  num_poly_points, num_poly_offsets, num_poly_ring_offsets)                                       \
-  CUSPATIAL_EXPECTS(num_poly_offsets > 0, "Polygon offsets must contain at least one (1) value"); \
-  CUSPATIAL_EXPECTS(num_poly_ring_offsets > 0,                                                    \
-                    "Polygon ring offsets must contain at least one (1) value");                  \
-  CUSPATIAL_EXPECTS(num_poly_ring_offsets >= num_poly_offsets,                                    \
-                    "Each polygon must have at least one (1) ring");                              \
-  CUSPATIAL_EXPECTS(num_poly_points >= 4 * (num_poly_ring_offsets - 1),                           \
-                    "Each ring must have at least four (4) vertices");
+#define CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES(                                               \
+  num_poly_points, num_poly_offsets, num_poly_ring_offsets)                                  \
+  CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_offsets > 0,                                        \
+                                "Polygon offsets must contain at least one (1) value");      \
+  CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_ring_offsets > 0,                                   \
+                                "Polygon ring offsets must contain at least one (1) value"); \
+  CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_ring_offsets >= num_poly_offsets,                   \
+                                "Each polygon must have at least one (1) ring");             \
+  CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_points >= 4 * (num_poly_ring_offsets - 1),          \
+                                "Each ring must have at least four (4) vertices");
 
 /**
  * @brief Macro for validating the data array sizes for a multipolygon.
@@ -116,8 +117,8 @@
  * [1]: https://github.com/geoarrow/geoarrow/blob/main/format.md
  * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout
  */
-#define CUSPATIAL_EXPECTS_VALID_MULTIPOLYGON_SIZES(                                \
-  num_poly_points, num_multipoly_offsets, num_poly_offsets, num_poly_ring_offsets) \
-  CUSPATIAL_EXPECTS(num_multipoly_offsets > 0,                                     \
-                    "Multipolygon offsets must contain at least one (1) value");   \
+#define CUSPATIAL_EXPECTS_VALID_MULTIPOLYGON_SIZES(                                          \
+  num_poly_points, num_multipoly_offsets, num_poly_offsets, num_poly_ring_offsets)           \
+  CUSPATIAL_HOST_DEVICE_EXPECTS(num_multipoly_offsets > 0,                                   \
+                                "Multipolygon offsets must contain at least one (1) value"); \
   CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES(num_poly_points, num_poly_offsets, num_poly_ring_offsets);

cpp/include/cuspatial/distance/linestring_polygon_distance.hpp

@@ -0,0 +1,49 @@
+/*
+ * 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.
+ */
+
+#pragma once
+
+#include <cuspatial/column/geometry_column_view.hpp>
+
+#include <cudf/column/column_view.hpp>
+
+#include <optional>
+
+namespace cuspatial {
+
+/**
+ * @ingroup distance
+ * @brief Compute pairwise (multi)linestring-to-(multi)polygon Cartesian distance
+ *
+ * @param multilinestrings Geometry column of multilinestrings
+ * @param multipolygons Geometry column of multipolygons
+ * @param mr Device memory resource used to allocate the returned column.
+ * @return Column of distances between each pair of input geometries, same type as input coordinate
+ * types.
+ *
+ * @throw cuspatial::logic_error if `multilinestrings` and `multipolygons` has different coordinate
+ * types.
+ * @throw cuspatial::logic_error if `multilinestrings` is not a linestring column and
+ * `multipolygons` is not a polygon column.
+ * @throw cuspatial::logic_error if input column sizes mismatch.
+ */
+
+std::unique_ptr<cudf::column> pairwise_linestring_polygon_distance(
+  geometry_column_view const& multilinestrings,
+  geometry_column_view const& multipolygons,
+  rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource());
+
+}  // namespace cuspatial

cpp/include/cuspatial/distance/polygon_distance.hpp

@@ -0,0 +1,40 @@
+/*
+ * 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 <cuspatial/column/geometry_column_view.hpp>
+
+#include <cudf/column/column.hpp>
+
+#include <rmm/mr/device/device_memory_resource.hpp>
+
+#include <memory>
+
+namespace cuspatial {
+
+/**
+ * @brief
+ *
+ * @param lhs
+ * @param rhs
+ * @param mr
+ * @return std::unique_ptr<cudf::column>
+ */
+std::unique_ptr<cudf::column> pairwise_polygon_distance(
+  geometry_column_view const& lhs,
+  geometry_column_view const& rhs,
+  rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource());
+
+}  // namespace cuspatial

cpp/include/cuspatial/error.hpp

@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2020-2022, 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.
@@ -16,6 +16,8 @@
 
 #pragma once
 
+#include <cuspatial/assert.cuh>
+
 #include <cuda_runtime_api.h>
 
 #include <stdexcept>
@@ -76,6 +78,31 @@ struct cuda_error : public std::runtime_error {
              : throw cuspatial::logic_error("cuSpatial failure at: " __FILE__ \
                                             ":" CUSPATIAL_STRINGIFY(__LINE__) ": " reason)
 
+/**---------------------------------------------------------------------------*
+ * @brief Macro for checking (pre-)conditions that throws an exception when
+ * a condition is violated.
+ *
+ * Example usage:
+ *
+ * @code
+ * CUSPATIAL_HOST_DEVICE_EXPECTS(lhs->dtype == rhs->dtype, "Column type mismatch");
+ * @endcode
+ *
+ * @param[in] cond Expression that evaluates to true or false
+ * @param[in] reason String literal description of the reason that cond is
+ * expected to be true
+ *
+ * (if on host)
+ * @throw cuspatial::logic_error if the condition evaluates to false.
+ * (if on device)
+ * program terminates and assertion error message is printed to stderr.
+ *---------------------------------------------------------------------------**/
+#ifndef __CUDA_ARCH__
+#define CUSPATIAL_HOST_DEVICE_EXPECTS(cond, reason) CUSPATIAL_EXPECTS(cond, reason)
+#else
+#define CUSPATIAL_HOST_DEVICE_EXPECTS(cond, reason) cuspatial_assert(cond&& reason)
+#endif
+
 /**---------------------------------------------------------------------------*
  * @brief Indicates that an erroneous code path has been taken.
  *

cpp/include/cuspatial/experimental/detail/algorithm/linestring_distance.cuh

@@ -0,0 +1,74 @@
+/*
+ * 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 <cuspatial/detail/utility/device_atomics.cuh>
+#include <cuspatial/detail/utility/linestring.cuh>
+
+#include <rmm/device_uvector.hpp>
+
+#include <thrust/optional.h>
+
+namespace cuspatial {
+namespace detail {
+
+/**
+ * @internal
+ * @brief The kernel to compute linestring to linestring distance
+ *
+ * Each thread of the kernel computes the distance between a segment in a linestring in pair 1 to a
+ * linestring in pair 2. For a segment in pair 1, the linestring index is looked up from the offset
+ * array and mapped to the linestring in the pair 2. The segment is then computed with all segments
+ * in the corresponding linestring in pair 2. This forms a local minima of the shortest distance,
+ * which is then combined with other segment results via an atomic operation to form the global
+ * minimum distance between the linestrings.
+ *
+ * `intersects` is an optional pointer to a boolean range where the `i`th element indicates the
+ * `i`th output should be set to 0 and bypass distance computation. This argument is optional, if
+ * set to nullopt, no distance computation will be bypassed.
+ */
+template <class MultiLinestringRange1, class MultiLinestringRange2, class OutputIt>
+__global__ void linestring_distance(MultiLinestringRange1 multilinestrings1,
+                                    MultiLinestringRange2 multilinestrings2,
+                                    thrust::optional<uint8_t*> intersects,
+                                    OutputIt distances_first)
+{
+  using T = typename MultiLinestringRange1::element_t;
+
+  for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < multilinestrings1.num_points();
+       idx += gridDim.x * blockDim.x) {
+    auto const part_idx = multilinestrings1.part_idx_from_point_idx(idx);
+    if (!multilinestrings1.is_valid_segment_id(idx, part_idx)) continue;
+    auto const geometry_idx = multilinestrings1.geometry_idx_from_part_idx(part_idx);
+
+    if (intersects.has_value() && intersects.value()[geometry_idx]) {
+      distances_first[geometry_idx] = 0;
+      continue;
+    }
+
+    auto [a, b]            = multilinestrings1.segment(idx);
+    T min_distance_squared = std::numeric_limits<T>::max();
+
+    for (auto const& linestring2 : multilinestrings2[geometry_idx]) {
+      for (auto [c, d] : linestring2) {
+        min_distance_squared = min(min_distance_squared, squared_segment_distance(a, b, c, d));
+      }
+    }
+    atomicMin(&distances_first[geometry_idx], static_cast<T>(sqrt(min_distance_squared)));
+  }
+}
+
+}  // namespace detail
+}  // namespace cuspatial

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

@@ -15,6 +15,7 @@
  */
 
 #include <cuspatial/detail/utility/zero_data.cuh>
+#include <cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh>
 #include <cuspatial/experimental/iterator_factory.cuh>
 #include <cuspatial/traits.hpp>
 #include <cuspatial/vec_2d.hpp>