Support Multilinestring in column API and python Bindings for `pair…

…wise_linestring_distance` (#786)

This is a port of the second half of #753, closes #706

Authors:
  - Michael Wang (https://github.com/isVoid)

Approvers:
  - Mark Harris (https://github.com/harrism)
  - H. Thomson Comer (https://github.com/thomcom)

URL: #786

cpp/include/cuspatial/distance/linestring_distance.hpp

@@ -25,23 +25,21 @@
 
 namespace cuspatial {
 
-/**
- * @addtogroup distance
- * @{
- */
-
 /**
  * @brief Compute shortest distance between pairs of linestrings
+ * @ingroup distance
  *
  * The shortest distance between two linestrings is defined as the shortest distance
  * between all pairs of segments of the two linestrings. If any of the segments intersect,
- * the distance is 0.
+ * the distance is 0. The shortest distance between two multilinestrings is defined as the
+ * the shortest distance between all pairs of linestrings of the two multilinestrings.
  *
- * The following example contains 4 pairs of linestrings.
+ * The following example contains 4 pairs of linestrings. The first array is a single linestring
+ * array and the second array is a multilinestring array.
  * ```
  * First pair:
  * (0, 1) -> (1, 0) -> (-1, 0)
- * (1, 1) -> (2, 1) -> (2, 0) -> (3, 0)
+ * {(1, 1) -> (2, 1) -> (2, 0) -> (3, 0)}
  *
  *     |
  *     *   #---#
@@ -57,65 +55,69 @@ namespace cuspatial {
  * Second pair:
  *
  * (0, 0) -> (0, 1)
- * (1, 0) -> (1, 1) -> (1, 2)
+ * {(1, 0) -> (1, 1) -> (1, 2), (1, -1) -> (1, -2) -> (1, -3)}
  *
- * These linestrings are parallel. Their distance is 1 (point (0, 0) to point (1, 0)).
+ * The linestrings in the pairs are parallel. Their distance is 1 (point (0, 0) to point (1, 0)).
  *
  * Third pair:
  *
  * (0, 0) -> (2, 2) -> (-2, 0)
- * (2, 0) -> (0, 2)
+ * {(2, 0) -> (0, 2), (0, 2) -> (-2, 0)}
  *
- * These linestrings intersect, so their distance is 0.
+ * The linestrings in the pairs intersect, so their distance is 0.
  *
  * Forth pair:
  *
  * (2, 2) -> (-2, -2)
- * (1, 1) -> (5, 5) -> (10, 0)
+ * {(1, 1) -> (5, 5) -> (10, 0), (-1, -1) -> (-5, -5) -> (-10, 0)}
  *
  * These linestrings contain colinear and overlapping sections, so
  * their distance is 0.
  *
  * The input of above example is:
- * linestring1_offsets:  {0, 3, 5, 8}
- * linestring1_points_x: {0, 1, -1, 0, 0, 0, 2, -2, 2, -2}
- * linestring1_points_y: {1, 0, 0, 0, 1, 0, 2, 0, 2, -2}
- * linestring2_offsets:  {0, 4, 7, 9}
- * linestring2_points_x: {1, 2, 2, 3, 1, 1, 1, 2, 0, 1, 5, 10}
- * linestring2_points_y: {1, 1, 0, 0, 0, 1, 2, 0, 2, 1, 5, 0}
+ * multilinestring1_geometry_offsets: nullopt
+ * linestring1_part_offsets:  {0, 3, 5, 8, 10}
+ * linestring1_points_xy:
+ * {0, 1, 1, 0, -1, 0, 0, 0, 0, 1, 0, 0, 2, 2, -2, 0, 2, 2, -2, -2}
+ *
+ * multilinestring2_geometry_offsets: {0, 1, 3, 5, 7}
+ * linestring2_offsets:  {0, 4, 7, 10, 12, 14, 17, 20}
+ * linestring2_points_xy: {1, 1, 2, 1, 2, 0, 3, 0, 1, 0, 1, 1, 1, 2, 1, -1, 1, -2, 1, -3, 2, 0, 0,
+ * 2, 0, 2, -2, 0, 1, 1, 5, 5, 10, 0, -1, -1, -5, -5, -10, 0}
  *
  * Result: {sqrt(2.0)/2, 1, 0, 0}
  * ```
  *
- * @param linestring1_offsets Indices of the first point of the first linestring of each pair
- * @param linestring1_points_x x-components of points in the first linestring of each pair
- * @param linestring1_points_y y-component of points in the first linestring of each pair
- * @param linestring2_offsets Indices of the first point of the second linestring of each pair
- * @param linestring2_points_x x-component of points in the second linestring of each pair
- * @param linestring2_points_y y-component of points in the second linestring of each pair
+ * @param multilinestring1_geometry_offsets Beginning and ending indices to each multilinestring in
+ * the first multilinestring array.
+ * @param linestring1_part_offsets Beginning and ending indices for each linestring in the point
+ * array. Because the coordinates are interleaved, the actual starting position for the coordinate
+ * of linestring `i` is `2*linestring_part_offsets[i]`.
+ * @param linestring1_points_xy Interleaved x, y-coordinates of linestring points.
+ * @param multilinestring2_geometry_offsets Beginning and ending indices to each multilinestring in
+ * the second multilinestring array.
+ * @param linestring2_part_offsets Beginning and ending indices for each linestring in the point
+ * array. Because the coordinates are interleaved, the actual starting position for the coordinate
+ * of linestring `i` is `2*linestring_part_offsets[i]`.
+ * @param linestring2_points_xy Interleaved x, y-coordinates of linestring points.
  * @param mr Device memory resource used to allocate the returned column's device memory
- * @return A column of shortest distances between each pair of linestrings
+ * @return A column of shortest distances between each pair of (multi)linestrings
  *
+ * @note If `multilinestring_geometry_offset` is std::nullopt, the input is a single linestring
+ * array.
  * @note If any of the linestring contains less than 2 points, the behavior is undefined.
  *
  * @throw cuspatial::logic_error if `linestring1_offsets.size() != linestring2_offsets.size()`
- * @throw cuspatial::logic_error if there is a size mismatch between the x- and y-coordinates of the
- * linestring points.
  * @throw cuspatial::logic_error if any of the point arrays have mismatched types.
  * @throw cuspatial::logic_error if any linestring has fewer than 2 points.
  *
  */
 std::unique_ptr<cudf::column> pairwise_linestring_distance(
-  cudf::device_span<cudf::size_type const> linestring1_offsets,
-  cudf::column_view const& linestring1_points_x,
-  cudf::column_view const& linestring1_points_y,
-  cudf::device_span<cudf::size_type const> linestring2_offsets,
-  cudf::column_view const& linestring2_points_x,
-  cudf::column_view const& linestring2_points_y,
+  std::optional<cudf::device_span<cudf::size_type const>> multilinestring1_geometry_offsets,
+  cudf::device_span<cudf::size_type const> linestring1_part_offsets,
+  cudf::column_view const& linestring1_points_xy,
+  std::optional<cudf::device_span<cudf::size_type const>> multilinestring2_geometry_offsets,
+  cudf::device_span<cudf::size_type const> linestring2_part_offsets,
+  cudf::column_view const& linestring2_points_xy,
   rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource());
-
-/**
- * @} // end of doxygen group
- */
-
 }  // namespace cuspatial

cpp/src/spatial/linestring_distance.cu

@@ -14,6 +14,9 @@
  * limitations under the License.
  */
 
+#include "../utility/double_boolean_dispatch.hpp"
+#include "../utility/iterator.hpp"
+
 #include <cuspatial/error.hpp>
 #include <cuspatial/experimental/iterator_factory.cuh>
 #include <cuspatial/experimental/linestring_distance.cuh>
@@ -23,59 +26,66 @@
 #include <cudf/column/column_factories.hpp>
 #include <cudf/column/column_view.hpp>
 #include <cudf/copying.hpp>
+#include <cudf/detail/utilities/device_atomics.cuh>
 #include <cudf/utilities/type_dispatcher.hpp>
 
 #include <rmm/cuda_stream_view.hpp>
 
-#include <thrust/iterator/counting_iterator.h>
-
 #include <limits>
 #include <memory>
 #include <type_traits>
-
 namespace cuspatial {
 namespace detail {
-struct pairwise_linestring_distance_functor {
+
+template <bool first_is_multilinestring, bool second_is_multilinestring>
+struct pairwise_linestring_distance_launch {
+  using SizeType = cudf::device_span<cudf::size_type const>::size_type;
+
   template <typename T>
   std::enable_if_t<std::is_floating_point<T>::value, std::unique_ptr<cudf::column>> operator()(
-    cudf::device_span<cudf::size_type const> linestring1_offsets,
-    cudf::column_view const& linestring1_points_x,
-    cudf::column_view const& linestring1_points_y,
-    cudf::device_span<cudf::size_type const> linestring2_offsets,
-    cudf::column_view const& linestring2_points_x,
-    cudf::column_view const& linestring2_points_y,
+    SizeType num_pairs,
+    std::optional<cudf::device_span<cudf::size_type const>> multilinestring1_geometry_offsets,
+    cudf::device_span<cudf::size_type const> linestring1_part_offsets,
+    cudf::column_view const& linestring1_points_xy,
+    std::optional<cudf::device_span<cudf::size_type const>> multilinestring2_geometry_offsets,
+    cudf::device_span<cudf::size_type const> linestring2_part_offsets,
+    cudf::column_view const& linestring2_points_xy,
     rmm::cuda_stream_view stream,
     rmm::mr::device_memory_resource* mr)
   {
-    auto const num_string_pairs = static_cast<cudf::size_type>(linestring1_offsets.size()) - 1;
-
-    auto distances = cudf::make_numeric_column(cudf::data_type{cudf::type_to_id<T>()},
-                                               num_string_pairs,
-                                               cudf::mask_state::UNALLOCATED,
-                                               stream,
-                                               mr);
-
-    auto linestring1_coords_it =
-      make_vec_2d_iterator(linestring1_points_x.begin<T>(), linestring1_points_y.begin<T>());
-    auto linestring2_coords_it =
-      make_vec_2d_iterator(linestring2_points_x.begin<T>(), linestring2_points_y.begin<T>());
-
-    auto multilinestrings1 = make_multilinestring_range(num_string_pairs,
-                                                        thrust::make_counting_iterator(0),
-                                                        num_string_pairs,
-                                                        linestring1_offsets.begin(),
-                                                        linestring1_points_x.size(),
-                                                        linestring1_coords_it);
-
-    auto multilinestrings2 = make_multilinestring_range(num_string_pairs,
-                                                        thrust::make_counting_iterator(0),
-                                                        num_string_pairs,
-                                                        linestring2_offsets.begin(),
-                                                        linestring2_points_x.size(),
-                                                        linestring2_coords_it);
+    auto const num_multilinestring1_parts =
+      static_cast<SizeType>(linestring1_part_offsets.size() - 1);
+    auto const num_multilinestring2_parts =
+      static_cast<SizeType>(linestring2_part_offsets.size() - 1);
+    auto const num_multilinestring1_points =
+      static_cast<SizeType>(linestring1_points_xy.size() / 2);
+    auto const num_multilinestring2_points =
+      static_cast<SizeType>(linestring2_points_xy.size() / 2);
+
+    auto distances = cudf::make_numeric_column(
+      cudf::data_type{cudf::type_to_id<T>()}, num_pairs, cudf::mask_state::UNALLOCATED, stream, mr);
+
+    auto linestring1_coords_it = make_vec_2d_iterator(linestring1_points_xy.begin<T>());
+    auto linestring2_coords_it = make_vec_2d_iterator(linestring2_points_xy.begin<T>());
+
+    auto multilinestrings1 = make_multilinestring_range(
+      num_pairs,
+      get_geometry_iterator_functor<first_is_multilinestring>{}(multilinestring1_geometry_offsets),
+      num_multilinestring1_parts,
+      linestring1_part_offsets.begin(),
+      num_multilinestring1_points,
+      linestring1_coords_it);
+
+    auto multilinestrings2 = make_multilinestring_range(
+      num_pairs,
+      get_geometry_iterator_functor<second_is_multilinestring>{}(multilinestring2_geometry_offsets),
+      num_multilinestring2_parts,
+      linestring2_part_offsets.begin(),
+      num_multilinestring2_points,
+      linestring2_coords_it);
 
     pairwise_linestring_distance(
-      multilinestrings1, multilinestrings2, distances->mutable_view().begin<T>());
+      multilinestrings1, multilinestrings2, distances->mutable_view().begin<T>(), stream);
 
     return distances;
   }
@@ -88,61 +98,72 @@ struct pairwise_linestring_distance_functor {
   }
 };
 
-std::unique_ptr<cudf::column> pairwise_linestring_distance(
-  cudf::device_span<cudf::size_type const> linestring1_offsets,
-  cudf::column_view const& linestring1_points_x,
-  cudf::column_view const& linestring1_points_y,
-  cudf::device_span<cudf::size_type const> linestring2_offsets,
-  cudf::column_view const& linestring2_points_x,
-  cudf::column_view const& linestring2_points_y,
-  rmm::cuda_stream_view stream,
-  rmm::mr::device_memory_resource* mr)
-{
-  CUSPATIAL_EXPECTS(linestring1_offsets.size() == linestring2_offsets.size(),
-                    "Mismatch number of linestrings in the linestring pair array.");
-
-  CUSPATIAL_EXPECTS(linestring1_points_x.size() == linestring1_points_y.size() and
-                      linestring2_points_x.size() == linestring2_points_y.size(),
-                    "The lengths of linestring coordinates arrays mismatch.");
-
-  CUSPATIAL_EXPECTS(linestring1_points_x.type() == linestring1_points_y.type() and
-                      linestring2_points_x.type() == linestring2_points_y.type() and
-                      linestring1_points_x.type() == linestring2_points_x.type(),
-                    "The types of linestring coordinates arrays mismatch.");
-
-  if (linestring1_offsets.size() - 1 == 0) { return cudf::empty_like(linestring1_points_x); }
-
-  return cudf::type_dispatcher(linestring1_points_x.type(),
-                               pairwise_linestring_distance_functor{},
-                               linestring1_offsets,
-                               linestring1_points_x,
-                               linestring1_points_y,
-                               linestring2_offsets,
-                               linestring2_points_x,
-                               linestring2_points_y,
-                               stream,
-                               mr);
-}
-
+template <bool first_is_multilinestring, bool second_is_multilinestring>
+struct pairwise_linestring_distance_functor {
+  std::unique_ptr<cudf::column> operator()(
+    std::optional<cudf::device_span<cudf::size_type const>> multilinestring1_geometry_offsets,
+    cudf::device_span<cudf::size_type const> linestring1_part_offsets,
+    cudf::column_view const& linestring1_points_xy,
+    std::optional<cudf::device_span<cudf::size_type const>> multilinestring2_geometry_offsets,
+    cudf::device_span<cudf::size_type const> linestring2_part_offsets,
+    cudf::column_view const& linestring2_points_xy,
+    rmm::cuda_stream_view stream,
+    rmm::mr::device_memory_resource* mr)
+  {
+    CUSPATIAL_EXPECTS(
+      linestring1_points_xy.size() % 2 == 0 && linestring2_points_xy.size() % 2 == 0,
+      "Points array must contain even number of coordinates.");
+
+    auto num_lhs = first_is_multilinestring ? multilinestring1_geometry_offsets.value().size() - 1
+                                            : linestring1_part_offsets.size() - 1;
+    auto num_rhs = second_is_multilinestring ? multilinestring2_geometry_offsets.value().size() - 1
+                                             : linestring2_part_offsets.size() - 1;
+
+    CUSPATIAL_EXPECTS(num_lhs == num_rhs, "Mismatch number of points and linestrings.");
+
+    CUSPATIAL_EXPECTS(linestring1_points_xy.type() == linestring2_points_xy.type(),
+                      "The types of linestring coordinates arrays mismatch.");
+
+    CUSPATIAL_EXPECTS(!(linestring1_points_xy.has_nulls() || linestring2_points_xy.has_nulls()),
+                      "All inputs must not have nulls.");
+
+    if (num_lhs == 0) { return cudf::empty_like(linestring1_points_xy); }
+
+    return cudf::type_dispatcher(
+      linestring1_points_xy.type(),
+      pairwise_linestring_distance_launch<first_is_multilinestring, second_is_multilinestring>{},
+      num_lhs,
+      multilinestring1_geometry_offsets,
+      linestring1_part_offsets,
+      linestring1_points_xy,
+      multilinestring2_geometry_offsets,
+      linestring2_part_offsets,
+      linestring2_points_xy,
+      stream,
+      mr);
+  }
+};
 }  // namespace detail
-
 std::unique_ptr<cudf::column> pairwise_linestring_distance(
-  cudf::device_span<cudf::size_type const> linestring1_offsets,
-  cudf::column_view const& linestring1_points_x,
-  cudf::column_view const& linestring1_points_y,
-  cudf::device_span<cudf::size_type const> linestring2_offsets,
-  cudf::column_view const& linestring2_points_x,
-  cudf::column_view const& linestring2_points_y,
+  std::optional<cudf::device_span<cudf::size_type const>> multilinestring1_geometry_offsets,
+  cudf::device_span<cudf::size_type const> linestring1_part_offsets,
+  cudf::column_view const& linestring1_points_xy,
+  std::optional<cudf::device_span<cudf::size_type const>> multilinestring2_geometry_offsets,
+  cudf::device_span<cudf::size_type const> linestring2_part_offsets,
+  cudf::column_view const& linestring2_points_xy,
   rmm::mr::device_memory_resource* mr)
 {
-  return detail::pairwise_linestring_distance(linestring1_offsets,
-                                              linestring1_points_x,
-                                              linestring1_points_y,
-                                              linestring2_offsets,
-                                              linestring2_points_x,
-                                              linestring2_points_y,
-                                              rmm::cuda_stream_default,
-                                              mr);
+  return double_boolean_dispatch<detail::pairwise_linestring_distance_functor>(
+    multilinestring1_geometry_offsets.has_value(),
+    multilinestring2_geometry_offsets.has_value(),
+    multilinestring1_geometry_offsets,
+    linestring1_part_offsets,
+    linestring1_points_xy,
+    multilinestring2_geometry_offsets,
+    linestring2_part_offsets,
+    linestring2_points_xy,
+    rmm::cuda_stream_default,
+    mr);
 }
 
 }  // namespace cuspatial

cpp/tests/CMakeLists.txt

@@ -84,7 +84,7 @@ ConfigureTest(POINT_LINESTRING_DISTANCE_TEST
     spatial/point_linestring_distance_test.cpp)
 
 ConfigureTest(LINESTRING_DISTANCE_TEST
-    spatial/linestring_distance_test.cpp spatial/linestring_distance_test_medium.cpp)
+    spatial/linestring_distance_test.cpp)
 
 ConfigureTest(POINT_LINESTRING_NEAREST_POINT_TEST
     spatial/point_linestring_nearest_points_test.cpp)
@@ -110,9 +110,6 @@ ConfigureTest(TRAJECTORY_BOUNDING_BOXES_TEST
 ConfigureTest(SPATIAL_WINDOW_POINT_TEST
     spatial_window/spatial_window_test.cpp)
 
-ConfigureTest(FLOAT_EQUIVALENT_UTILITY_TEST
-    utility_test/test_float_equivalent.cu)
-
 # Experimental API
 ConfigureTest(HAVERSINE_TEST_EXP
     experimental/spatial/haversine_test.cu)
@@ -126,7 +123,8 @@ ConfigureTest(HAUSDORFF_TEST_EXP
     experimental/spatial/hausdorff_test.cu)
 
 ConfigureTest(LINESTRING_DISTANCE_TEST_EXP
-    experimental/spatial/linestring_distance_test.cu)
+    experimental/spatial/linestring_distance_test.cu
+    experimental/spatial/linestring_distance_test_medium.cu)
 
 ConfigureTest(POINT_LINESTRING_NEAREST_POINT_TEST_EXP
     experimental/spatial/point_linestring_nearest_points_test.cu)