diff --git a/cpp/benchmarks/pairwise_linestring_distance.cu b/cpp/benchmarks/pairwise_linestring_distance.cu index 4e059355b..e70c72cee 100644 --- a/cpp/benchmarks/pairwise_linestring_distance.cu +++ b/cpp/benchmarks/pairwise_linestring_distance.cu @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -104,13 +105,21 @@ void pairwise_linestring_distance_benchmark(nvbench::state& state, nvbench::type auto [ls2, ls2_offset] = generate_linestring(num_string_pairs, num_segments_per_string, 1, {100, 100}); - auto ls1_offset_begin = ls1_offset.begin(); - auto ls2_offset_begin = ls2_offset.begin(); - auto ls1_points_begin = ls1.begin(); - auto ls2_points_begin = ls2.begin(); - auto distances = rmm::device_vector(ls1.size()); - auto out_it = distances.begin(); - + auto distances = rmm::device_vector(ls1.size()); + auto out_it = distances.begin(); + + auto multilinestrings1 = make_multilinestring_range(1, + thrust::make_counting_iterator(0), + num_string_pairs, + ls1_offset.begin(), + ls1.size(), + ls1.begin()); + auto multilinestrings2 = make_multilinestring_range(1, + thrust::make_counting_iterator(0), + num_string_pairs, + ls2_offset.begin(), + ls2.size(), + ls2.begin()); auto const total_points = ls1.size() + ls2.size(); state.add_element_count(num_string_pairs, "LineStringPairs"); @@ -120,22 +129,8 @@ void pairwise_linestring_distance_benchmark(nvbench::state& state, nvbench::type state.add_global_memory_writes(num_string_pairs); state.exec(nvbench::exec_tag::sync, - [&ls1_offset_begin, - &num_string_pairs, - &ls1_points_begin, - ls1_size = ls1.size(), - &ls2_offset_begin, - &ls2_points_begin, - ls2_size = ls2.size(), - &out_it](nvbench::launch& launch) { - pairwise_linestring_distance(ls1_offset_begin, - ls1_offset_begin + num_string_pairs, - ls1_points_begin, - ls1_points_begin + ls1_size, - ls2_offset_begin, - ls2_points_begin, - ls2_points_begin + ls2_size, - out_it); + [&multilinestrings1, &multilinestrings2, &out_it](nvbench::launch& launch) { + pairwise_linestring_distance(multilinestrings1, multilinestrings2, out_it); }); } diff --git a/cpp/include/cuspatial/detail/iterator.hpp b/cpp/include/cuspatial/detail/iterator.hpp index 677ed14ce..05b8a7aa2 100644 --- a/cpp/include/cuspatial/detail/iterator.hpp +++ b/cpp/include/cuspatial/detail/iterator.hpp @@ -15,8 +15,7 @@ */ #pragma once - -#include +#include #include #include @@ -24,8 +23,8 @@ namespace cuspatial { namespace detail { -template -inline auto make_counting_transform_iterator(cudf::size_type start, UnaryFunction f) +template +inline CUSPATIAL_HOST_DEVICE auto make_counting_transform_iterator(IndexType start, UnaryFunction f) { return thrust::make_transform_iterator(thrust::make_counting_iterator(start), f); } diff --git a/cpp/include/cuspatial/detail/utility/linestring.cuh b/cpp/include/cuspatial/detail/utility/linestring.cuh index fec4b9cc2..32fca397c 100644 --- a/cpp/include/cuspatial/detail/utility/linestring.cuh +++ b/cpp/include/cuspatial/detail/utility/linestring.cuh @@ -92,10 +92,9 @@ __forceinline__ T __device__ segment_distance_no_intersect_or_colinear(vec_2d vec_2d const& c, vec_2d const& d) { - auto dist_sqr = std::min(std::min(point_to_segment_distance_squared(a, c, d), - point_to_segment_distance_squared(b, c, d)), - std::min(point_to_segment_distance_squared(c, a, b), - point_to_segment_distance_squared(d, a, b))); + auto dist_sqr = min( + min(point_to_segment_distance_squared(a, c, d), point_to_segment_distance_squared(b, c, d)), + min(point_to_segment_distance_squared(c, a, b), point_to_segment_distance_squared(d, a, b))); return dist_sqr; } diff --git a/cpp/include/cuspatial/experimental/detail/geometry/linestring_ref.cuh b/cpp/include/cuspatial/experimental/detail/geometry/linestring_ref.cuh new file mode 100644 index 000000000..ebef9a460 --- /dev/null +++ b/cpp/include/cuspatial/experimental/detail/geometry/linestring_ref.cuh @@ -0,0 +1,71 @@ +/* + * 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 +#include +#include + +#include +#include + +namespace cuspatial { + +template +struct to_segment_functor { + using element_t = iterator_vec_base_type; + using difference_type = typename thrust::iterator_difference::type; + VecIterator _point_begin; + + CUSPATIAL_HOST_DEVICE + to_segment_functor(VecIterator point_begin) : _point_begin(point_begin) {} + + CUSPATIAL_HOST_DEVICE + thrust::pair, vec_2d> operator()(difference_type i) + { + return {_point_begin[i], _point_begin[i + 1]}; + } +}; + +template +CUSPATIAL_HOST_DEVICE linestring_ref::linestring_ref(VecIterator begin, + VecIterator end) + : _point_begin(begin), _point_end(end) +{ + using T = iterator_vec_base_type; + static_assert(is_same, iterator_value_type>(), "must be vec2d type"); +} + +template +CUSPATIAL_HOST_DEVICE auto linestring_ref::num_segments() const +{ + // The number of segment equals the number of points minus 1. And the number of points + // is thrust::distance(_point_begin, _point_end) - 1. + return thrust::distance(_point_begin, _point_end) - 1; +} + +template +CUSPATIAL_HOST_DEVICE auto linestring_ref::segment_begin() const +{ + return detail::make_counting_transform_iterator(0, to_segment_functor{_point_begin}); +} + +template +CUSPATIAL_HOST_DEVICE auto linestring_ref::segment_end() const +{ + return segment_begin() + num_segments(); +} + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/detail/geometry_collection/multilinestring_ref.cuh b/cpp/include/cuspatial/experimental/detail/geometry_collection/multilinestring_ref.cuh new file mode 100644 index 000000000..6ca4fe4a6 --- /dev/null +++ b/cpp/include/cuspatial/experimental/detail/geometry_collection/multilinestring_ref.cuh @@ -0,0 +1,70 @@ +#pragma once + +#include "cuspatial/cuda_utils.hpp" +#include +#include + +#include + +namespace cuspatial { + +template +struct to_linestring_functor { + using difference_type = typename thrust::iterator_difference::type; + PartIterator part_begin; + VecIterator point_begin; + + CUSPATIAL_HOST_DEVICE + to_linestring_functor(PartIterator part_begin, VecIterator point_begin) + : part_begin(part_begin), point_begin(point_begin) + { + } + + CUSPATIAL_HOST_DEVICE auto operator()(difference_type i) + { + return linestring_ref{point_begin + part_begin[i], point_begin + part_begin[i + 1]}; + } +}; + +template +class multilinestring_ref; + +template +CUSPATIAL_HOST_DEVICE multilinestring_ref::multilinestring_ref( + PartIterator part_begin, PartIterator part_end, VecIterator point_begin, VecIterator point_end) + : _part_begin(part_begin), _part_end(part_end), _point_begin(point_begin), _point_end(point_end) +{ +} + +template +CUSPATIAL_HOST_DEVICE auto multilinestring_ref::num_linestrings() const +{ + return thrust::distance(_part_begin, _part_end) - 1; +} + +template +CUSPATIAL_HOST_DEVICE auto multilinestring_ref::part_begin() const +{ + return detail::make_counting_transform_iterator(0, + to_linestring_functor{_part_begin, _point_begin}); +} + +template +CUSPATIAL_HOST_DEVICE auto multilinestring_ref::part_end() const +{ + return part_begin() + num_linestrings(); +} + +template +CUSPATIAL_HOST_DEVICE auto multilinestring_ref::point_begin() const +{ + return _point_begin; +} + +template +CUSPATIAL_HOST_DEVICE auto multilinestring_ref::point_end() const +{ + return _point_end; +} + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/detail/linestring_distance.cuh b/cpp/include/cuspatial/experimental/detail/linestring_distance.cuh index 575115838..97648d74b 100644 --- a/cpp/include/cuspatial/experimental/detail/linestring_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/linestring_distance.cuh @@ -25,17 +25,9 @@ #include #include -#include -#include -#include #include -#include -#include -#include -#include #include -#include #include namespace cuspatial { @@ -43,7 +35,7 @@ namespace detail { /** * @internal - * @brief The kernel to compute point to linestring distance + * @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 @@ -51,142 +43,66 @@ namespace detail { * 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 globally * minimum distance between the linestrings. - * - * @tparam Cart2dItA Iterator to 2d cartesian coordinates. Must meet requirements of - * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam Cart2dItB Iterator to 2d cartesian coordinates. Must meet requirements of - * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam OffsetIterator Iterator to linestring offsets. Must meet requirements of - * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam OutputIterator Iterator to output distances. Must meet requirements of - * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible and mutable. - * - * @param[in] linestring1_offsets_begin Iterator to the begin of the range of linestring offsets in - * pair 1. - * @param[in] linestring1_offsets_end Iterator to the end of the range of linestring offsets in - * pair 1. - * @param[in] linestring1_points_xs_begin Iterator to the begin of the range of x coordinates of - * points in pair 1. - * @param[in] linestring1_points_xs_end Iterator to the end of the range of x coordinates of points - * in pair 1. - * @param[in] linestring2_offsets_begin Iterator to the begin of the range of linestring offsets - * in pair 2. - * @param[in] linestring2_points_xs_begin Iterator to the begin of the range of x coordinates of - * points in pair 2. - * @param[in] linestring2_points_xs_end Iterator to the end of the range of x coordinates of points - * in pair 2. - * @param[out] distances Iterator to the output range of shortest distances between pairs. - * - * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator - * "LegacyRandomAccessIterator" */ -template -void __global__ pairwise_linestring_distance_kernel(OffsetIterator linestring1_offsets_begin, - OffsetIterator linestring1_offsets_end, - Cart2dItA linestring1_points_begin, - Cart2dItA linestring1_points_end, - OffsetIterator linestring2_offsets_begin, - Cart2dItB linestring2_points_begin, - Cart2dItB linestring2_points_end, - OutputIterator distances) +template +__global__ void pairwise_linestring_distance_kernel(MultiLinestringRange1 multilinestrings1, + MultiLinestringRange2 multilinestrings2, + OutputIt distances_first) { - using T = typename std::iterator_traits::value_type::value_type; - - auto const p1_idx = threadIdx.x + blockIdx.x * blockDim.x; - std::size_t const num_linestrings = - thrust::distance(linestring1_offsets_begin, linestring1_offsets_end); - - std::size_t const linestring1_num_points = - thrust::distance(linestring1_points_begin, linestring1_points_end); - std::size_t const linestring2_num_points = - thrust::distance(linestring2_points_begin, linestring2_points_end); - - if (p1_idx >= linestring1_num_points) { return; } - - auto linestring_it = - thrust::upper_bound(thrust::seq, linestring1_offsets_begin, linestring1_offsets_end, p1_idx); - std::size_t const linestring_idx = - thrust::distance(linestring1_offsets_begin, thrust::prev(linestring_it)); - - auto const ls1_end = endpoint_index_of_linestring( - linestring_idx, linestring1_offsets_begin, num_linestrings, linestring1_num_points); - - if (p1_idx == ls1_end) { - // Current point is the end point of the line string. - return; + 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); + auto [a, b] = multilinestrings1.segment(idx); + T min_distance_squared = std::numeric_limits::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(sqrt(min_distance_squared))); } - - auto const ls2_start = *(linestring2_offsets_begin + linestring_idx); - auto const ls2_end = endpoint_index_of_linestring( - linestring_idx, linestring2_offsets_begin, num_linestrings, linestring2_num_points); - - auto const& A = thrust::raw_reference_cast(linestring1_points_begin[p1_idx]); - auto const& B = thrust::raw_reference_cast(linestring1_points_begin[p1_idx + 1]); - - auto min_squared_distance = std::numeric_limits::max(); - for (auto p2_idx = ls2_start; p2_idx < ls2_end; p2_idx++) { - auto const& C = thrust::raw_reference_cast(linestring2_points_begin[p2_idx]); - auto const& D = thrust::raw_reference_cast(linestring2_points_begin[p2_idx + 1]); - min_squared_distance = std::min(min_squared_distance, squared_segment_distance(A, B, C, D)); - } - - atomicMin(&thrust::raw_reference_cast(*(distances + linestring_idx)), - static_cast(std::sqrt(min_squared_distance))); } } // namespace detail -template -OutputIt pairwise_linestring_distance(OffsetIterator linestring1_offsets_first, - OffsetIterator linestring1_offsets_last, - Cart2dItA linestring1_points_first, - Cart2dItA linestring1_points_last, - OffsetIterator linestring2_offsets_first, - Cart2dItB linestring2_points_first, - Cart2dItB linestring2_points_last, +template +OutputIt pairwise_linestring_distance(MultiLinestringRange1 multilinestrings1, + MultiLinestringRange2 multilinestrings2, OutputIt distances_first, rmm::cuda_stream_view stream) { - using T = typename cuspatial::iterator_vec_base_type; + using T = typename MultiLinestringRange1::element_t; - static_assert(is_same_floating_point, - typename cuspatial::iterator_value_type>(), + static_assert(is_same_floating_point(), "Inputs and output must have the same floating point value type."); static_assert(is_same, - typename cuspatial::iterator_value_type, - typename cuspatial::iterator_value_type>(), + typename MultiLinestringRange1::point_t, + typename MultiLinestringRange2::point_t>(), "All input types must be cuspatial::vec_2d with the same value type"); - auto const num_linestring_pairs = - thrust::distance(linestring1_offsets_first, linestring1_offsets_last) - 1; - auto const num_linestring1_points = - thrust::distance(linestring1_points_first, linestring1_points_last); - auto const num_linestring2_points = - thrust::distance(linestring2_points_first, linestring2_points_last); + CUSPATIAL_EXPECTS(multilinestrings1.size() == multilinestrings2.size(), + "Inputs must have the same number of rows."); thrust::fill(rmm::exec_policy(stream), distances_first, - distances_first + num_linestring_pairs, + distances_first + multilinestrings1.size(), std::numeric_limits::max()); - std::size_t constexpr threads_per_block = 64; + std::size_t constexpr threads_per_block = 256; std::size_t const num_blocks = - (num_linestring1_points + threads_per_block - 1) / threads_per_block; + (multilinestrings1.num_points() + threads_per_block - 1) / threads_per_block; detail::pairwise_linestring_distance_kernel<<>>( - linestring1_offsets_first, - linestring1_offsets_last, - linestring1_points_first, - linestring1_points_last, - linestring2_offsets_first, - linestring2_points_first, - linestring2_points_last, - distances_first); + multilinestrings1, multilinestrings2, distances_first); CUSPATIAL_CUDA_TRY(cudaGetLastError()); - return distances_first + num_linestring_pairs; + return distances_first + multilinestrings1.size(); } } // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.cuh index e6b664699..3188491d7 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.cuh @@ -23,7 +23,7 @@ #include #include -#include +#include #include #include @@ -31,7 +31,36 @@ namespace cuspatial { -using namespace cuspatial::detail; +using namespace detail; + +template +struct to_multilinestring_functor { + using difference_type = typename thrust::iterator_difference::type; + GeometryIterator _geometry_begin; + PartIterator _part_begin; + VecIterator _point_begin; + VecIterator _point_end; + + CUSPATIAL_HOST_DEVICE + to_multilinestring_functor(GeometryIterator geometry_begin, + PartIterator part_begin, + VecIterator point_begin, + VecIterator point_end) + : _geometry_begin(geometry_begin), + _part_begin(part_begin), + _point_begin(point_begin), + _point_end(point_end) + { + } + + CUSPATIAL_HOST_DEVICE auto operator()(difference_type i) + { + return multilinestring_ref{_part_begin + _geometry_begin[i], + thrust::next(_part_begin + _geometry_begin[i + 1]), + _point_begin, + _point_end}; + } +}; template class multilinestring_range; @@ -42,14 +71,14 @@ multilinestring_range::multilinestr GeometryIterator geometry_end, PartIterator part_begin, PartIterator part_end, - VecIterator points_begin, - VecIterator points_end) - : geometry_begin(geometry_begin), - geometry_end(geometry_end), - part_begin(part_begin), - part_end(part_end), - points_begin(points_begin), - points_end(points_end) + VecIterator point_begin, + VecIterator point_end) + : _geometry_begin(geometry_begin), + _geometry_end(geometry_end), + _part_begin(part_begin), + _part_end(part_end), + _point_begin(point_begin), + _point_end(point_end) { } @@ -57,21 +86,36 @@ template ::num_multilinestrings() { - return thrust::distance(geometry_begin, geometry_end) - 1; + return thrust::distance(_geometry_begin, _geometry_end) - 1; } template CUSPATIAL_HOST_DEVICE auto multilinestring_range::num_linestrings() { - return thrust::distance(part_begin, part_end) - 1; + return thrust::distance(_part_begin, _part_end) - 1; } template CUSPATIAL_HOST_DEVICE auto multilinestring_range::num_points() { - return thrust::distance(points_begin, points_end); + return thrust::distance(_point_begin, _point_end); +} + +template +CUSPATIAL_HOST_DEVICE auto +multilinestring_range::multilinestring_begin() +{ + return detail::make_counting_transform_iterator( + 0, to_multilinestring_functor{_geometry_begin, _part_begin, _point_begin, _point_end}); +} + +template +CUSPATIAL_HOST_DEVICE auto +multilinestring_range::multilinestring_end() +{ + return multilinestring_begin() + num_multilinestrings(); } template @@ -80,8 +124,8 @@ CUSPATIAL_HOST_DEVICE auto multilinestring_range::part_idx_from_point_idx( IndexType point_idx) { - auto part_it = thrust::upper_bound(thrust::seq, part_begin, part_end, point_idx); - return thrust::distance(part_begin, thrust::prev(part_it)); + auto part_it = thrust::upper_bound(thrust::seq, _part_begin, _part_end, point_idx); + return thrust::distance(_part_begin, thrust::prev(part_it)); } template @@ -90,8 +134,8 @@ CUSPATIAL_HOST_DEVICE auto multilinestring_range::geometry_idx_from_part_idx( IndexType part_idx) { - auto geom_it = thrust::upper_bound(thrust::seq, geometry_begin, geometry_end, part_idx); - return thrust::distance(geometry_begin, thrust::prev(geom_it)); + auto geom_it = thrust::upper_bound(thrust::seq, _geometry_begin, _geometry_end, part_idx); + return thrust::distance(_geometry_begin, thrust::prev(geom_it)); } template @@ -110,9 +154,9 @@ multilinestring_range::is_valid_seg IndexType1 segment_idx, IndexType2 part_idx) { if constexpr (std::is_signed_v) - return segment_idx >= 0 && segment_idx < (part_begin[part_idx + 1] - 1); + return segment_idx >= 0 && segment_idx < (_part_begin[part_idx + 1] - 1); else - return segment_idx < (part_begin[part_idx + 1] - 1); + return segment_idx < (_part_begin[part_idx + 1] - 1); } template @@ -122,7 +166,16 @@ CUSPATIAL_HOST_DEVICE thrust::pair< vec_2d::element_t>> multilinestring_range::segment(IndexType segment_idx) { - return thrust::make_pair(points_begin[segment_idx], points_begin[segment_idx + 1]); + return thrust::make_pair(_point_begin[segment_idx], _point_begin[segment_idx + 1]); +} + +template +template +CUSPATIAL_HOST_DEVICE auto +multilinestring_range::operator[]( + IndexType multilinestring_idx) +{ + return multilinestring_begin()[multilinestring_idx]; } } // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/geometry/linestring_ref.cuh b/cpp/include/cuspatial/experimental/geometry/linestring_ref.cuh new file mode 100644 index 000000000..b6135e433 --- /dev/null +++ b/cpp/include/cuspatial/experimental/geometry/linestring_ref.cuh @@ -0,0 +1,49 @@ +/* + * 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 + +namespace cuspatial { + +/** + * @brief Represent a reference to a linestring stored in a structure of arrays. + * + * @tparam VecIterator type of iterator to the underlying point array. + */ +template +class linestring_ref { + public: + CUSPATIAL_HOST_DEVICE linestring_ref(VecIterator begin, VecIterator end); + + CUSPATIAL_HOST_DEVICE auto num_segments() const; + + /// Return iterator to the first segment of the linestring + CUSPATIAL_HOST_DEVICE auto segment_begin() const; + /// Return iterator to one past the last segment + CUSPATIAL_HOST_DEVICE auto segment_end() const; + + /// Return iterator to the first segment of the linestring + CUSPATIAL_HOST_DEVICE auto begin() const { return segment_begin(); } + /// Return iterator to one past the last segment + CUSPATIAL_HOST_DEVICE auto end() const { return segment_end(); } + + protected: + VecIterator _point_begin; + VecIterator _point_end; +}; + +} // namespace cuspatial +#include diff --git a/cpp/include/cuspatial/experimental/geometry_collection/multilinestring_ref.cuh b/cpp/include/cuspatial/experimental/geometry_collection/multilinestring_ref.cuh new file mode 100644 index 000000000..d48830668 --- /dev/null +++ b/cpp/include/cuspatial/experimental/geometry_collection/multilinestring_ref.cuh @@ -0,0 +1,63 @@ +/* + * 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 + +namespace cuspatial { + +/** + * @brief Represent a reference to a multilinestring stored in a structure of arrays. + * + * @tparam PartIterator type of iterator to the part offset array. + * @tparam VecIterator type of iterator to the underlying point array. + */ +template +class multilinestring_ref { + public: + CUSPATIAL_HOST_DEVICE multilinestring_ref(PartIterator part_begin, + PartIterator part_end, + VecIterator point_begin, + VecIterator point_end); + /// Return the number of linestrings in the multilinestring. + CUSPATIAL_HOST_DEVICE auto num_linestrings() const; + /// Return the number of linestrings in the multilinestring. + CUSPATIAL_HOST_DEVICE auto size() const { return num_linestrings(); } + + /// Return iterator to the first linestring. + CUSPATIAL_HOST_DEVICE auto part_begin() const; + /// Return iterator to one past the last linestring. + CUSPATIAL_HOST_DEVICE auto part_end() const; + + /// Return iterator to the first point of the multilinestring. + CUSPATIAL_HOST_DEVICE auto point_begin() const; + /// Return iterator to one past the last point of the multilinestring. + CUSPATIAL_HOST_DEVICE auto point_end() const; + + /// Return iterator to the first linestring of the multilinestring. + CUSPATIAL_HOST_DEVICE auto begin() const { return part_begin(); } + /// Return iterator to one past the last linestring of the multilinestring. + CUSPATIAL_HOST_DEVICE auto end() const { return part_end(); } + + protected: + PartIterator _part_begin; + PartIterator _part_end; + VecIterator _point_begin; + VecIterator _point_end; +}; + +} // namespace cuspatial + +#include diff --git a/cpp/include/cuspatial/experimental/geometry_collection/multipoint_ref.cuh b/cpp/include/cuspatial/experimental/geometry_collection/multipoint_ref.cuh index 880b80499..4acb55a63 100644 --- a/cpp/include/cuspatial/experimental/geometry_collection/multipoint_ref.cuh +++ b/cpp/include/cuspatial/experimental/geometry_collection/multipoint_ref.cuh @@ -19,7 +19,7 @@ namespace cuspatial { /** - * @brief Represent a multipoint stored in structure of array on memory. + * @brief Represent a reference to multipoint stored in a structure of arrays. * * @tparam VecIterator type of iterator to the underlying point array. */ diff --git a/cpp/include/cuspatial/experimental/linestring_distance.cuh b/cpp/include/cuspatial/experimental/linestring_distance.cuh index 684e5696e..31c7eaa6b 100644 --- a/cpp/include/cuspatial/experimental/linestring_distance.cuh +++ b/cpp/include/cuspatial/experimental/linestring_distance.cuh @@ -28,48 +28,23 @@ namespace cuspatial { * between all pairs of segments of the two linestrings. If any of the segments intersect, * the distance is 0. * - * @tparam Cart2dItA iterator type for point array of the first linestring of each pair. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam Cart2dItB iterator type for point array of the second linestring of each pair. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam OffsetIterator iterator type for offset array. Must meet the requirements of - * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam OutputIt iterator type for output array. Must meet the requirements of - * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. + * @tparam MultiLinestringRange an instance of template type `multilinestring_range` + * @tparam OutputIt iterator type for output array. Must meet the requirements of [LRAI](LinkLRAI) + * and be device writable. * - * @param linestring1_offsets_first beginning of range of the offsets to the first linestring of - * each pair - * @param linestring1_offsets_last end of range of the offsets to the first linestring of each pair - * @param linestring1_points_first beginning of range of the point of the first linestring of each - * pair - * @param linestring1_points_last end of range of the point of the first linestring of each pair - * @param linestring2_offsets_first beginning of range of the offsets to the second linestring of - * each pair - * @param linestring2_points_first beginning of range of the point of the second linestring of each - * pair - * @param linestring2_points_last end of range of the point of the second linestring of each pair - * @param distances_first beginning iterator to output + * @param multilinestrings1 Range object of the lhs multilinestring array + * @param multilinestrings2 Range object of the rhs multilinestring array * @param stream The CUDA stream to use for device memory operations and kernel launches. - * @return Output iterator to one past the last element in the output range - * - * @pre all input iterators for coordinates must have `cuspatial::vec_2d` type. - * @pre all scalar types must be floating point types, and must be the same type for all input - * iterators and output iterators. + * @return Output iterator to the element past the last distance computed. * * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator * "LegacyRandomAccessIterator" */ -template -OutputIt pairwise_linestring_distance(OffsetIterator linestring1_offsets_first, - OffsetIterator linestring1_offsets_last, - Cart2dItA linestring1_points_first, - Cart2dItA linestring1_points_last, - OffsetIterator linestring2_offsets_first, - Cart2dItB linestring2_points_first, - Cart2dItB linestring2_points_last, +template +OutputIt pairwise_linestring_distance(MultiLinestringRange1 multilinestrings1, + MultiLinstringRange2 multilinestrings2, OutputIt distances_first, rmm::cuda_stream_view stream = rmm::cuda_stream_default); - } // namespace cuspatial #include diff --git a/cpp/include/cuspatial/experimental/ranges/multilinestring_range.cuh b/cpp/include/cuspatial/experimental/ranges/multilinestring_range.cuh index e1560c5e9..09d96ca71 100644 --- a/cpp/include/cuspatial/experimental/ranges/multilinestring_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multilinestring_range.cuh @@ -63,68 +63,66 @@ class multilinestring_range { VecIterator points_begin, VecIterator points_end); - /** - * @brief Return the number of multilinestrings in the array. - */ + /// Return the number of multilinestrings in the array. CUSPATIAL_HOST_DEVICE auto size() { return num_multilinestrings(); } - /** - * @brief Return the number of multilinestrings in the array. - */ + /// Return the number of multilinestrings in the array. CUSPATIAL_HOST_DEVICE auto num_multilinestrings(); - /** - * @brief Return the total number of linestrings in the array. - */ + /// Return the total number of linestrings in the array. CUSPATIAL_HOST_DEVICE auto num_linestrings(); - /** - * @brief Return the total number of points in the array. - */ + /// Return the total number of points in the array. CUSPATIAL_HOST_DEVICE auto num_points(); - /** - * @brief Given the index of a point, return the part (linestring) index where the point locates. - */ + /// Return the iterator to the first multilinestring in the range. + CUSPATIAL_HOST_DEVICE auto multilinestring_begin(); + + /// Return the iterator to the one past the last multilinestring in the range. + CUSPATIAL_HOST_DEVICE auto multilinestring_end(); + + /// Return the iterator to the first multilinestring in the range. + CUSPATIAL_HOST_DEVICE auto begin() { return multilinestring_begin(); } + + /// Return the iterator to the one past the last multilinestring in the range. + CUSPATIAL_HOST_DEVICE auto end() { return multilinestring_end(); } + + /// Given the index of a point, return the part (linestring) index where the point locates. template CUSPATIAL_HOST_DEVICE auto part_idx_from_point_idx(IndexType point_idx); - /** - * @brief Given the index of a part (linestring), return the geometry (multilinestring) index - * where the linestring locates. - */ + /// Given the index of a part (linestring), return the geometry (multilinestring) index + /// where the linestring locates. template CUSPATIAL_HOST_DEVICE auto geometry_idx_from_part_idx(IndexType part_idx); - /** - * @brief Given the index of a point, return the geometry (multilinestring) index where the - * point locates. - */ + /// Given the index of a point, return the geometry (multilinestring) index where the + /// point locates. template CUSPATIAL_HOST_DEVICE auto geometry_idx_from_point_idx(IndexType point_idx); - /** - * @brief Given an index of a segment, returns true if the index is valid. - * The index of a segment is the same as the index to the starting point of the segment. - * Thus, the index to the last point of a linestring is an invalid segment index. - */ + /// Given an index of a segment, returns true if the index is valid. + /// The index of a segment is the same as the index to the starting point of the segment. + /// Thus, the index to the last point of a linestring is an invalid segment index. template CUSPATIAL_HOST_DEVICE bool is_valid_segment_id(IndexType1 segment_idx, IndexType2 part_idx); - /** - * @brief Returns the segment given a segment index. - */ + /// Returns the segment given a segment index. template CUSPATIAL_HOST_DEVICE thrust::pair, vec_2d> segment( IndexType segment_idx); + /// Returns the `multilinestring_idx`th multilinestring in the range. + template + CUSPATIAL_HOST_DEVICE auto operator[](IndexType multilinestring_idx); + protected: - GeometryIterator geometry_begin; - GeometryIterator geometry_end; - PartIterator part_begin; - PartIterator part_end; - VecIterator points_begin; - VecIterator points_end; + GeometryIterator _geometry_begin; + GeometryIterator _geometry_end; + PartIterator _part_begin; + PartIterator _part_end; + VecIterator _point_begin; + VecIterator _point_end; }; /** diff --git a/cpp/include/cuspatial_test/vector_factories.cuh b/cpp/include/cuspatial_test/vector_factories.cuh new file mode 100644 index 000000000..8cacef3df --- /dev/null +++ b/cpp/include/cuspatial_test/vector_factories.cuh @@ -0,0 +1,25 @@ +/* + * 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 + +#include + +template +auto make_device_vector(std::initializer_list inl) +{ + return rmm::device_vector(inl.begin(), inl.end()); +} diff --git a/cpp/src/spatial/linestring_distance.cu b/cpp/src/spatial/linestring_distance.cu index 44baae1f0..a3d503e77 100644 --- a/cpp/src/spatial/linestring_distance.cu +++ b/cpp/src/spatial/linestring_distance.cu @@ -17,24 +17,17 @@ #include #include #include +#include #include #include #include #include -#include -#include #include #include -#include -#include -#include -#include -#include #include -#include #include #include @@ -67,15 +60,22 @@ struct pairwise_linestring_distance_functor { auto linestring2_coords_it = make_vec_2d_iterator(linestring2_points_x.begin(), linestring2_points_y.begin()); - pairwise_linestring_distance(linestring1_offsets.begin(), - linestring1_offsets.end(), - linestring1_coords_it, - linestring1_coords_it + linestring1_points_x.size(), - linestring2_offsets.begin(), - linestring2_coords_it, - linestring2_coords_it + linestring2_points_x.size(), - distances->mutable_view().begin(), - stream); + 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); + + pairwise_linestring_distance( + multilinestrings1, multilinestrings2, distances->mutable_view().begin()); return distances; } diff --git a/cpp/tests/experimental/spatial/linestring_distance_test.cu b/cpp/tests/experimental/spatial/linestring_distance_test.cu index 4bc59ee82..b4ccb77e5 100644 --- a/cpp/tests/experimental/spatial/linestring_distance_test.cu +++ b/cpp/tests/experimental/spatial/linestring_distance_test.cu @@ -14,28 +14,28 @@ * limitations under the License. */ -#include "tests/utility/vector_equality.hpp" +#include + +#include #include #include #include +#include #include - -#include -#include - #include #include namespace cuspatial { namespace test { -using namespace cudf; - template struct PairwiseLinestringDistanceTest : public ::testing::Test { }; +struct PairwiseLinestringDistanceTestUntyped : public ::testing::Test { +}; + // float and double are logically the same but would require seperate tests due to precision. using TestTypes = ::testing::Types; TYPED_TEST_CASE(PairwiseLinestringDistanceTest, TestTypes); @@ -45,26 +45,34 @@ TYPED_TEST(PairwiseLinestringDistanceTest, FromSeparateArrayInputs) using T = TypeParam; using CartVec = std::vector>; + auto constexpr num_pairs = 1; + auto a_cart2d = rmm::device_vector>{ CartVec({{0.0f, 0.0f}, {1.0f, 0.0f}, {2.0f, 0.0f}, {3.0f, 0.0f}, {4.0f, 0.0f}})}; auto b_cart2d = rmm::device_vector>{ CartVec({{0.0f, 1.0f}, {1.0f, 1.0f}, {2.0f, 1.0f}, {3.0f, 1.0f}, {4.0f, 1.0f}})}; - auto offset = rmm::device_vector{std::vector{0, 5}}; - - auto distance = rmm::device_vector{1}; + auto offset = make_device_vector({0, 5}); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + offset.size() - 1, + offset.begin(), + a_cart2d.size(), + a_cart2d.begin()); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + offset.size() - 1, + offset.begin(), + b_cart2d.size(), + b_cart2d.begin()); + + auto got = rmm::device_vector(num_pairs); auto expected = rmm::device_vector{std::vector{1.0}}; - auto ret = pairwise_linestring_distance(offset.begin(), - offset.end(), - a_cart2d.begin(), - a_cart2d.end(), - offset.begin(), - b_cart2d.begin(), - b_cart2d.end(), - distance.begin()); + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); - test::expect_vector_equivalent(expected, distance); - EXPECT_EQ(offset.size() - 1, std::distance(distance.begin(), ret)); + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); } TYPED_TEST(PairwiseLinestringDistanceTest, FromSamePointArrayInput) @@ -72,30 +80,34 @@ TYPED_TEST(PairwiseLinestringDistanceTest, FromSamePointArrayInput) using T = TypeParam; using CartVec = std::vector>; - auto cart2ds = rmm::device_vector>{ - CartVec({{0.0f, 0.0f}, {1.0f, 0.0f}, {2.0f, 0.0f}, {3.0f, 0.0f}, {4.0f, 0.0f}})}; - auto offset_a = rmm::device_vector{std::vector{0, 3}}; - auto offset_b = rmm::device_vector{std::vector{0, 4}}; + auto constexpr num_pairs = 1; - auto a_begin = cart2ds.begin(); - auto a_end = cart2ds.begin() + 3; - auto b_begin = cart2ds.begin() + 1; - auto b_end = cart2ds.end(); + auto cart2ds = make_device_vector>( + {{0.0f, 0.0f}, {1.0f, 0.0f}, {2.0f, 0.0f}, {3.0f, 0.0f}, {4.0f, 0.0f}}); + auto offset_a = make_device_vector({0, 3}); + auto offset_b = make_device_vector({0, 4}); - auto distance = rmm::device_vector{1}; - auto expected = rmm::device_vector{std::vector{0.0}}; + auto got = rmm::device_vector(1); + auto expected = make_device_vector({0.0}); - auto ret = pairwise_linestring_distance(offset_a.begin(), - offset_a.end(), - a_begin, - a_end, - offset_a.begin(), - b_begin, - b_end, - distance.begin()); + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + offset_a.size() - 1, + offset_a.begin(), + 3, + cart2ds.begin()); - test::expect_vector_equivalent(expected, distance); - EXPECT_EQ(offset_a.size() - 1, std::distance(distance.begin(), ret)); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + offset_b.size() - 1, + offset_b.begin(), + 4, + thrust::next(cart2ds.begin())); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); } TYPED_TEST(PairwiseLinestringDistanceTest, FromTransformIterator) @@ -103,28 +115,41 @@ TYPED_TEST(PairwiseLinestringDistanceTest, FromTransformIterator) using T = TypeParam; using CartVec = std::vector>; + auto constexpr num_pairs = 1; + auto a_cart2d_x = rmm::device_vector{std::vector{0.0, 1.0, 2.0, 3.0, 4.0}}; auto a_cart2d_y = rmm::device_vector(5, 0.0); auto a_begin = make_vec_2d_iterator(a_cart2d_x.begin(), a_cart2d_y.begin()); - auto a_end = a_begin + a_cart2d_x.size(); auto b_cart2d_x = rmm::device_vector{std::vector{0.0, 1.0, 2.0, 3.0, 4.0}}; auto b_cart2d_y = rmm::device_vector(5, 1.0); auto b_begin = make_vec_2d_iterator(b_cart2d_x.begin(), b_cart2d_y.begin()); - auto b_end = b_begin + b_cart2d_x.size(); auto offset = rmm::device_vector{std::vector{0, 5}}; - auto distance = rmm::device_vector{1}; + auto got = rmm::device_vector{1}; auto expected = rmm::device_vector{std::vector{1.0}}; - auto ret = pairwise_linestring_distance( - offset.begin(), offset.end(), a_begin, a_end, offset.begin(), b_begin, b_end, distance.begin()); + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + offset.size() - 1, + offset.begin(), + a_cart2d_x.size(), + a_begin); + + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + offset.size() - 1, + offset.begin(), + b_cart2d_x.size(), + b_begin); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); - test::expect_vector_equivalent(expected, distance); - EXPECT_EQ(offset.size() - 1, std::distance(distance.begin(), ret)); + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); } TYPED_TEST(PairwiseLinestringDistanceTest, FromMixedIterator) @@ -132,6 +157,8 @@ TYPED_TEST(PairwiseLinestringDistanceTest, FromMixedIterator) using T = TypeParam; using CartVec = std::vector>; + auto constexpr num_pairs = 1; + auto a_cart2d = rmm::device_vector>{ CartVec({{0.0f, 0.0f}, {1.0f, 0.0f}, {2.0f, 0.0f}, {3.0f, 0.0f}, {4.0f, 0.0f}})}; @@ -139,25 +166,31 @@ TYPED_TEST(PairwiseLinestringDistanceTest, FromMixedIterator) auto b_cart2d_y = rmm::device_vector(5, 1.0); auto b_begin = make_vec_2d_iterator(b_cart2d_x.begin(), b_cart2d_y.begin()); - auto b_end = b_begin + b_cart2d_x.size(); auto offset_a = rmm::device_vector{std::vector{0, 5}}; auto offset_b = rmm::device_vector{std::vector{0, 5}}; - auto distance = rmm::device_vector{1}; + auto got = rmm::device_vector{1}; auto expected = rmm::device_vector{std::vector{1.0}}; - auto ret = pairwise_linestring_distance(offset_a.begin(), - offset_a.end(), - a_cart2d.begin(), - a_cart2d.end(), - offset_b.begin(), - b_begin, - b_end, - distance.begin()); + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + offset_a.size() - 1, + offset_a.begin(), + a_cart2d.size(), + a_cart2d.begin()); + + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + offset_b.size() - 1, + offset_b.begin(), + b_cart2d_x.size(), + b_begin); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); - test::expect_vector_equivalent(expected, distance); - EXPECT_EQ(offset_a.size() - 1, std::distance(distance.begin(), ret)); + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); } TYPED_TEST(PairwiseLinestringDistanceTest, FromLongInputs) @@ -165,35 +198,756 @@ TYPED_TEST(PairwiseLinestringDistanceTest, FromLongInputs) using T = TypeParam; using CartVec = std::vector>; - auto num_points = 1000; + auto constexpr num_pairs = 5; + auto constexpr num_points = 1000; auto a_cart2d_x_begin = thrust::make_constant_iterator(T{0.0}); auto a_cart2d_y_begin = thrust::make_counting_iterator(T{0.0}); auto a_cart2d_begin = make_vec_2d_iterator(a_cart2d_x_begin, a_cart2d_y_begin); - auto a_cart2d_end = a_cart2d_begin + num_points; auto b_cart2d_x_begin = thrust::make_constant_iterator(T{42.0}); auto b_cart2d_y_begin = thrust::make_counting_iterator(T{0.0}); auto b_cart2d_begin = make_vec_2d_iterator(b_cart2d_x_begin, b_cart2d_y_begin); - auto b_cart2d_end = b_cart2d_begin + num_points; auto offset = rmm::device_vector{std::vector{0, 100, 200, 300, 400, num_points}}; - auto distance = rmm::device_vector{5}; + auto got = rmm::device_vector{num_pairs}; auto expected = rmm::device_vector{std::vector{42.0, 42.0, 42.0, 42.0, 42.0}}; - auto ret = pairwise_linestring_distance(offset.begin(), - offset.end(), - a_cart2d_begin, - a_cart2d_end, - offset.begin(), - b_cart2d_begin, - b_cart2d_end, - distance.begin()); - - test::expect_vector_equivalent(expected, distance); - EXPECT_EQ(offset.size() - 1, std::distance(distance.begin(), ret)); + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + offset.size() - 1, + offset.begin(), + num_points, + a_cart2d_begin); + + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + offset.size() - 1, + offset.begin(), + num_points, + b_cart2d_begin); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairLinestringParallel) +{ + using T = TypeParam; + // Linestring 1: (0.0, 0.0), (1.0, 1.0) + // Linestring 2: (1.0, 0.0), (2.0, 1.0) + + int32_t constexpr num_pairs = 1; + + auto linestring1_offsets = make_device_vector({0, 2}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 1.0, 1.0}); + auto linestring2_offsets = make_device_vector({0, 2}); + auto linestring2_points_xy = make_device_vector({1.0, 0.0, 2.0, 1.0}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({0.7071067811865476}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairLinestringEndpointsDistance) +{ + using T = TypeParam; + // Linestring 1: (0.0, 0.0), (1.0, 1.0), (2.0, 2.0) + // Linestring 2: (2.0, 0.0), (1.0, -1.0), (0.0, -1.0) + + auto constexpr num_pairs = 1; + + auto linestring1_offsets = make_device_vector({0, 3}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 1.0, 1.0, 2.0, 2.0}); + auto linestring2_offsets = make_device_vector({0, 3}); + auto linestring2_points_xy = make_device_vector({2.0, 0.0, 1.0, -1.0, 0.0, -1.0}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({1.0}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairLinestringProjectionNotOnLine) +{ + using T = TypeParam; + + auto constexpr num_pairs = 1; + + // Linestring 1: (0.0, 0.0), (1.0, 1.0) + // Linestring 2: (3.0, 1.5), (3.0, 2.0) + auto linestring1_offsets = make_device_vector({0, 2}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 1.0, 1.0}); + auto linestring2_offsets = make_device_vector({0, 2}); + auto linestring2_points_xy = make_device_vector({3.0, 1.5, 3.0, 2.0}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({2.0615528128088303}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairLinestringPerpendicular) +{ + using T = TypeParam; + + auto constexpr num_pairs = 1; + + // Linestring 1: (0.0, 0.0), (2.0, 0.0) + // Linestring 2: (1.0, 1.0), (1.0, 2.0) + auto linestring1_offsets = make_device_vector({0, 2}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 2.0, 0.0}); + auto linestring2_offsets = make_device_vector({0, 2}); + auto linestring2_points_xy = make_device_vector({1.0, 1.0, 1.0, 2.0}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({1.0}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairLinestringIntersects) +{ + using T = TypeParam; + + auto constexpr num_pairs = 1; + + // Linestring 1: (0.0, 0.0), (1.0, 1.0) + // Linestring 2: (0.0, 1.0), (1.0, 0.0) + auto linestring1_offsets = make_device_vector({0, 2}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 1.0, 1.0}); + auto linestring2_offsets = make_device_vector({0, 2}); + auto linestring2_points_xy = make_device_vector({0.0, 1.0, 1.0, 0.0}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({0.0}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairLinestringSharedVertex) +{ + using T = TypeParam; + + auto constexpr num_pairs = 1; + + // Linestring 1: (0.0, 0.0), (0.0, 2.0), (2.0, 2.0) + // Linestring 2: (2.0, 2.0), (2.0, 1.0), (1.0, 1.0), (2.5, 0.0) + auto linestring1_offsets = make_device_vector({0, 3}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 0.0, 2.0, 2.0, 2.0}); + auto linestring2_offsets = make_device_vector({0, 4}); + auto linestring2_points_xy = make_device_vector({2.0, 2.0, 2.0, 1.0, 1.0, 1.0, 2.5, 0.0}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({0.0}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairLinestringCoincide) +{ + using T = TypeParam; + + auto constexpr num_pairs = 1; + // Linestring 1: (0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0) + // Linestring 2: (2.0, 1.0), (1.0, 1.0), (1.0, 0.0), (2.0, 0.0), (2.0, 0.5) + auto linestring1_offsets = make_device_vector({0, 4}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0}); + auto linestring2_offsets = make_device_vector({0, 5}); + auto linestring2_points_xy = + make_device_vector({2.0, 1.0, 1.0, 1.0, 1.0, 0.0, 2.0, 0.0, 2.0, 0.5}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({0.0}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairDegenerateCollinearNoIntersect) +{ + using T = TypeParam; + + auto constexpr num_pairs = 1; + + // Linestring1: (0.0, 0.0) -> (0.0, 1.0) + // Linestring2: (0.0, 2.0) -> (0.0, 3.0) + auto linestring1_offsets = make_device_vector({0, 2}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 0.0, 1.0}); + auto linestring2_offsets = make_device_vector({0, 2}); + auto linestring2_points_xy = make_device_vector({0.0, 2.0, 0.0, 3.0}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({1.0}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairCollinearNoIntersect) +{ + using T = TypeParam; + + auto constexpr num_pairs = 1; + + // Linestring1: (0.0, 0.0) -> (1.0, 1.0) + // Linestring2: (2.0, 2.0) -> (3.0, 3.0) + auto linestring1_offsets = make_device_vector({0, 2}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 1.0, 1.0}); + auto linestring2_offsets = make_device_vector({0, 2}); + auto linestring2_points_xy = make_device_vector({2.0, 2.0, 3.0, 3.0}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({1.4142135623730951}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairDegenerateCollinearIntersect) +{ + using T = TypeParam; + + auto constexpr num_pairs = 1; + + auto linestring1_offsets = make_device_vector({0, 2}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 2.0, 2.0}); + auto linestring2_offsets = make_device_vector({0, 2}); + auto linestring2_points_xy = make_device_vector({1.0, 1.0, 3.0, 3.0}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({0.0}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TEST_F(PairwiseLinestringDistanceTestUntyped, OnePairDeterminantDoublePrecisionDenormalized) +{ + // Vector ab: (1e-155, 2e-155) + // Vector cd: (2e-155, 1e-155) + // determinant of matrix [a, b] = -3e-310, a denormalized number + + auto constexpr num_pairs = 1; + + auto linestring1_offsets = make_device_vector({0, 2}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 1e-155, 2e-155}); + auto linestring2_offsets = make_device_vector({0, 2}); + auto linestring2_points_xy = make_device_vector({4e-155, 5e-155, 6e-155, 6e-155}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({4.24264068711929e-155}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TEST_F(PairwiseLinestringDistanceTestUntyped, OnePairDeterminantSinglePrecisionDenormalized) +{ + // Vector ab: (1e-20, 2e-20) + // Vector cd: (2e-20, 1e-20) + // determinant of matrix [ab, cd] = -3e-40, a denormalized number + + auto constexpr num_pairs = 1; + + auto linestring1_offsets = make_device_vector({0, 2}); + auto linestring1_points_xy = make_device_vector({0.0, 0.0, 1e-20, 2e-20}); + auto linestring2_offsets = make_device_vector({0, 2}); + auto linestring2_points_xy = make_device_vector({4e-20, 5e-20, 6e-20, 6e-20}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({4.2426405524813e-20}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + // Expect a slightly greater floating point error compared to 4 ULP + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got, 1e-25f); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairRandom1) +{ + using T = TypeParam; + + auto constexpr num_pairs = 1; + auto linestring1_offsets = make_device_vector({0, 3}); + auto linestring1_points_xy = make_device_vector({-22556.235212018168, + 41094.0501840996, + -16375.655690574613, + 42992.319790050366, + -20082.724633593425, + 33759.13529113619}); + auto linestring2_offsets = make_device_vector({0, 2}); + auto linestring2_points_xy = make_device_vector( + {4365.496374409238, -59857.47177852941, 1671.0269165650761, -54931.9723439855}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({91319.97744223749}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairIntersectFromRealData1) +{ + // Example extracted from a pair of trajectry in geolife dataset + using T = TypeParam; + auto constexpr num_pairs = 1; + + auto linestring1_offsets = make_device_vector({0, 5}); + auto linestring1_points_xy = make_device_vector({39.97551667, + 116.33028333, + 39.97585, + 116.3304, + 39.97598333, + 116.33046667, + 39.9761, + 116.3305, + 39.97623333, + 116.33056667}); + + auto linestring2_offsets = make_device_vector({0, 55}); + auto linestring2_points_xy = make_device_vector( + {39.97381667, 116.34211667, 39.97341667, 116.34215, 39.9731, 116.34218333, + 39.97293333, 116.34221667, 39.97233333, 116.34225, 39.97218333, 116.34243333, + 39.97218333, 116.34296667, 39.97215, 116.34478333, 39.97168333, 116.34486667, + 39.97093333, 116.34485, 39.97073333, 116.34468333, 39.9705, 116.34461667, + 39.96991667, 116.34465, 39.96961667, 116.34465, 39.96918333, 116.34466667, + 39.96891667, 116.34465, 39.97531667, 116.33036667, 39.97533333, 116.32961667, + 39.97535, 116.3292, 39.97515, 116.32903333, 39.97506667, 116.32985, + 39.97508333, 116.33128333, 39.9751, 116.33195, 39.97513333, 116.33618333, + 39.97511667, 116.33668333, 39.97503333, 116.33818333, 39.97513333, 116.34, + 39.97523333, 116.34045, 39.97521667, 116.34183333, 39.97503333, 116.342, + 39.97463333, 116.34203333, 39.97443333, 116.3422, 39.96838333, 116.3445, + 39.96808333, 116.34451667, 39.96771667, 116.3445, 39.96745, 116.34453333, + 39.96735, 116.34493333, 39.9673, 116.34506667, 39.96718333, 116.3451, + 39.96751667, 116.34483333, 39.9678, 116.3448, 39.9676, 116.3449, + 39.96741667, 116.345, 39.9672, 116.34506667, 39.97646667, 116.33006667, + 39.9764, 116.33015, 39.97625, 116.33026667, 39.9762, 116.33038333, + 39.97603333, 116.33036667, 39.97581667, 116.3303, 39.9757, 116.33033333, + 39.97551667, 116.33035, 39.97535, 116.3304, 39.97543333, 116.33078333, + 39.97538333, 116.33066667}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({0.0}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, OnePairFromRealData) +{ + // Example extracted from a pair of trajectry in geolife dataset + using T = TypeParam; + + auto constexpr num_pairs = 1; + auto linestring1_offsets = make_device_vector({0, 2}); + auto linestring1_points_xy = + make_device_vector({39.9752666666667, 116.334316666667, 39.9752666666667, 116.334533333333}); + auto linestring2_offsets = make_device_vector({0, 2}); + auto linestring2_points_xy = + make_device_vector({39.9752666666667, 116.323966666667, 39.9752666666667, 116.3236}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + // Steps to reproduce: + // Create a float32/float64 numpy array with the literal inputs as above. + // Construct a shapely.geometry.LineString object and compute the result. + // Cast the result to np.float32/np.float64. + auto expected = + make_device_vector({std::is_same_v ? 0.010353088f : 0.010349999999988313}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, TwoPairsSingleLineString) +{ + using T = TypeParam; + auto constexpr num_pairs = 2; + + // First pair lhs: + // (41658.902315589876, 14694.11814724456)->(46600.70359801489, 8771.431887804214)-> + // (47079.510547637154, 10199.68027155776)->(51498.48049880379, 17049.62665643919) + // First pair rhs: + // (24046.170375947084, 27878.56737867571)->(20614.007047185743, 26489.74880629428) + + // Second pair lhs: + // (-27429.917796286478, -33240.8339287343)->(-21764.269974046114, -37974.45515744517)-> + // (-14460.71813363161, -31333.481529957502)->(-18226.13032712476, -30181.03842467982) + // Second pair rhs: + // (48381.39607717942, -8366.313156569413)->(53346.77764665915, -2066.3869793077383) + + auto linestring1_offsets = make_device_vector({0, 4, 8}); + auto linestring1_points_xy = make_device_vector({41658.902315589876, + 14694.11814724456, + 46600.70359801489, + 8771.431887804214, + 47079.510547637154, + 10199.68027155776, + 51498.48049880379, + 17049.62665643919, + -27429.917796286478, + -33240.8339287343, + -21764.269974046114, + -37974.45515744517, + -14460.71813363161, + -31333.481529957502, + -18226.13032712476, + -30181.03842467982}); + + auto linestring2_offsets = make_device_vector({0, 2, 4}); + auto linestring2_points_xy = make_device_vector({ + 24046.170375947084, + 27878.56737867571, + 20614.007047185743, + 26489.74880629428, + 48381.39607717942, + -8366.313156569413, + 53346.77764665915, + -2066.3869793077383, + }); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({22000.86425379464, 66907.56415814416}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); +} + +TYPED_TEST(PairwiseLinestringDistanceTest, FourPairsSingleLineString) +{ + using T = TypeParam; + auto constexpr num_pairs = 4; + + auto linestring1_offsets = make_device_vector({0, 3, 5, 8, 10}); + + auto linestring1_points_xy = + make_device_vector({0, 1, 1, 0, -1, 0, 0, 0, 0, 1, 0, 0, 2, 2, -2, 0, 2, 2, -2, -2}); + + auto linestring2_offsets = make_device_vector({0, 4, 7, 9, 12}); + auto linestring2_points_xy = make_device_vector( + {1, 1, 2, 1, 2, 0, 3, 0, 1, 0, 1, 1, 1, 2, 2, 0, 0, 2, 1, 1, 5, 5, 10, 0}); + + auto linestring1_points_it = make_vec_2d_iterator(linestring1_points_xy.begin()); + auto linestring2_points_it = make_vec_2d_iterator(linestring2_points_xy.begin()); + + auto expected = make_device_vector({static_cast(std::sqrt(2.0) * 0.5), 1.0, 0.0, 0.0}); + auto got = rmm::device_vector(expected.size()); + + auto mlinestrings1 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring1_offsets.size() - 1, + linestring1_offsets.begin(), + linestring1_points_xy.size() / 2, + linestring1_points_it); + auto mlinestrings2 = make_multilinestring_range(num_pairs, + thrust::make_counting_iterator(0), + linestring2_offsets.size() - 1, + linestring2_offsets.begin(), + linestring2_points_xy.size() / 2, + linestring2_points_it); + + auto ret = pairwise_linestring_distance(mlinestrings1, mlinestrings2, got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(num_pairs, std::distance(got.begin(), ret)); } } // namespace test