Refactor haversine_distance to a header-only API

CHANGELOG.md

@@ -1,3 +1,7 @@
+# cuSpatial 22.04.00 (Date TBD)
+
+Please see https://github.com/rapidsai/cuspatial/releases/tag/v22.04.00a for the latest changes to this development branch.
+
 # cuSpatial 22.02.00 (Date TBD)
 
 Please see https://github.com/rapidsai/cuspatial/releases/tag/v22.02.00a for the latest changes to this development branch.

conda/environments/cuspatial_dev_cuda11.0.yml

@@ -7,7 +7,7 @@ channels:
 dependencies:
   - clang=11.1.0
   - clang-tools=11.1.0
-  - cudf=22.02.*
+  - cudf=22.04.*
   - cudatoolkit=11.0
   - gdal>=3.0.2
   - geopandas>=0.9.0,<0.10.0a0

conda/environments/cuspatial_dev_cuda11.1.yml

@@ -7,7 +7,7 @@ channels:
 dependencies:
   - clang=11.1.0
   - clang-tools=11.1.0
-  - cudf=22.02.*
+  - cudf=22.04.*
   - cudatoolkit=11.1
   - gdal>=3.0.2
   - geopandas>=0.9.0,<0.10.0a0

conda/environments/cuspatial_dev_cuda11.2.yml

@@ -7,7 +7,7 @@ channels:
 dependencies:
   - clang=11.1.0
   - clang-tools=11.1.0
-  - cudf=22.02.*
+  - cudf=22.04.*
   - cudatoolkit=11.2
   - gdal>=3.0.2
   - geopandas>=0.9.0,<0.10.0a0

cpp/CMakeLists.txt

@@ -28,7 +28,7 @@ elseif(CMAKE_CUDA_ARCHITECTURES STREQUAL "")
   set(CUSPATIAL_BUILD_FOR_DETECTED_ARCHS TRUE)
 endif()
 
-project(CUSPATIAL VERSION 22.02.00 LANGUAGES C CXX)
+project(CUSPATIAL VERSION 22.04.00 LANGUAGES C CXX)
 
 # Needed because GoogleBenchmark changes the state of FindThreads.cmake,
 # causing subsequent runs to have different values for the `Threads::Threads` target.

cpp/include/cuspatial/experimental/detail/haversine.hpp

@@ -0,0 +1,125 @@
+/*
+ * 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.
+ */
+
+#pragma once
+
+#include <cuspatial/constants.hpp>
+#include <cuspatial/error.hpp>
+
+#include <rmm/cuda_stream_view.hpp>
+#include <rmm/exec_policy.hpp>
+#include <rmm/mr/device/device_memory_resource.hpp>
+#include <rmm/mr/device/per_device_resource.hpp>
+
+#include <thrust/iterator/constant_iterator.h>
+#include <thrust/iterator/zip_iterator.h>
+#include <thrust/transform.h>
+#include <thrust/tuple.h>
+
+#include <type_traits>
+
+namespace cuspatial {
+
+namespace detail {
+
+template <typename T>
+__device__ T calculate_haversine_distance(T radius, T a_lon, T a_lat, T b_lon, T b_lat)
+{
+  auto ax = a_lon * DEGREE_TO_RADIAN;
+  auto ay = a_lat * DEGREE_TO_RADIAN;
+  auto bx = b_lon * DEGREE_TO_RADIAN;
+  auto by = b_lat * DEGREE_TO_RADIAN;
+
+  // haversine formula
+  auto x        = (bx - ax) / 2;
+  auto y        = (by - ay) / 2;
+  auto sinysqrd = sin(y) * sin(y);
+  auto sinxsqrd = sin(x) * sin(x);
+  auto scale    = cos(ay) * cos(by);
+
+  return 2 * radius * asin(sqrt(sinysqrd + sinxsqrd * scale));
+};
+
+template <class LonItA,
+          class LatItA,
+          class LonItB,
+          class LatItB,
+          class OutputIt,
+          class T = typename std::iterator_traits<LonItA>::value_type>
+OutputIt haversine_distance(LonItA a_lon_first,
+                            LonItA a_lon_last,
+                            LatItA a_lat_first,
+                            LonItB b_lon_first,
+                            LatItB b_lat_first,
+                            OutputIt distance_first,
+                            T const radius               = EARTH_RADIUS_KM,
+                            rmm::cuda_stream_view stream = rmm::cuda_stream_default)
+{
+  static_assert(
+    std::conjunction_v<std::is_floating_point<T>,
+                       std::is_floating_point<typename std::iterator_traits<LonItA>::value_type>,
+                       std::is_floating_point<typename std::iterator_traits<LatItA>::value_type>,
+                       std::is_floating_point<typename std::iterator_traits<LonItB>::value_type>,
+                       std::is_floating_point<typename std::iterator_traits<LatItB>::value_type>,
+                       std::is_floating_point<typename std::iterator_traits<OutputIt>::value_type>>,
+    "Haversine distance supports only floating-point coordinates.");
+
+  CUSPATIAL_EXPECTS(radius > 0, "radius must be positive.");
+
+  auto input_tuple = thrust::make_tuple(thrust::make_constant_iterator(static_cast<T>(radius)),
+                                        a_lon_first,
+                                        a_lat_first,
+                                        b_lon_first,
+                                        b_lat_first);
+
+  auto input_iter = thrust::make_zip_iterator(input_tuple);
+
+  auto output_size = std::distance(a_lon_first, a_lon_last);
+
+  return thrust::transform(rmm::exec_policy(stream),
+                           input_iter,
+                           input_iter + output_size,
+                           distance_first,
+                           [] __device__(auto inputs) {
+                             return calculate_haversine_distance(thrust::get<0>(inputs),
+                                                                 thrust::get<1>(inputs),
+                                                                 thrust::get<2>(inputs),
+                                                                 thrust::get<3>(inputs),
+                                                                 thrust::get<4>(inputs));
+                           });
+}
+}  // namespace detail
+
+template <class LonItA, class LatItA, class LonItB, class LatItB, class OutputIt, class T>
+OutputIt haversine_distance(LonItA a_lon_first,
+                            LonItA a_lon_last,
+                            LatItA a_lat_first,
+                            LonItB b_lon_first,
+                            LatItB b_lat_first,
+                            OutputIt distance_first,
+                            T const radius)
+{
+  return detail::haversine_distance(a_lon_first,
+                                    a_lon_last,
+                                    a_lat_first,
+                                    b_lon_first,
+                                    b_lat_first,
+                                    distance_first,
+                                    radius,
+                                    rmm::cuda_stream_default);
+}
+
+}  // namespace cuspatial

cpp/include/cuspatial/experimental/haversine.hpp

@@ -0,0 +1,75 @@
+/*
+ * 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.
+ */
+
+#pragma once
+
+#include <cuspatial/constants.hpp>
+
+#include <rmm/mr/device/device_memory_resource.hpp>
+#include <rmm/mr/device/per_device_resource.hpp>
+
+namespace cuspatial {
+
+/**
+ * brief Compute haversine distances between points in set A to the corresponding points in set B.
+ *
+ * Computes N haversine distances, where N is `std::distance(a_lon_first, a_lon_last)`. One distance
+ * is computed for each (a_lon[i], a_lat[i]) and (b_lon[i], b_lat[i]) point pair. `distance_first`
+ * must be an iterator to output storage allocated for N distances.
+ *
+ * All iterators must have the same floating-point `value_type`.
+ *
+ * https://en.wikipedia.org/wiki/Haversine_formula
+ *
+ * @param[in]  a_lon_first: beginning of range of longitude of points in set A
+ * @param[in]  a_lon_last: end of the range of longitude of points in set A
+ * @param[in]  b_lon_first: beginning of range of longitude of points in set B
+ * @param[in]  b_lat_first: beginning of range of latitude of points in set B
+ * @param[out] distance_first: beginning of output range of haversine distances
+ * @param[in] radius: radius of the sphere on which the points reside. default: 6371.0
+ *            (approximate radius of Earth in km)
+ *
+ * @tparam LonItA must meet the requirements of [LegacyForwardIterator][LinkLFI] and be
+ * device-accessible.
+ * @tparam LatItA must meet the requirements of [LegacyForwardIterator][LinkLFI] and be
+ * device-accessible.
+ * @tparam LonItB must meet the requirements of [LegacyForwardIterator][LinkLFI] and be
+ * device-accessible.
+ * @tparam LatItB must meet the requirements of [LegacyForwardIterator][LinkLFI] and be
+ * device-accessible.
+ * @tparam OutputIt must meet the requirements of [LegacyForwardIterator][LinkLFI] and be
+ * device-accessible.
+ *
+ * @return Output iterator to the element past the last distance computed.
+ *
+ * [LinkLFI]: https://en.cppreference.com/w/cpp/named_req/ForwardIterator "LegacyForwardIterator"
+ */
+template <class LonItA,
+          class LatItA,
+          class LonItB,
+          class LatItB,
+          class OutputIt,
+          class T = typename std::iterator_traits<LonItA>::value_type>
+OutputIt haversine_distance(LonItA a_lon_first,
+                            LonItA a_lon_last,
+                            LatItA a_lat_first,
+                            LonItB b_lon_first,
+                            LatItB b_lat_first,
+                            OutputIt distance_first,
+                            T const radius = EARTH_RADIUS_KM);
+}  // namespace cuspatial
+
+#include <cuspatial/experimental/detail/haversine.hpp>

cpp/include/cuspatial/haversine.hpp

@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2020, NVIDIA CORPORATION.
+ * 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.
@@ -44,5 +44,4 @@ std::unique_ptr<cudf::column> haversine_distance(
   cudf::column_view const& b_lat,
   double const radius                 = EARTH_RADIUS_KM,
   rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource());
-
 }  // namespace cuspatial

cpp/src/spatial/haversine.cu

@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2020, NVIDIA CORPORATION.
+ * 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.
@@ -16,41 +16,22 @@
 
 #include <cuspatial/constants.hpp>
 #include <cuspatial/error.hpp>
+#include <cuspatial/experimental/haversine.hpp>
 
 #include <cudf/column/column.hpp>
-#include <cudf/column/column_factories.hpp>
 #include <cudf/column/column_view.hpp>
 #include <cudf/copying.hpp>
 #include <cudf/types.hpp>
 #include <cudf/utilities/type_dispatcher.hpp>
 
 #include <rmm/cuda_stream_view.hpp>
-#include <rmm/exec_policy.hpp>
 #include <rmm/mr/device/device_memory_resource.hpp>
 
 #include <memory>
 #include <type_traits>
 
 namespace {
 
-template <typename T>
-__device__ T calculate_haversine_distance(T radius, T a_lon, T a_lat, T b_lon, T b_lat)
-{
-  auto ax = a_lon * DEGREE_TO_RADIAN;
-  auto ay = a_lat * DEGREE_TO_RADIAN;
-  auto bx = b_lon * DEGREE_TO_RADIAN;
-  auto by = b_lat * DEGREE_TO_RADIAN;
-
-  // haversine formula
-  auto x        = (bx - ax) / 2;
-  auto y        = (by - ay) / 2;
-  auto sinysqrd = sin(y) * sin(y);
-  auto sinxsqrd = sin(x) * sin(x);
-  auto scale    = cos(ay) * cos(by);
-
-  return 2 * radius * asin(sqrt(sinysqrd + sinxsqrd * scale));
-};
-
 struct haversine_functor {
   template <typename T, typename... Args>
   std::enable_if_t<not std::is_floating_point<T>::value, std::unique_ptr<cudf::column>> operator()(
@@ -65,7 +46,7 @@ struct haversine_functor {
     cudf::column_view const& a_lat,
     cudf::column_view const& b_lon,
     cudf::column_view const& b_lat,
-    double radius,
+    T radius,
     rmm::cuda_stream_view stream,
     rmm::mr::device_memory_resource* mr)
   {
@@ -74,25 +55,15 @@ struct haversine_functor {
     auto mask_policy = cudf::mask_allocation_policy::NEVER;
     auto result      = cudf::allocate_like(a_lon, a_lon.size(), mask_policy, mr);
 
-    auto input_tuple = thrust::make_tuple(thrust::make_constant_iterator(static_cast<T>(radius)),
-                                          a_lon.begin<T>(),
-                                          a_lat.begin<T>(),
-                                          b_lon.begin<T>(),
-                                          b_lat.begin<T>());
-
-    auto input_iter = thrust::make_zip_iterator(input_tuple);
-
-    thrust::transform(rmm::exec_policy(stream),
-                      input_iter,
-                      input_iter + result->size(),
-                      result->mutable_view().begin<T>(),
-                      [] __device__(auto inputs) {
-                        return calculate_haversine_distance(thrust::get<0>(inputs),
-                                                            thrust::get<1>(inputs),
-                                                            thrust::get<2>(inputs),
-                                                            thrust::get<3>(inputs),
-                                                            thrust::get<4>(inputs));
-                      });
+    cuspatial::detail::haversine_distance(
+      a_lon.begin<T>(),
+      a_lon.end<T>(),
+      a_lat.begin<T>(),
+      b_lon.begin<T>(),
+      b_lat.begin<T>(),
+      static_cast<cudf::mutable_column_view>(*result).begin<T>(),
+      T{radius},
+      stream);
 
     return result;
   }

cpp/tests/CMakeLists.txt

@@ -1,5 +1,5 @@
 #=============================================================================
-# Copyright (c) 2019-2021, NVIDIA CORPORATION.
+# Copyright (c) 2019-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.
@@ -88,3 +88,7 @@ ConfigureTest(TRAJECTORY_BOUNDING_BOXES_TEST
 
 ConfigureTest(SPATIAL_WINDOW_POINT_TEST
     spatial_window/spatial_window_test.cpp)
+
+# Experimental API
+ConfigureTest(HAVERSINE_TEST_EXP
+    experimental/spatial/haversine_test.cu)