From ac872e7eab3c07f6555d2fdb8e70c1d456e8b9c9 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Thu, 11 May 2023 15:53:10 -0700 Subject: [PATCH 01/27] [skip ci] initial segment_counter working --- cpp/include/cuspatial/detail/functors.cuh | 2 +- .../detail/method/segment_method.cuh | 183 ++++++++++++++++++ .../detail/method/segment_method_view.cuh | 92 +++++++++ .../detail/range/multilinestring_range.cuh | 9 + .../cuspatial/range/multilinestring_range.cuh | 5 + cpp/tests/range/multilinestring_range_test.cu | 29 ++- cpp/tests/range/multipolygon_range_test.cu | 4 +- 7 files changed, 318 insertions(+), 6 deletions(-) create mode 100644 cpp/include/cuspatial/detail/method/segment_method.cuh create mode 100644 cpp/include/cuspatial/detail/method/segment_method_view.cuh diff --git a/cpp/include/cuspatial/detail/functors.cuh b/cpp/include/cuspatial/detail/functors.cuh index 1f390753d..75d60061a 100644 --- a/cpp/include/cuspatial/detail/functors.cuh +++ b/cpp/include/cuspatial/detail/functors.cuh @@ -49,7 +49,7 @@ struct offset_pair_to_count_functor { /** * @brief Convert counts of points to counts of segments in a linestring. * - * A Multilinestring is composed of a series of Linestrings. Each Linestring is composed of a + * A Multilinestring is composed of a series of Linestrings. Each Linestring is composed of a series of * segments. The number of segments in a multilinestring is the number of points in the * multilinestring minus the number of linestrings. * diff --git a/cpp/include/cuspatial/detail/method/segment_method.cuh b/cpp/include/cuspatial/detail/method/segment_method.cuh new file mode 100644 index 000000000..8351a924d --- /dev/null +++ b/cpp/include/cuspatial/detail/method/segment_method.cuh @@ -0,0 +1,183 @@ +/* + * Copyright (c) 2023, NVIDIA CORPORATION. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + + +#pragma once + +#include + +#include "segment_method_view.cuh" + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace cuspatial { +namespace detail { + +template +struct equals { + __device__ bool operator()(IndexType x) const { + return x == value; + } +}; + +template +struct greater_than_zero_functor { + __device__ IndexType operator()(IndexType x) const { + return x > 0; + } +}; + +template +class segment_method { + +using index_t = iterator_value_type; + +public: + // segment_methods is always internal use, thus memory consumed is always temporary, + // therefore always use default device memory resource. + segment_method(ParentRange parent_range, rmm::cuda_stream_view stream) : + _range(parent_range), + _non_empty_geometry_prefix_sum(0, stream) { + + auto offset_range = range(_range.part_offset_begin(), _range.part_offset_end()); + auto count_begin = thrust::make_transform_iterator( + thrust::make_zip_iterator(offset_range.begin(), thrust::next(offset_range.begin())), + offset_pair_to_count_functor{}); + + // Preemptive test: does the given range contain any empty ring/linestring? + _contains_empty_geom = thrust::any_of( + rmm::exec_policy(stream), + count_begin, + count_begin + _range.num_linestrings(), + equals{} + ); + + std::cout << std::boolalpha << "contains empty geometry: " << _contains_empty_geom << std::endl; + + if (_contains_empty_geom) + { + + auto count_greater_than_zero = thrust::make_transform_iterator( + count_begin, + greater_than_zero_functor{} + ); + + thrust::device_vector count_greater_than_zero_truth( + count_greater_than_zero, + count_greater_than_zero + _range.num_linestrings()); + + test::print_device_vector( + count_greater_than_zero_truth, + "count_greater_than_zero_truth: " + ); + + // Compute the number of empty linestrings + _non_empty_geometry_prefix_sum.resize( + _range.num_multilinestrings() + 1, stream + ); + + auto key_begin = make_geometry_id_iterator( + _range.geometry_offsets_begin(), + _range.geometry_offsets_end() + ); + + thrust::device_vector key_truth( + key_begin, + key_begin + _range.num_linestrings()); + + test::print_device_vector( + key_truth, + "key_truth: " + ); + + zero_data_async( + _non_empty_geometry_prefix_sum.begin(), + _non_empty_geometry_prefix_sum.end(), + stream + ); + thrust::reduce_by_key( + rmm::exec_policy(stream), + key_begin, + key_begin + _range.num_linestrings(), + count_greater_than_zero, + thrust::make_discard_iterator(), + thrust::next(_non_empty_geometry_prefix_sum.begin()), + thrust::equal_to{}, + thrust::plus{} + ); + + test::print_device_vector( + _non_empty_geometry_prefix_sum, + "non_empty_linestrings: "); + + thrust::inclusive_scan( + rmm::exec_policy(stream), + thrust::next(_non_empty_geometry_prefix_sum.begin()), + _non_empty_geometry_prefix_sum.end(), + thrust::next(_non_empty_geometry_prefix_sum.begin()) + ); + + _num_segments = _range.num_points() - _non_empty_geometry_prefix_sum.element( + _non_empty_geometry_prefix_sum.size() - 1, stream + ); + } else { + _num_segments = _range.num_points() - _range.num_multilinestrings(); + } + + test::print_device_vector( + _non_empty_geometry_prefix_sum, + "non_empty_geometry_prefix_sum: "); + } + + auto view() { + auto index_range = range( + _non_empty_geometry_prefix_sum.begin(), + _non_empty_geometry_prefix_sum.end() + ); + return segment_method_view{ + _range, + index_range, + _num_segments, + _contains_empty_geom + }; + } + +private: + ParentRange _range; + bool _contains_empty_geom; + index_t _num_segments; + rmm::device_uvector _non_empty_geometry_prefix_sum; +}; + +} // namespace detail + + + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/method/segment_method_view.cuh b/cpp/include/cuspatial/detail/method/segment_method_view.cuh new file mode 100644 index 000000000..03d28612b --- /dev/null +++ b/cpp/include/cuspatial/detail/method/segment_method_view.cuh @@ -0,0 +1,92 @@ +/* + * Copyright (c) 2023, NVIDIA CORPORATION. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + + +#pragma once + +#include +#include + +#include +#include + +#include + +namespace cuspatial { +namespace detail { + + +template +class segment_method_view { + +using index_t = typename IndexRange::value_type; + +public: + segment_method_view( + ParentRange range, IndexRange non_empty_geometry_prefix_sum, + index_t num_segments, bool contains_empty_ring + ) : _range(range), + _non_empty_geometry_prefix_sum(non_empty_geometry_prefix_sum), + _num_segments(num_segments), + _contains_empty_ring(contains_empty_ring) + {} + + auto non_empty_linestring_count_begin() { + auto begin = _non_empty_geometry_prefix_sum.begin(); + auto paired_begin = thrust::make_zip_iterator( + begin, thrust::next(begin) + ); + return thrust::make_transform_iterator( + paired_begin, offset_pair_to_count_functor{} + ); + } + + auto segment_count_begin() { + + auto num_geoms = _range.num_multilinestrings(); + auto num_points_begin = _range.multilinestring_point_count_begin(); + auto n_point_linestring_pair_it = thrust::make_zip_iterator( + num_points_begin, this->non_empty_linestring_count_begin() + ); + + return thrust::make_transform_iterator( + n_point_linestring_pair_it, point_count_to_segment_count_functor{} + ); + } + + auto segment_count_end() { + return thrust::next(this->segment_count_begin(), _num_segments); + } + + index_t num_segments() { + return _num_segments; + } + +private: + ParentRange _range; + IndexRange _non_empty_geometry_prefix_sum; + index_t _num_segments; + bool _contains_empty_ring; +}; + +template +segment_method_view( + ParentRange, IndexRange, typename IndexRange::value_type, bool +) -> segment_method_view; + +} // namespace detail + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh index 5983efc00..92123bf5b 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh @@ -30,6 +30,7 @@ #include #include #include +#include #include @@ -280,6 +281,14 @@ multilinestring_range::segment_end( return segment_begin() + num_segments(); } +template +CUSPATIAL_HOST_DEVICE auto +multilinestring_range::segment_methods(rmm::cuda_stream_view stream) +{ + return segment_method{*this, stream}; +} + + template CUSPATIAL_HOST_DEVICE auto multilinestring_range::segment_offset_begin() diff --git a/cpp/include/cuspatial/range/multilinestring_range.cuh b/cpp/include/cuspatial/range/multilinestring_range.cuh index b79b2ae77..1489d688b 100644 --- a/cpp/include/cuspatial/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/range/multilinestring_range.cuh @@ -23,6 +23,8 @@ #include #include +#include + #include namespace cuspatial { @@ -172,6 +174,9 @@ class multilinestring_range { /// Returns an iterator to the end of the segment CUSPATIAL_HOST_DEVICE auto segment_end(); + /// Returns an iterator to the end of the segment + CUSPATIAL_HOST_DEVICE auto segment_methods(rmm::cuda_stream_view); + /// Returns the `multilinestring_idx`th multilinestring in the range. template CUSPATIAL_HOST_DEVICE auto operator[](IndexType multilinestring_idx); diff --git a/cpp/tests/range/multilinestring_range_test.cu b/cpp/tests/range/multilinestring_range_test.cu index 0d00d5c5f..a1cc6e03d 100644 --- a/cpp/tests/range/multilinestring_range_test.cu +++ b/cpp/tests/range/multilinestring_range_test.cu @@ -76,10 +76,10 @@ struct MultilinestringRangeTest : public BaseFixture { { auto multilinestring_array = make_multilinestring_array(geometry_offset, part_offset, coordinates); + auto rng = multilinestring_array.range(); auto d_expected = thrust::device_vector(expected.begin(), expected.end()); - auto rng = multilinestring_array.range(); rmm::device_uvector got(rng.num_multilinestrings(), stream()); thrust::copy(rmm::exec_policy(stream()), @@ -90,6 +90,29 @@ struct MultilinestringRangeTest : public BaseFixture { CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); } + void run_multilinestring_segment_method_count_test(std::initializer_list geometry_offset, + std::initializer_list part_offset, + std::initializer_list> coordinates, + std::initializer_list expected) + { + auto multilinestring_array = + make_multilinestring_array(geometry_offset, part_offset, coordinates); + auto rng = multilinestring_array.range(); + auto multilinestring_with_segment_methods = rng.segment_methods(stream()); + auto methods_view = multilinestring_with_segment_methods.view(); + + auto d_expected = thrust::device_vector(expected.begin(), expected.end()); + + rmm::device_uvector got(rng.num_multilinestrings(), stream()); + thrust::copy(rmm::exec_policy(stream()), + methods_view.segment_count_begin(), + methods_view.segment_count_end(), + got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); + } + + void run_multilinestring_linestring_count_test(std::initializer_list geometry_offset, std::initializer_list part_offset, std::initializer_list> coordinates, @@ -312,13 +335,13 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest4) } // FIXME: contains empty linestring -TYPED_TEST(MultilinestringRangeTest, DISABLED_MultilinestringSegmentCountTest5) +TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest5) { using T = TypeParam; using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, {0, 1, 2, 3}, {0, 3, 3, 6}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, diff --git a/cpp/tests/range/multipolygon_range_test.cu b/cpp/tests/range/multipolygon_range_test.cu index 0f1b08809..6929c5136 100644 --- a/cpp/tests/range/multipolygon_range_test.cu +++ b/cpp/tests/range/multipolygon_range_test.cu @@ -345,7 +345,7 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount4) {9, 3}); } -// FIXME: multipolygon doesn't constructor doesn't allow empty rings, should it? +// FIXME: multipolygon constructor doesn't allow empty rings, should it? TYPED_TEST(MultipolygonRangeTest, DISABLED_MultipolygonSegmentCount_ConatainsEmptyRing) { CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, @@ -367,7 +367,7 @@ TYPED_TEST(MultipolygonRangeTest, DISABLED_MultipolygonSegmentCount_ConatainsEmp {6, 3}); } -// FIXME: multipolygon doesn't constructor doesn't allow empty rings, should it? +// FIXME: multipolygon constructor doesn't allow empty rings, should it? TYPED_TEST(MultipolygonRangeTest, DISABLED_MultipolygonSegmentCount_ConatainsEmptyPart) { CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, From 84f317c52b3ac3121cdfb5565dec9d045634fa2f Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 16 May 2023 11:19:55 -0700 Subject: [PATCH 02/27] segment iterator works with empty linestrings --- cpp/include/cuspatial/detail/functors.cuh | 42 +++--- .../detail/method/segment_method_view.cuh | 129 ++++++++++-------- .../detail/range/multilinestring_range.cuh | 7 +- .../detail/range/multipolygon_range.cuh | 6 +- cpp/tests/range/multilinestring_range_test.cu | 28 +++- 5 files changed, 128 insertions(+), 84 deletions(-) diff --git a/cpp/include/cuspatial/detail/functors.cuh b/cpp/include/cuspatial/detail/functors.cuh index 75d60061a..e48a1082e 100644 --- a/cpp/include/cuspatial/detail/functors.cuh +++ b/cpp/include/cuspatial/detail/functors.cuh @@ -75,8 +75,12 @@ struct point_count_to_segment_count_functor { }; /** - * @brief Given an offset iterator it, returns an iterator of the distance between it and an input - * index i + * @brief Given an offset iterator it that partitions a point range, return an offset iterator that + * partitions the segment range made from the same point range. + * + * One partition to a point range introduces one invalid segment, except empty partitions. + * Therefore, the offsets that partitions the segment range is the offset that partitions the point + * range subtracts the number of *non-empty* point partitions that precedes the current point range. * * @tparam OffsetIterator Iterator type to the offset * @@ -87,20 +91,21 @@ struct point_count_to_segment_count_functor { * * Used to create iterator to segment offsets, such as `segment_offset_begin`. */ -template -struct to_distance_iterator { - OffsetIterator begin; +template +struct to_segment_offset_iterator { + OffsetIterator point_partition_begin; + CountIterator non_empty_partitions_begin; template CUSPATIAL_HOST_DEVICE auto operator()(IndexType i) { - return begin[i] - i; + return point_partition_begin[i] - non_empty_partitions_begin[i]; } }; /// Deduction guide for to_distance_iterator -template -to_distance_iterator(OffsetIterator) -> to_distance_iterator; +template +to_segment_offset_iterator(OffsetIterator, CountIterator) -> to_segment_offset_iterator; /** * @brief Return a segment from the a partitioned range of points @@ -114,29 +119,30 @@ to_distance_iterator(OffsetIterator) -> to_distance_iterator; * @tparam OffsetIterator the iterator type indicating partitions of the point range. * @tparam CoordinateIterator the iterator type to the point range. */ -template +template struct to_valid_segment_functor { using element_t = iterator_vec_base_type; - OffsetIterator begin; - OffsetIterator end; + OffsetIterator segment_offset_begin; + OffsetIterator segment_offset_end; + CountIterator non_empty_partitions_begin; CoordinateIterator point_begin; template - CUSPATIAL_HOST_DEVICE segment operator()(IndexType i) + CUSPATIAL_HOST_DEVICE segment operator()(IndexType sid) { - auto kit = thrust::upper_bound(thrust::seq, begin, end, i); - auto k = thrust::distance(begin, kit); - auto pid = i + k - 1; + auto kit = thrust::upper_bound(thrust::seq, segment_offset_begin, segment_offset_end, sid); + auto k = thrust::distance(segment_offset_begin, kit); + auto pid = non_empty_partitions_begin[sid] + k - 1; return segment{point_begin[pid], point_begin[pid + 1]}; } }; /// Deduction guide for to_valid_segment_functor -template -to_valid_segment_functor(OffsetIterator, OffsetIterator, CoordinateIterator) - -> to_valid_segment_functor; +template +to_valid_segment_functor(OffsetIterator, OffsetIterator, CountIterator, CoordinateIterator) + -> to_valid_segment_functor; } // namespace detail } // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/method/segment_method_view.cuh b/cpp/include/cuspatial/detail/method/segment_method_view.cuh index 03d28612b..d2d0a7d89 100644 --- a/cpp/include/cuspatial/detail/method/segment_method_view.cuh +++ b/cpp/include/cuspatial/detail/method/segment_method_view.cuh @@ -14,78 +14,93 @@ * limitations under the License. */ - #pragma once -#include +#include + #include +#include +#include #include #include +#include #include namespace cuspatial { namespace detail { - -template +template class segment_method_view { - -using index_t = typename IndexRange::value_type; - -public: - segment_method_view( - ParentRange range, IndexRange non_empty_geometry_prefix_sum, - index_t num_segments, bool contains_empty_ring - ) : _range(range), - _non_empty_geometry_prefix_sum(non_empty_geometry_prefix_sum), - _num_segments(num_segments), - _contains_empty_ring(contains_empty_ring) - {} - - auto non_empty_linestring_count_begin() { - auto begin = _non_empty_geometry_prefix_sum.begin(); - auto paired_begin = thrust::make_zip_iterator( - begin, thrust::next(begin) - ); - return thrust::make_transform_iterator( - paired_begin, offset_pair_to_count_functor{} - ); - } - - auto segment_count_begin() { - - auto num_geoms = _range.num_multilinestrings(); - auto num_points_begin = _range.multilinestring_point_count_begin(); - auto n_point_linestring_pair_it = thrust::make_zip_iterator( - num_points_begin, this->non_empty_linestring_count_begin() - ); - - return thrust::make_transform_iterator( - n_point_linestring_pair_it, point_count_to_segment_count_functor{} - ); - } - - auto segment_count_end() { - return thrust::next(this->segment_count_begin(), _num_segments); - } - - index_t num_segments() { - return _num_segments; - } - -private: - ParentRange _range; - IndexRange _non_empty_geometry_prefix_sum; - index_t _num_segments; - bool _contains_empty_ring; + using index_t = typename IndexRange::value_type; + + public: + segment_method_view(ParentRange range, + IndexRange non_empty_geometry_prefix_sum, + index_t num_segments, + bool contains_empty_ring) + : _range(range), + _non_empty_geometry_prefix_sum(non_empty_geometry_prefix_sum), + _num_segments(num_segments), + _contains_empty_ring(contains_empty_ring) + { + } + + auto non_empty_linestring_count_begin() + { + auto begin = _non_empty_geometry_prefix_sum.begin(); + auto paired_begin = thrust::make_zip_iterator(begin, thrust::next(begin)); + return thrust::make_transform_iterator(paired_begin, offset_pair_to_count_functor{}); + } + + auto segment_count_begin() + { + auto num_geoms = _range.num_multilinestrings(); + auto num_points_begin = _range.multilinestring_point_count_begin(); + auto n_point_linestring_pair_it = + thrust::make_zip_iterator(num_points_begin, this->non_empty_linestring_count_begin()); + + return thrust::make_transform_iterator(n_point_linestring_pair_it, + point_count_to_segment_count_functor{}); + } + + auto segment_count_end() { return thrust::next(this->segment_count_begin(), _num_segments); } + + index_t num_segments() { return _num_segments; } + + auto segment_offset_begin() + { + return make_counting_transform_iterator( + 0, + to_segment_offset_iterator{_range.part_offset_begin(), + _non_empty_geometry_prefix_sum.begin()}); + } + + auto segment_offset_end() { return segment_offset_begin() + num_segments(); } + + auto segment_begin() + { + return make_counting_transform_iterator( + 0, + to_valid_segment_functor{segment_offset_begin(), + segment_offset_end(), + _non_empty_geometry_prefix_sum.begin(), + _range.point_begin()}); + } + + auto segment_end() { return segment_begin() + _num_segments; } + + private: + ParentRange _range; + IndexRange _non_empty_geometry_prefix_sum; + index_t _num_segments; + bool _contains_empty_ring; }; -template -segment_method_view( - ParentRange, IndexRange, typename IndexRange::value_type, bool -) -> segment_method_view; +template +segment_method_view(ParentRange, IndexRange, typename IndexRange::value_type, bool) + -> segment_method_view; } // namespace detail diff --git a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh index 92123bf5b..c0ac4d5e2 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh @@ -16,6 +16,7 @@ #pragma once +#include "thrust/iterator/counting_iterator.h" #include #include #include @@ -271,7 +272,7 @@ multilinestring_range::segment_begi return detail::make_counting_transform_iterator( 0, detail::to_valid_segment_functor{ - this->segment_offset_begin(), this->segment_offset_end(), _point_begin}); + this->segment_offset_begin(), this->segment_offset_end(), thrust::make_counting_iterator(0), _point_begin}); } template @@ -288,12 +289,12 @@ multilinestring_range::segment_meth return segment_method{*this, stream}; } - template CUSPATIAL_HOST_DEVICE auto multilinestring_range::segment_offset_begin() { - return detail::make_counting_transform_iterator(0, detail::to_distance_iterator{_part_begin}); + return detail::make_counting_transform_iterator(0, detail::to_segment_offset_iterator{ + _part_begin, thrust::make_counting_iterator(0)}); } template diff --git a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh index 85781e55b..ffc04f12e 100644 --- a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh @@ -16,6 +16,7 @@ #pragma once +#include "thrust/iterator/counting_iterator.h" #include #include #include @@ -397,7 +398,7 @@ multipolygon_range::s return detail::make_counting_transform_iterator( 0, detail::to_valid_segment_functor{ - this->subtracted_ring_begin(), this->subtracted_ring_end(), _point_begin}); + this->subtracted_ring_begin(), this->subtracted_ring_end(), thrust::make_counting_iterator(0), _point_begin}); } template :: subtracted_ring_begin() { - return detail::make_counting_transform_iterator(0, detail::to_distance_iterator{_ring_begin}); + return detail::make_counting_transform_iterator(0, detail::to_segment_offset_iterator{ + _ring_begin, thrust::make_counting_iterator(0)}); } template geometry_offset, + std::initializer_list part_offset, + std::initializer_list> coordinates, + std::initializer_list> expected) + { + auto multilinestring_array = + make_multilinestring_array(geometry_offset, part_offset, coordinates); + auto rng = multilinestring_array.range(); + auto segment_methods = rng.segment_methods(stream()); + auto segment_view = segment_methods.view(); + + rmm::device_uvector> got(segment_view.num_segments(), stream()); + thrust::copy(rmm::exec_policy(stream()), segment_view.segment_begin(), segment_view.segment_end(), got.begin()); + + auto d_expected = thrust::device_vector>(expected.begin(), expected.end()); + CUSPATIAL_EXPECT_VEC2D_PAIRS_EQUIVALENT(d_expected, got); + } + + + void run_multilinestring_point_count_test(std::initializer_list geometry_offset, std::initializer_list part_offset, std::initializer_list> coordinates, @@ -244,14 +264,14 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorManyPairTest) } /// FIXME: Currently, segment iterator doesn't handle empty linestrings. -TYPED_TEST(MultilinestringRangeTest, DISABLED_SegmentIteratorWithEmptyLineTest) +TYPED_TEST(MultilinestringRangeTest, SegmentIteratorWithEmptyLineTest) { using T = TypeParam; using P = vec_2d; using S = segment; CUSPATIAL_RUN_TEST( - this->run_segment_test_single, + this->run_segment_test_with_method_single, {0, 1, 2, 3}, {0, 3, 3, 6}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, @@ -295,13 +315,13 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest) } // FIXME: contains empty linestring -TYPED_TEST(MultilinestringRangeTest, DISABLED_MultilinestringSegmentCountTest2) +TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest2) { using T = TypeParam; using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, {0, 1, 3}, {0, 3, 3, 6}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, From ce28b3c7e13a5eb9ddc76d13619a4824af0ce114 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 16 May 2023 14:26:12 -0700 Subject: [PATCH 03/27] enable polygon segment count test --- cpp/tests/range/multipolygon_range_test.cu | 40 +++++++++++++++++++--- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/cpp/tests/range/multipolygon_range_test.cu b/cpp/tests/range/multipolygon_range_test.cu index 6929c5136..5a6ea0c00 100644 --- a/cpp/tests/range/multipolygon_range_test.cu +++ b/cpp/tests/range/multipolygon_range_test.cu @@ -100,6 +100,32 @@ struct MultipolygonRangeTest : public BaseFixture { CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(got, d_expected); } + void run_multipolygon_segment_method_count_single( + std::initializer_list geometry_offset, + std::initializer_list part_offset, + std::initializer_list ring_offset, + std::initializer_list> coordinates, + std::initializer_list expected_segment_counts) + { + auto multipolygon_array = + make_multipolygon_array(geometry_offset, part_offset, ring_offset, coordinates); + auto rng = multipolygon_array.range(); + auto segment_methods = rng.segment_methods(stream()); + auto segment_view = segment_methods.view(); + + auto got = rmm::device_uvector(rng.num_multipolygons(), stream()); + + thrust::copy(rmm::exec_policy(stream()), + segment_view.segment_count_begin(), + segment_view.segment_count_end(), + got.begin()); + + auto d_expected = thrust::device_vector(expected_segment_counts.begin(), + expected_segment_counts.end()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(got, d_expected); + } + void test_multipolygon_as_multilinestring( std::initializer_list multipolygon_geometry_offset, std::initializer_list multipolygon_part_offset, @@ -346,15 +372,18 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount4) } // FIXME: multipolygon constructor doesn't allow empty rings, should it? -TYPED_TEST(MultipolygonRangeTest, DISABLED_MultipolygonSegmentCount_ConatainsEmptyRing) +TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount_ConatainsEmptyRing) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_count_single, {0, 2, 3}, {0, 2, 3, 4}, - {0, 4, 4, 8, 12}, + {0, 7, 7, 11, 18}, {{0, 0}, {1, 0}, {1, 1}, + {0.5, 1.5}, + {0, 1.0}, + {0.5, 0.5}, {0, 0}, {0.2, 0.2}, {0.2, 0.3}, @@ -363,8 +392,11 @@ TYPED_TEST(MultipolygonRangeTest, DISABLED_MultipolygonSegmentCount_ConatainsEmp {0, 0}, {1, 0}, {1, 1}, + {0.5, 1.5}, + {0, 1.0}, + {0.5, 0.5}, {0, 0}}, - {6, 3}); + {9, 6}); } // FIXME: multipolygon constructor doesn't allow empty rings, should it? From cd2d229c45b0d83220eadc15d82793ada6ff79f1 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 16 May 2023 14:26:51 -0700 Subject: [PATCH 04/27] add polygon segment methods --- cpp/include/cuspatial/range/multipolygon_range.cuh | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/cpp/include/cuspatial/range/multipolygon_range.cuh b/cpp/include/cuspatial/range/multipolygon_range.cuh index 1b8947345..3619e8902 100644 --- a/cpp/include/cuspatial/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/range/multipolygon_range.cuh @@ -16,7 +16,6 @@ #pragma once -#include #include #include @@ -24,6 +23,11 @@ #include #include +#include + +#include + + namespace cuspatial { /** @@ -117,6 +121,12 @@ class multipolygon_range { /// Return the iterator to the one past the last point in the range. CUSPATIAL_HOST_DEVICE auto point_end(); + /// Return the iterator to the first geometry offset in the range. + CUSPATIAL_HOST_DEVICE auto geometry_offsets_begin() { return _part_begin; } + + /// Return the iterator to the one past the last geometry offset in the range. + CUSPATIAL_HOST_DEVICE auto geometry_offsets_end() { return _part_end; } + /// Return the iterator to the first part offset in the range. CUSPATIAL_HOST_DEVICE auto part_offset_begin() { return _part_begin; } @@ -186,6 +196,8 @@ class multipolygon_range { /// Returns an iterator to the end of the segment CUSPATIAL_HOST_DEVICE auto segment_end(); + CUSPATIAL_HOST_DEVICE auto segment_methods(rmm::cuda_stream_view); + /// Range Casting /// Cast the range of multipolygons as a range of multipoints, ignoring all edge connections and From fcbe105dabaf3bca3df064754ee66c9d5332deba Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 17 May 2023 09:20:47 -0700 Subject: [PATCH 05/27] remove incomplete validations and add implicit assumptions in documentation --- .../cuspatial/detail/utility/validation.hpp | 28 ++++++++----------- 1 file changed, 12 insertions(+), 16 deletions(-) diff --git a/cpp/include/cuspatial/detail/utility/validation.hpp b/cpp/include/cuspatial/detail/utility/validation.hpp index 1a4bd2001..e66d8dd08 100644 --- a/cpp/include/cuspatial/detail/utility/validation.hpp +++ b/cpp/include/cuspatial/detail/utility/validation.hpp @@ -36,9 +36,7 @@ */ #define CUSPATIAL_EXPECTS_VALID_LINESTRING_SIZES(num_linestring_points, num_linestring_offsets) \ CUSPATIAL_HOST_DEVICE_EXPECTS(num_linestring_offsets > 0, \ - "Polygon offsets must contain at least one (1) value"); \ - CUSPATIAL_HOST_DEVICE_EXPECTS(num_linestring_points >= 2 * (num_linestring_offsets - 1), \ - "Each linestring must have at least two vertices"); + "Polygon offsets must contain at least one (1) value"); /** * @brief Macro for validating the data array sizes for multilinestrings. @@ -69,8 +67,6 @@ * Raises an exception if any of the following are false: * - The number of polygon offsets is greater than zero. * - The number of ring offsets is greater than zero. - * - There is at least one ring offset per polygon offset. - * - There are at least four vertices per ring offset. * * Polygons follow [GeoArrow data layout][1]. Offsets arrays (polygons and rings) have one more * element than the number of items in the array. The last offset is always the sum of the previous @@ -78,8 +74,11 @@ * last ring offset plus the number of rings in the last polygon. See * [Arrow Variable-Size Binary layout](2). Note that an empty list still has one offset: {0}. * - * Rings are assumed to be closed (closed means the first and last vertices of - * each ring are equal). Therefore rings must have at least 4 vertices. + * The following are not explicitly checked in this macro, but is always assumed in cuspatial: + * + * 1. Polygon can contain zero (0) or more rings. A polygon with zero rings is an empty polygon. + * 2. Rings are assumed to be closed (closed means the first and last vertices of each ring are equal). + * Rings can also be empty. Therefore each ring must contain zero (0) or four (4) or more coordinates. * * [1]: https://github.com/geoarrow/geoarrow/blob/main/format.md * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout @@ -89,11 +88,7 @@ CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_offsets > 0, \ "Polygon offsets must contain at least one (1) value"); \ CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_ring_offsets > 0, \ - "Polygon ring offsets must contain at least one (1) value"); \ - CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_ring_offsets >= num_poly_offsets, \ - "Each polygon must have at least one (1) ring"); \ - CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_points >= 4 * (num_poly_ring_offsets - 1), \ - "Each ring must have at least four (4) vertices"); + "Polygon ring offsets must contain at least one (1) value"); /** * @brief Macro for validating the data array sizes for a multipolygon. @@ -102,8 +97,6 @@ * - The number of multipolygon offsets is greater than zero. * - The number of polygon offsets is greater than zero. * - The number of ring offsets is greater than zero. - * - There is at least one ring offset per polygon offset. - * - There are at least four vertices per ring offset. * * MultiPolygons follow [GeoArrow data layout][1]. Offsets arrays (polygons and rings) have one more * element than the number of items in the array. The last offset is always the sum of the previous @@ -111,8 +104,11 @@ * last ring offset plus the number of rings in the last polygon. See * [Arrow Variable-Size Binary layout](2). Note that an empty list still has one offset: {0}. * - * Rings are assumed to be closed (closed means the first and last vertices of - * each ring are equal). Therefore rings must have at least 4 vertices. + * The following are not explicitly checked in this macro, but is always assumed in cuspatial: + * + * 1. Polygon can contain zero (0) or more rings. A polygon with zero rings is an empty polygon. + * 2. Rings are assumed to be closed (closed means the first and last vertices of each ring are equal). + * Rings can also be empty. Therefore each ring must contain zero (0) or four (4) or more coordinates. * * [1]: https://github.com/geoarrow/geoarrow/blob/main/format.md * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout From ee3576ed9acd6a34b77979e8517a754a2a330173 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 17 May 2023 09:21:10 -0700 Subject: [PATCH 06/27] enable multipolygon range tests --- .../detail/method/segment_method.cuh | 28 ++++++++----------- .../detail/method/segment_method_view.cuh | 5 ++-- .../detail/range/multipolygon_range.cuh | 11 ++++++++ cpp/tests/range/multipolygon_range_test.cu | 8 +++--- 4 files changed, 29 insertions(+), 23 deletions(-) diff --git a/cpp/include/cuspatial/detail/method/segment_method.cuh b/cpp/include/cuspatial/detail/method/segment_method.cuh index 8351a924d..5b4770bf7 100644 --- a/cpp/include/cuspatial/detail/method/segment_method.cuh +++ b/cpp/include/cuspatial/detail/method/segment_method.cuh @@ -62,7 +62,7 @@ using index_t = iterator_value_type; public: // segment_methods is always internal use, thus memory consumed is always temporary, // therefore always use default device memory resource. - segment_method(ParentRange parent_range, rmm::cuda_stream_view stream) : + segment_method(ParentRange parent_range, rmm::cuda_stream_view stream) : _range(parent_range), _non_empty_geometry_prefix_sum(0, stream) { @@ -71,18 +71,15 @@ public: thrust::make_zip_iterator(offset_range.begin(), thrust::next(offset_range.begin())), offset_pair_to_count_functor{}); - // Preemptive test: does the given range contain any empty ring/linestring? - _contains_empty_geom = thrust::any_of( - rmm::exec_policy(stream), - count_begin, - count_begin + _range.num_linestrings(), - equals{} - ); - - std::cout << std::boolalpha << "contains empty geometry: " << _contains_empty_geom << std::endl; + // // Preemptive test: does the given range contain any empty ring/linestring? + // _contains_empty_geom = thrust::any_of( + // rmm::exec_policy(stream), + // count_begin, + // count_begin + _range.num_linestrings(), + // equals{} + // ); - if (_contains_empty_geom) - { + // std::cout << std::boolalpha << "contains empty geometry: " << _contains_empty_geom << std::endl; auto count_greater_than_zero = thrust::make_transform_iterator( count_begin, @@ -102,7 +99,7 @@ public: _non_empty_geometry_prefix_sum.resize( _range.num_multilinestrings() + 1, stream ); - + auto key_begin = make_geometry_id_iterator( _range.geometry_offsets_begin(), _range.geometry_offsets_end() @@ -125,7 +122,7 @@ public: thrust::reduce_by_key( rmm::exec_policy(stream), key_begin, - key_begin + _range.num_linestrings(), + key_begin + _range.num_linestrings(), count_greater_than_zero, thrust::make_discard_iterator(), thrust::next(_non_empty_geometry_prefix_sum.begin()), @@ -147,9 +144,6 @@ public: _num_segments = _range.num_points() - _non_empty_geometry_prefix_sum.element( _non_empty_geometry_prefix_sum.size() - 1, stream ); - } else { - _num_segments = _range.num_points() - _range.num_multilinestrings(); - } test::print_device_vector( _non_empty_geometry_prefix_sum, diff --git a/cpp/include/cuspatial/detail/method/segment_method_view.cuh b/cpp/include/cuspatial/detail/method/segment_method_view.cuh index d2d0a7d89..0454b1ae4 100644 --- a/cpp/include/cuspatial/detail/method/segment_method_view.cuh +++ b/cpp/include/cuspatial/detail/method/segment_method_view.cuh @@ -56,7 +56,6 @@ class segment_method_view { auto segment_count_begin() { - auto num_geoms = _range.num_multilinestrings(); auto num_points_begin = _range.multilinestring_point_count_begin(); auto n_point_linestring_pair_it = thrust::make_zip_iterator(num_points_begin, this->non_empty_linestring_count_begin()); @@ -65,7 +64,9 @@ class segment_method_view { point_count_to_segment_count_functor{}); } - auto segment_count_end() { return thrust::next(this->segment_count_begin(), _num_segments); } + auto segment_count_end() { + std::cout << "num multilinestrings: " << _range.num_multilinestrings() << std::endl; + return thrust::next(this->segment_count_begin(), _range.num_multilinestrings()); } index_t num_segments() { return _num_segments; } diff --git a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh index ffc04f12e..6f1ae66eb 100644 --- a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh @@ -28,6 +28,7 @@ #include #include +#include #include #include #include @@ -433,6 +434,16 @@ multipolygon_range::s return subtracted_ring_begin() + thrust::distance(_ring_begin, _ring_end); } +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::segment_methods(rmm::cuda_stream_view stream){ + auto multilinestring_range = this->as_multilinestring_range(); + return segment_method{multilinestring_range, stream}; +} + template run_multipolygon_segment_method_count_single, {0, 2, 3}, @@ -400,11 +400,11 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount_ConatainsEmptyRing) } // FIXME: multipolygon constructor doesn't allow empty rings, should it? -TYPED_TEST(MultipolygonRangeTest, DISABLED_MultipolygonSegmentCount_ConatainsEmptyPart) +TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount_ContainsEmptyPart) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_count_single, {0, 3, 4}, - {0, 2, 2, 3, 4}, + {0, 1, 1, 2, 3}, {0, 4, 8, 12}, {{0, 0}, {1, 0}, From 7a53d057d557748bc6f2ab06dfe11a52441d59d9 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Thu, 18 May 2023 10:04:25 -0700 Subject: [PATCH 07/27] improve doc of validators --- .../cuspatial/detail/utility/validation.hpp | 28 +++++++++--------- cpp/tests/range/multilinestring_range_test.cu | 29 ++----------------- 2 files changed, 18 insertions(+), 39 deletions(-) diff --git a/cpp/include/cuspatial/detail/utility/validation.hpp b/cpp/include/cuspatial/detail/utility/validation.hpp index e66d8dd08..9ce142317 100644 --- a/cpp/include/cuspatial/detail/utility/validation.hpp +++ b/cpp/include/cuspatial/detail/utility/validation.hpp @@ -23,7 +23,6 @@ * * Raises an exception if any of the following are false: * - The number of linestring offsets is greater than zero. - * - There are at least two vertices per linestring offset. * * Linestrings follow [GeoArrow data layout][1]. Offsets arrays have one more element than the * number of items in the array. The last offset is always the sum of the previous offset and the @@ -31,6 +30,9 @@ * last linestring offset plus one. See [Arrow Variable-Size Binary layout](2). Note that an * empty list still has one offset: {0}. * + * The following is not explicitly checked in this macro, but is always assumed in cuspatial: + * 1. LineString can contain zero or two or more vertices. + * * [1]: https://github.com/geoarrow/geoarrow/blob/main/format.md * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout */ @@ -44,7 +46,6 @@ * Raises an exception if any of the following are false: * - The number of multilinestring offsets is greater than zero. * - The number of linestring offsets is greater than zero. - * - There are at least two vertices per linestring offset. * * Multilinestrings follow [GeoArrow data layout][1]. Offsets arrays have one more element than the * number of items in the array. The last offset is always the sum of the previous offset and the @@ -76,18 +77,18 @@ * * The following are not explicitly checked in this macro, but is always assumed in cuspatial: * - * 1. Polygon can contain zero (0) or more rings. A polygon with zero rings is an empty polygon. - * 2. Rings are assumed to be closed (closed means the first and last vertices of each ring are equal). - * Rings can also be empty. Therefore each ring must contain zero (0) or four (4) or more coordinates. + * 1. Polygon can contain zero or more rings. A polygon with zero rings is an empty polygon. + * 2. Rings are assumed to be closed (closed means the first and last vertices of each ring are + * equal). Rings can also be empty. Therefore each ring must contain zero or four or more vertices. * * [1]: https://github.com/geoarrow/geoarrow/blob/main/format.md * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout */ -#define CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES( \ - num_poly_points, num_poly_offsets, num_poly_ring_offsets) \ - CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_offsets > 0, \ - "Polygon offsets must contain at least one (1) value"); \ - CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_ring_offsets > 0, \ +#define CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES( \ + num_poly_points, num_poly_offsets, num_poly_ring_offsets) \ + CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_offsets > 0, \ + "Polygon offsets must contain at least one (1) value"); \ + CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_ring_offsets > 0, \ "Polygon ring offsets must contain at least one (1) value"); /** @@ -106,9 +107,10 @@ * * The following are not explicitly checked in this macro, but is always assumed in cuspatial: * - * 1. Polygon can contain zero (0) or more rings. A polygon with zero rings is an empty polygon. - * 2. Rings are assumed to be closed (closed means the first and last vertices of each ring are equal). - * Rings can also be empty. Therefore each ring must contain zero (0) or four (4) or more coordinates. + * 1. Polygon can contain zero or more rings. A polygon with zero rings is an empty polygon. + * 2. Rings are assumed to be closed (closed means the first and last vertices of each ring are + * equal). Rings can also be empty. Therefore each ring must contain zero or four or more + * coordinates. * * [1]: https://github.com/geoarrow/geoarrow/blob/main/format.md * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout diff --git a/cpp/tests/range/multilinestring_range_test.cu b/cpp/tests/range/multilinestring_range_test.cu index 0119043b4..78b1554a3 100644 --- a/cpp/tests/range/multilinestring_range_test.cu +++ b/cpp/tests/range/multilinestring_range_test.cu @@ -89,27 +89,6 @@ struct MultilinestringRangeTest : public BaseFixture { CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); } - void run_multilinestring_segment_count_test(std::initializer_list geometry_offset, - std::initializer_list part_offset, - std::initializer_list> coordinates, - std::initializer_list expected) - { - auto multilinestring_array = - make_multilinestring_array(geometry_offset, part_offset, coordinates); - auto rng = multilinestring_array.range(); - - auto d_expected = thrust::device_vector(expected.begin(), expected.end()); - - - rmm::device_uvector got(rng.num_multilinestrings(), stream()); - thrust::copy(rmm::exec_policy(stream()), - rng.multilinestring_segment_count_begin(), - rng.multilinestring_segment_count_end(), - got.begin()); - - CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); - } - void run_multilinestring_segment_method_count_test(std::initializer_list geometry_offset, std::initializer_list part_offset, std::initializer_list> coordinates, @@ -263,7 +242,6 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorManyPairTest) S{P{21, 21}, P{22, 22}}}); } -/// FIXME: Currently, segment iterator doesn't handle empty linestrings. TYPED_TEST(MultilinestringRangeTest, SegmentIteratorWithEmptyLineTest) { using T = TypeParam; @@ -311,10 +289,9 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest) using S = segment; CUSPATIAL_RUN_TEST( - this->run_multilinestring_segment_count_test, {0, 1}, {0, 3}, {P{0, 0}, P{1, 1}, P{2, 2}}, {2}); + this->run_multilinestring_segment_method_count_test, {0, 1}, {0, 3}, {P{0, 0}, P{1, 1}, P{2, 2}}, {2}); } -// FIXME: contains empty linestring TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest2) { using T = TypeParam; @@ -334,7 +311,7 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest3) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, {0, 1, 2}, {0, 3, 6}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, @@ -347,7 +324,7 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest4) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, {0, 1, 3}, {0, 3, 5, 7}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{20, 20}, P{21, 21}}, From 2181e966f63743fec26f0cebaa3a725aabde9e3c Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Thu, 18 May 2023 10:08:12 -0700 Subject: [PATCH 08/27] replace all multilinestring range tests with segment method test --- .../detail/method/segment_method.cuh | 212 +++++++----------- .../detail/method/segment_method_view.cuh | 70 ++++-- cpp/tests/range/multilinestring_range_test.cu | 93 ++++---- 3 files changed, 188 insertions(+), 187 deletions(-) diff --git a/cpp/include/cuspatial/detail/method/segment_method.cuh b/cpp/include/cuspatial/detail/method/segment_method.cuh index 5b4770bf7..e35c09b45 100644 --- a/cpp/include/cuspatial/detail/method/segment_method.cuh +++ b/cpp/include/cuspatial/detail/method/segment_method.cuh @@ -14,164 +14,124 @@ * limitations under the License. */ - #pragma once #include #include "segment_method_view.cuh" -#include -#include #include #include #include +#include +#include -#include #include +#include #include #include #include -#include #include +#include #include namespace cuspatial { namespace detail { -template +template struct equals { - __device__ bool operator()(IndexType x) const { - return x == value; - } + __device__ bool operator()(IndexType x) const { return x == value; } }; -template +template struct greater_than_zero_functor { - __device__ IndexType operator()(IndexType x) const { - return x > 0; - } + __device__ IndexType operator()(IndexType x) const { return x > 0; } }; -template +template class segment_method { - -using index_t = iterator_value_type; - -public: - // segment_methods is always internal use, thus memory consumed is always temporary, - // therefore always use default device memory resource. - segment_method(ParentRange parent_range, rmm::cuda_stream_view stream) : - _range(parent_range), - _non_empty_geometry_prefix_sum(0, stream) { - - auto offset_range = range(_range.part_offset_begin(), _range.part_offset_end()); - auto count_begin = thrust::make_transform_iterator( - thrust::make_zip_iterator(offset_range.begin(), thrust::next(offset_range.begin())), - offset_pair_to_count_functor{}); - - // // Preemptive test: does the given range contain any empty ring/linestring? - // _contains_empty_geom = thrust::any_of( - // rmm::exec_policy(stream), - // count_begin, - // count_begin + _range.num_linestrings(), - // equals{} - // ); - - // std::cout << std::boolalpha << "contains empty geometry: " << _contains_empty_geom << std::endl; - - auto count_greater_than_zero = thrust::make_transform_iterator( - count_begin, - greater_than_zero_functor{} - ); - - thrust::device_vector count_greater_than_zero_truth( - count_greater_than_zero, - count_greater_than_zero + _range.num_linestrings()); - - test::print_device_vector( - count_greater_than_zero_truth, - "count_greater_than_zero_truth: " - ); - - // Compute the number of empty linestrings - _non_empty_geometry_prefix_sum.resize( - _range.num_multilinestrings() + 1, stream - ); - - auto key_begin = make_geometry_id_iterator( - _range.geometry_offsets_begin(), - _range.geometry_offsets_end() - ); - - thrust::device_vector key_truth( - key_begin, - key_begin + _range.num_linestrings()); - - test::print_device_vector( - key_truth, - "key_truth: " - ); - - zero_data_async( - _non_empty_geometry_prefix_sum.begin(), - _non_empty_geometry_prefix_sum.end(), - stream - ); - thrust::reduce_by_key( - rmm::exec_policy(stream), - key_begin, - key_begin + _range.num_linestrings(), - count_greater_than_zero, - thrust::make_discard_iterator(), - thrust::next(_non_empty_geometry_prefix_sum.begin()), - thrust::equal_to{}, - thrust::plus{} - ); - - test::print_device_vector( - _non_empty_geometry_prefix_sum, - "non_empty_linestrings: "); - - thrust::inclusive_scan( - rmm::exec_policy(stream), - thrust::next(_non_empty_geometry_prefix_sum.begin()), - _non_empty_geometry_prefix_sum.end(), - thrust::next(_non_empty_geometry_prefix_sum.begin()) - ); - - _num_segments = _range.num_points() - _non_empty_geometry_prefix_sum.element( - _non_empty_geometry_prefix_sum.size() - 1, stream - ); - - test::print_device_vector( - _non_empty_geometry_prefix_sum, - "non_empty_geometry_prefix_sum: "); + using index_t = iterator_value_type; + + public: + // segment_methods is always internal use, thus memory consumed is always temporary, + // therefore always use default device memory resource. + segment_method(ParentRange parent_range, rmm::cuda_stream_view stream) + : _range(parent_range), + _non_empty_linestring_prefix_sum(parent_range.num_linestrings() + 1, stream) + { + auto offset_range = range(_range.part_offset_begin(), _range.part_offset_end()); + auto count_begin = thrust::make_transform_iterator( + thrust::make_zip_iterator(offset_range.begin(), thrust::next(offset_range.begin())), + offset_pair_to_count_functor{}); + + // // Preemptive test: does the given range contain any empty ring/linestring? + // _contains_empty_geom = thrust::any_of( + // rmm::exec_policy(stream), + // count_begin, + // count_begin + _range.num_linestrings(), + // equals{} + // ); + + // std::cout << std::boolalpha << "contains empty geometry: " << _contains_empty_geom << + // std::endl; + + auto count_greater_than_zero = + thrust::make_transform_iterator(count_begin, greater_than_zero_functor{}); + + { + thrust::device_vector count_greater_than_zero_truth( + count_greater_than_zero, count_greater_than_zero + _range.num_linestrings()); + + test::print_device_vector(count_greater_than_zero_truth, "count_greater_than_zero_truth: "); } - auto view() { - auto index_range = range( - _non_empty_geometry_prefix_sum.begin(), - _non_empty_geometry_prefix_sum.end() - ); - return segment_method_view{ - _range, - index_range, - _num_segments, - _contains_empty_geom - }; - } - -private: - ParentRange _range; - bool _contains_empty_geom; - index_t _num_segments; - rmm::device_uvector _non_empty_geometry_prefix_sum; + // Compute the number of empty linestrings + // auto key_begin = make_geometry_id_iterator(_range.geometry_offsets_begin(), + // _range.geometry_offsets_end()); + // { + // thrust::device_vector key_truth(key_begin, key_begin + _range.num_linestrings()); + + // test::print_device_vector(key_truth, "key_truth: "); + // } + + zero_data_async( + _non_empty_linestring_prefix_sum.begin(), _non_empty_linestring_prefix_sum.end(), stream); + // thrust::reduce_by_key(rmm::exec_policy(stream), + // key_begin, + // key_begin + _range.num_linestrings(), + // count_greater_than_zero, + // thrust::make_discard_iterator(), + // thrust::next(_non_empty_linestring_prefix_sum.begin()), + // thrust::equal_to{}, + // thrust::plus{}); + + thrust::inclusive_scan(rmm::exec_policy(stream), + count_greater_than_zero, + count_greater_than_zero + parent_range.num_linestrings(), + thrust::next(_non_empty_linestring_prefix_sum.begin())); + + _num_segments = _range.num_points() - _non_empty_linestring_prefix_sum.element( + _non_empty_linestring_prefix_sum.size() - 1, stream); + + test::print_device_vector(_non_empty_linestring_prefix_sum, "non_empty_geometry_prefix_sum: "); + } + + auto view() + { + auto index_range = + range(_non_empty_linestring_prefix_sum.begin(), _non_empty_linestring_prefix_sum.end()); + return segment_method_view{ + _range, index_range, _num_segments, _contains_empty_geom}; + } + + private: + ParentRange _range; + bool _contains_empty_geom; + index_t _num_segments; + rmm::device_uvector _non_empty_linestring_prefix_sum; }; } // namespace detail - - } // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/method/segment_method_view.cuh b/cpp/include/cuspatial/detail/method/segment_method_view.cuh index 0454b1ae4..0ffeec18c 100644 --- a/cpp/include/cuspatial/detail/method/segment_method_view.cuh +++ b/cpp/include/cuspatial/detail/method/segment_method_view.cuh @@ -16,8 +16,10 @@ #pragma once +#include "thrust/iterator/permutation_iterator.h" #include +#include #include #include #include @@ -45,32 +47,28 @@ class segment_method_view { _num_segments(num_segments), _contains_empty_ring(contains_empty_ring) { + thrust::device_vector soffsets(segment_offset_begin(), segment_offset_end()); + test::print_device_vector(soffsets, "Segment offsets: "); } - auto non_empty_linestring_count_begin() - { - auto begin = _non_empty_geometry_prefix_sum.begin(); - auto paired_begin = thrust::make_zip_iterator(begin, thrust::next(begin)); - return thrust::make_transform_iterator(paired_begin, offset_pair_to_count_functor{}); - } - - auto segment_count_begin() - { - auto num_points_begin = _range.multilinestring_point_count_begin(); - auto n_point_linestring_pair_it = - thrust::make_zip_iterator(num_points_begin, this->non_empty_linestring_count_begin()); + // CUSPATIAL_HOST_DEVICE auto segment_count_begin() + // { + // auto num_points_begin = _range.multilinestring_point_count_begin(); + // auto n_point_linestring_pair_it = + // thrust::make_zip_iterator(num_points_begin, this->_non_empty_linestring_count_begin()); - return thrust::make_transform_iterator(n_point_linestring_pair_it, - point_count_to_segment_count_functor{}); - } + // return thrust::make_transform_iterator(n_point_linestring_pair_it, + // point_count_to_segment_count_functor{}); + // } - auto segment_count_end() { - std::cout << "num multilinestrings: " << _range.num_multilinestrings() << std::endl; - return thrust::next(this->segment_count_begin(), _range.num_multilinestrings()); } + // CUSPATIAL_HOST_DEVICE auto segment_count_end() + // { + // return thrust::next(this->segment_count_begin(), _range.num_multilinestrings()); + // } - index_t num_segments() { return _num_segments; } + CUSPATIAL_HOST_DEVICE index_t num_segments() { return _num_segments; } - auto segment_offset_begin() + CUSPATIAL_HOST_DEVICE auto segment_offset_begin() { return make_counting_transform_iterator( 0, @@ -78,9 +76,28 @@ class segment_method_view { _non_empty_geometry_prefix_sum.begin()}); } - auto segment_offset_end() { return segment_offset_begin() + num_segments(); } + CUSPATIAL_HOST_DEVICE auto segment_offset_end() + { + return segment_offset_begin() + _non_empty_geometry_prefix_sum.size(); + } + + CUSPATIAL_HOST_DEVICE auto segment_count_begin() + { + auto permuted_offsets_it = + thrust::make_permutation_iterator(segment_offset_begin(), _range.geometry_offsets_begin()); + + auto zipped_offset_it = + thrust::make_zip_iterator(permuted_offsets_it, thrust::next(permuted_offsets_it)); + + return thrust::make_transform_iterator(zipped_offset_it, offset_pair_to_count_functor{}); + } + + CUSPATIAL_HOST_DEVICE auto segment_count_end() + { + return segment_count_begin() + _range.num_multilinestrings(); + } - auto segment_begin() + CUSPATIAL_HOST_DEVICE auto segment_begin() { return make_counting_transform_iterator( 0, @@ -90,13 +107,20 @@ class segment_method_view { _range.point_begin()}); } - auto segment_end() { return segment_begin() + _num_segments; } + CUSPATIAL_HOST_DEVICE auto segment_end() { return segment_begin() + _num_segments; } private: ParentRange _range; IndexRange _non_empty_geometry_prefix_sum; index_t _num_segments; bool _contains_empty_ring; + + CUSPATIAL_HOST_DEVICE auto _non_empty_linestring_count_begin() + { + auto begin = _non_empty_geometry_prefix_sum.begin(); + auto paired_begin = thrust::make_zip_iterator(begin, thrust::next(begin)); + return thrust::make_transform_iterator(paired_begin, offset_pair_to_count_functor{}); + } }; template diff --git a/cpp/tests/range/multilinestring_range_test.cu b/cpp/tests/range/multilinestring_range_test.cu index 78b1554a3..88d5dace4 100644 --- a/cpp/tests/range/multilinestring_range_test.cu +++ b/cpp/tests/range/multilinestring_range_test.cu @@ -13,6 +13,7 @@ * See the License for the specific language governing permissions and * limitations under the License. */ +#include #include #include @@ -49,25 +50,31 @@ struct MultilinestringRangeTest : public BaseFixture { } void run_segment_test_with_method_single(std::initializer_list geometry_offset, - std::initializer_list part_offset, - std::initializer_list> coordinates, - std::initializer_list> expected) + std::initializer_list part_offset, + std::initializer_list> coordinates, + std::initializer_list> expected) { auto multilinestring_array = make_multilinestring_array(geometry_offset, part_offset, coordinates); - auto rng = multilinestring_array.range(); + auto rng = multilinestring_array.range(); auto segment_methods = rng.segment_methods(stream()); - auto segment_view = segment_methods.view(); + auto segment_view = segment_methods.view(); + + auto segment_offsets = thrust::device_vector(segment_view.segment_offset_begin(), + segment_view.segment_offset_end()); + + test::print_device_vector(segment_offsets, "segment_offsets: "); rmm::device_uvector> got(segment_view.num_segments(), stream()); - thrust::copy(rmm::exec_policy(stream()), segment_view.segment_begin(), segment_view.segment_end(), got.begin()); + thrust::copy(rmm::exec_policy(stream()), + segment_view.segment_begin(), + segment_view.segment_end(), + got.begin()); auto d_expected = thrust::device_vector>(expected.begin(), expected.end()); CUSPATIAL_EXPECT_VEC2D_PAIRS_EQUIVALENT(d_expected, got); } - - void run_multilinestring_point_count_test(std::initializer_list geometry_offset, std::initializer_list part_offset, std::initializer_list> coordinates, @@ -89,16 +96,17 @@ struct MultilinestringRangeTest : public BaseFixture { CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); } - void run_multilinestring_segment_method_count_test(std::initializer_list geometry_offset, - std::initializer_list part_offset, - std::initializer_list> coordinates, - std::initializer_list expected) + void run_multilinestring_segment_method_count_test( + std::initializer_list geometry_offset, + std::initializer_list part_offset, + std::initializer_list> coordinates, + std::initializer_list expected) { auto multilinestring_array = make_multilinestring_array(geometry_offset, part_offset, coordinates); - auto rng = multilinestring_array.range(); + auto rng = multilinestring_array.range(); auto multilinestring_with_segment_methods = rng.segment_methods(stream()); - auto methods_view = multilinestring_with_segment_methods.view(); + auto methods_view = multilinestring_with_segment_methods.view(); auto d_expected = thrust::device_vector(expected.begin(), expected.end()); @@ -111,7 +119,6 @@ struct MultilinestringRangeTest : public BaseFixture { CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); } - void run_multilinestring_linestring_count_test(std::initializer_list geometry_offset, std::initializer_list part_offset, std::initializer_list> coordinates, @@ -173,8 +180,11 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorOneSegmentTest) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST( - this->run_segment_test_single, {0, 1}, {0, 2}, {P{0, 0}, P{1, 1}}, {S{P{0, 0}, P{1, 1}}}); + CUSPATIAL_RUN_TEST(this->run_segment_test_with_method_single, + {0, 1}, + {0, 2}, + {P{0, 0}, P{1, 1}}, + {S{P{0, 0}, P{1, 1}}}); } TYPED_TEST(MultilinestringRangeTest, SegmentIteratorTwoSegmentTest) @@ -183,7 +193,7 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorTwoSegmentTest) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_segment_test_single, + CUSPATIAL_RUN_TEST(this->run_segment_test_with_method_single, {0, 2}, {0, 2, 4}, {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}}, @@ -196,7 +206,7 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorTwoSegmentTest2) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_segment_test_single, + CUSPATIAL_RUN_TEST(this->run_segment_test_with_method_single, {0, 1, 2}, {0, 2, 4}, {P{0, 0}, P{1, 1}, P{0, 0}, P{1, 1}}, @@ -209,25 +219,29 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorManyPairTest) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_segment_test_single, + CUSPATIAL_RUN_TEST(this->run_segment_test_with_method_single, {0, 1, 2, 3}, {0, 6, 11, 14}, - {P{0, 0}, - P{1, 1}, - P{2, 2}, - P{3, 3}, - P{4, 4}, - P{5, 5}, - - P{10, 10}, - P{11, 11}, - P{12, 12}, - P{13, 13}, - P{14, 14}, - - P{20, 20}, - P{21, 21}, - P{22, 22}}, + { + + P{0, 0}, + P{1, 1}, + P{2, 2}, + P{3, 3}, + P{4, 4}, + P{5, 5}, + + P{10, 10}, + P{11, 11}, + P{12, 12}, + P{13, 13}, + P{14, 14}, + + P{20, 20}, + P{21, 21}, + P{22, 22} + + }, {S{P{0, 0}, P{1, 1}}, S{P{1, 1}, P{2, 2}}, @@ -288,8 +302,11 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST( - this->run_multilinestring_segment_method_count_test, {0, 1}, {0, 3}, {P{0, 0}, P{1, 1}, P{2, 2}}, {2}); + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, + {0, 1}, + {0, 3}, + {P{0, 0}, P{1, 1}, P{2, 2}}, + {2}); } TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest2) From 5e3a45aef2176f541bdda8ce2b4f2b99588187ce Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Thu, 18 May 2023 10:18:48 -0700 Subject: [PATCH 09/27] updates multipolygon tests to make use the new segment methods view --- cpp/tests/range/multipolygon_range_test.cu | 66 ++++++++-------------- 1 file changed, 24 insertions(+), 42 deletions(-) diff --git a/cpp/tests/range/multipolygon_range_test.cu b/cpp/tests/range/multipolygon_range_test.cu index ba76526d5..1862f4179 100644 --- a/cpp/tests/range/multipolygon_range_test.cu +++ b/cpp/tests/range/multipolygon_range_test.cu @@ -33,19 +33,25 @@ using namespace cuspatial::test; template struct MultipolygonRangeTest : public BaseFixture { - void run_multipolygon_segment_iterator_single(std::initializer_list geometry_offset, - std::initializer_list part_offset, - std::initializer_list ring_offset, - std::initializer_list> coordinates, - std::initializer_list> expected) + void run_multipolygon_segment_method_iterator_single( + std::initializer_list geometry_offset, + std::initializer_list part_offset, + std::initializer_list ring_offset, + std::initializer_list> coordinates, + std::initializer_list> expected) { auto multipolygon_array = make_multipolygon_array(geometry_offset, part_offset, ring_offset, coordinates); - auto rng = multipolygon_array.range(); + auto rng = multipolygon_array.range(); + auto segment_methods = rng.segment_methods(stream()); + auto segment_view = segment_methods.view(); - auto got = rmm::device_uvector>(rng.num_segments(), stream()); + auto got = rmm::device_uvector>(segment_view.num_segments(), stream()); - thrust::copy(rmm::exec_policy(stream()), rng.segment_begin(), rng.segment_end(), got.begin()); + thrust::copy(rmm::exec_policy(stream()), + segment_view.segment_begin(), + segment_view.segment_end(), + got.begin()); auto d_expected = thrust::device_vector>(expected.begin(), expected.end()); @@ -76,30 +82,6 @@ struct MultipolygonRangeTest : public BaseFixture { CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(got, d_expected); } - void run_multipolygon_segment_count_single( - std::initializer_list geometry_offset, - std::initializer_list part_offset, - std::initializer_list ring_offset, - std::initializer_list> coordinates, - std::initializer_list expected_segment_counts) - { - auto multipolygon_array = - make_multipolygon_array(geometry_offset, part_offset, ring_offset, coordinates); - auto rng = multipolygon_array.range(); - - auto got = rmm::device_uvector(rng.num_multipolygons(), stream()); - - thrust::copy(rmm::exec_policy(stream()), - rng.multipolygon_segment_count_begin(), - rng.multipolygon_segment_count_end(), - got.begin()); - - auto d_expected = thrust::device_vector(expected_segment_counts.begin(), - expected_segment_counts.end()); - - CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(got, d_expected); - } - void run_multipolygon_segment_method_count_single( std::initializer_list geometry_offset, std::initializer_list part_offset, @@ -109,9 +91,9 @@ struct MultipolygonRangeTest : public BaseFixture { { auto multipolygon_array = make_multipolygon_array(geometry_offset, part_offset, ring_offset, coordinates); - auto rng = multipolygon_array.range(); + auto rng = multipolygon_array.range(); auto segment_methods = rng.segment_methods(stream()); - auto segment_view = segment_methods.view(); + auto segment_view = segment_methods.view(); auto got = rmm::device_uvector(rng.num_multipolygons(), stream()); @@ -184,7 +166,7 @@ TYPED_TEST(MultipolygonRangeTest, SegmentIterators) using T = TypeParam; using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_iterator_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, {0, 1}, {0, 1}, {0, 4}, @@ -194,7 +176,7 @@ TYPED_TEST(MultipolygonRangeTest, SegmentIterators) TYPED_TEST(MultipolygonRangeTest, SegmentIterators2) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_iterator_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, {0, 1}, {0, 2}, {0, 4, 8}, @@ -209,7 +191,7 @@ TYPED_TEST(MultipolygonRangeTest, SegmentIterators2) TYPED_TEST(MultipolygonRangeTest, SegmentIterators3) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_iterator_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, {0, 2}, {0, 1, 2}, {0, 4, 8}, @@ -224,7 +206,7 @@ TYPED_TEST(MultipolygonRangeTest, SegmentIterators3) TYPED_TEST(MultipolygonRangeTest, SegmentIterators4) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_iterator_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, {0, 1, 2}, {0, 1, 2}, {0, 4, 8}, @@ -306,7 +288,7 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonCountIterator4) TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_count_single, {0, 1}, {0, 1}, {0, 4}, @@ -317,7 +299,7 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount) TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount2) { CUSPATIAL_RUN_TEST( - this->run_multipolygon_segment_count_single, + this->run_multipolygon_segment_method_count_single, {0, 1}, {0, 2}, {0, 4, 8}, @@ -327,7 +309,7 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount2) TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount3) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_count_single, {0, 2}, {0, 2, 3}, {0, 4, 8, 12}, @@ -348,7 +330,7 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount3) TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount4) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_count_single, {0, 2, 3}, {0, 2, 3, 4}, {0, 4, 8, 12, 16}, From 53e3ce03f25427e475b7a063b9cbd258b62ec06c Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Thu, 18 May 2023 10:32:13 -0700 Subject: [PATCH 10/27] remove segment methods from ranges and update all usage to `segment_methods` --- .../distance/linestring_polygon_distance.cuh | 56 ++++++++++++----- cpp/include/cuspatial/detail/functors.cuh | 20 ++++-- .../detail/range/multilinestring_range.cuh | 19 +----- .../detail/range/multipolygon_range.cuh | 62 +++---------------- .../cuspatial/range/multilinestring_range.cuh | 12 ++-- .../cuspatial/range/multipolygon_range.cuh | 15 +---- 6 files changed, 67 insertions(+), 117 deletions(-) diff --git a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh index aa0c7cc11..2ee28c74c 100644 --- a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh +++ b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh @@ -60,12 +60,16 @@ namespace detail { * @param distances Output range of distances, pre-filled with std::numerical_limits::max() */ template void __global__ pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestrings, + MultiLinestringSegmentRange multilinestring_segments, MultiPolygonRange multipolygons, + MultiPolygonSegmentRange multipolygon_segments, IndexRange thread_bounds, IndexRange multilinestrings_segment_offsets, IndexRange multipolygons_segment_offsets, @@ -88,9 +92,13 @@ pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestring continue; } + printf("geometry_id: %d\n", static_cast(geometry_id)); + printf("segment_count_range size: %d\n", + static_cast(thrust::distance(multilinestring_segments.segment_count_begin(), + multilinestring_segments.segment_count_end()))); // Retrieve the number of segments in multilinestrings[geometry_id] auto num_segment_this_multilinestring = - multilinestrings.multilinestring_segment_count_begin()[geometry_id]; + multilinestring_segments.segment_count_begin()[geometry_id]; // The segment id from the multilinestring this thread is computing (local_id + global_offset) auto multilinestring_segment_id = local_idx % num_segment_this_multilinestring + multilinestrings_segment_offsets[geometry_id]; @@ -98,8 +106,13 @@ pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestring auto multipolygon_segment_id = local_idx / num_segment_this_multilinestring + multipolygons_segment_offsets[geometry_id]; - auto [a, b] = multilinestrings.segment_begin()[multilinestring_segment_id]; - auto [c, d] = multipolygons.segment_begin()[multipolygon_segment_id]; + printf("multilinestring_segment_id: %d\n", static_cast(multilinestring_segment_id)); + printf("num segments: %d\n", + static_cast(thrust::distance(multilinestring_segments.segment_begin(), + multilinestring_segments.segment_end()))); + auto [a, b] = multilinestring_segments.segment_begin()[multilinestring_segment_id]; + printf("Here!\n"); + auto [c, d] = multipolygon_segments.segment_begin()[multipolygon_segment_id]; atomicMin(&distances[geometry_id], sqrt(squared_segment_distance(a, b, c, d))); } @@ -125,12 +138,22 @@ OutputIt pairwise_linestring_polygon_distance(MultiLinestringRange multilinestri auto multipoints = multilinestrings.as_multipoint_range(); auto multipoint_intersects = point_polygon_intersects(multipoints, multipolygons, stream); + // Make views to the segments in the multilinestring + auto multilinestring_segments = multilinestrings.segment_methods(stream); + auto multilinestring_segments_range = multilinestring_segments.view(); + auto multilinestring_segment_count_begin = multilinestring_segments_range.segment_count_begin(); + + // Make views to the segments in the multilinestring + auto multipolygon_segments = multipolygons.segment_methods(stream); + auto multipolygon_segments_range = multipolygon_segments.view(); + auto multipolygon_segment_count_begin = multipolygon_segments_range.segment_count_begin(); + // Compute the "boundary" of threads. Threads are partitioned based on the number of linestrings // times the number of polygons in a multipoint-multipolygon pair. - auto segment_count_product_it = thrust::make_transform_iterator( - thrust::make_zip_iterator(multilinestrings.multilinestring_segment_count_begin(), - multipolygons.multipolygon_segment_count_begin()), - thrust::make_zip_function(thrust::multiplies{})); + auto segment_count_product_it = + thrust::make_transform_iterator(thrust::make_zip_iterator(multilinestring_segment_count_begin, + multipolygon_segment_count_begin), + thrust::make_zip_function(thrust::multiplies{})); // Computes the "thread boundary" of each pair. This array partitions the thread range by // geometries. E.g. threadIdx within [thread_bounds[i], thread_bounds[i+1]) computes distances of @@ -154,17 +177,16 @@ OutputIt pairwise_linestring_polygon_distance(MultiLinestringRange multilinestri detail::zero_data_async( multipolygon_segment_offsets.begin(), multipolygon_segment_offsets.end(), stream); - thrust::inclusive_scan(rmm::exec_policy(stream), - multilinestrings.multilinestring_segment_count_begin(), - multilinestrings.multilinestring_segment_count_begin() + - multilinestrings.num_multilinestrings(), - thrust::next(multilinestring_segment_offsets.begin())); - thrust::inclusive_scan( rmm::exec_policy(stream), - multipolygons.multipolygon_segment_count_begin(), - multipolygons.multipolygon_segment_count_begin() + multipolygons.num_multipolygons(), - thrust::next(multipolygon_segment_offsets.begin())); + multilinestring_segment_count_begin, + multilinestring_segment_count_begin + multilinestrings.num_multilinestrings(), + thrust::next(multilinestring_segment_offsets.begin())); + + thrust::inclusive_scan(rmm::exec_policy(stream), + multipolygon_segment_count_begin, + multipolygon_segment_count_begin + multipolygons.num_multipolygons(), + thrust::next(multipolygon_segment_offsets.begin())); // Initialize output range thrust::fill(rmm::exec_policy(stream), @@ -177,7 +199,9 @@ OutputIt pairwise_linestring_polygon_distance(MultiLinestringRange multilinestri detail::pairwise_linestring_polygon_distance_kernel<<>>( multilinestrings, + multilinestring_segments_range, multipolygons, + multipolygon_segments_range, range{thread_bounds.begin(), thread_bounds.end()}, range{multilinestring_segment_offsets.begin(), multilinestring_segment_offsets.end()}, range{multipolygon_segment_offsets.begin(), multipolygon_segment_offsets.end()}, diff --git a/cpp/include/cuspatial/detail/functors.cuh b/cpp/include/cuspatial/detail/functors.cuh index e48a1082e..6a121a3fa 100644 --- a/cpp/include/cuspatial/detail/functors.cuh +++ b/cpp/include/cuspatial/detail/functors.cuh @@ -49,8 +49,8 @@ struct offset_pair_to_count_functor { /** * @brief Convert counts of points to counts of segments in a linestring. * - * A Multilinestring is composed of a series of Linestrings. Each Linestring is composed of a series of - * segments. The number of segments in a multilinestring is the number of points in the + * A Multilinestring is composed of a series of Linestrings. Each Linestring is composed of a series + * of segments. The number of segments in a multilinestring is the number of points in the * multilinestring minus the number of linestrings. * * Caveats: This has a strong assumption that the Multilinestring does not contain empty @@ -105,7 +105,8 @@ struct to_segment_offset_iterator { /// Deduction guide for to_distance_iterator template -to_segment_offset_iterator(OffsetIterator, CountIterator) -> to_segment_offset_iterator; +to_segment_offset_iterator(OffsetIterator, CountIterator) + -> to_segment_offset_iterator; /** * @brief Return a segment from the a partitioned range of points @@ -131,9 +132,16 @@ struct to_valid_segment_functor { template CUSPATIAL_HOST_DEVICE segment operator()(IndexType sid) { - auto kit = thrust::upper_bound(thrust::seq, segment_offset_begin, segment_offset_end, sid); - auto k = thrust::distance(segment_offset_begin, kit); - auto pid = non_empty_partitions_begin[sid] + k - 1; + auto kit = + thrust::prev(thrust::upper_bound(thrust::seq, segment_offset_begin, segment_offset_end, sid)); + auto geometry_id = thrust::distance(segment_offset_begin, kit); + auto preceding_non_empty_linestrings = non_empty_partitions_begin[geometry_id]; + auto pid = sid + preceding_non_empty_linestrings; + + printf("sid: %d geometry_id: %d pid: %d\n", + static_cast(sid), + static_cast(geometry_id), + static_cast(pid)); return segment{point_begin[pid], point_begin[pid + 1]}; } diff --git a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh index c0ac4d5e2..8dd834307 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh @@ -233,23 +233,6 @@ CUSPATIAL_HOST_DEVICE auto multilinestring_range -CUSPATIAL_HOST_DEVICE auto multilinestring_range:: - multilinestring_segment_count_begin() -{ - auto n_point_linestring_pair_it = thrust::make_zip_iterator( - multilinestring_point_count_begin(), multilinestring_linestring_count_begin()); - return thrust::make_transform_iterator(n_point_linestring_pair_it, - detail::point_count_to_segment_count_functor{}); -} - -template -CUSPATIAL_HOST_DEVICE auto multilinestring_range:: - multilinestring_segment_count_end() -{ - return multilinestring_segment_count_begin() + num_multilinestrings(); -} - template CUSPATIAL_HOST_DEVICE auto multilinestring_range:: multilinestring_linestring_count_begin() @@ -283,7 +266,7 @@ multilinestring_range::segment_end( } template -CUSPATIAL_HOST_DEVICE auto +auto multilinestring_range::segment_methods(rmm::cuda_stream_view stream) { return segment_method{*this, stream}; diff --git a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh index 6f1ae66eb..53e6c9043 100644 --- a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh @@ -324,32 +324,6 @@ multipolygon_range:: return multipolygon_ring_count_begin() + num_multipolygons(); } -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range:: - multipolygon_segment_count_begin() -{ - auto multipolygon_point_ring_count_it = - thrust::make_zip_iterator(multipolygon_point_count_begin(), multipolygon_ring_count_begin()); - - return thrust::make_transform_iterator(multipolygon_point_ring_count_it, - detail::point_count_to_segment_count_functor{}); -} - -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range:: - multipolygon_segment_count_end() -{ - return multipolygon_segment_count_begin() + num_multipolygons(); -} - template ::g return segment{_point_begin[segment_idx], _point_begin[segment_idx + 1]}; } -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range::segment_begin() -{ - return detail::make_counting_transform_iterator( - 0, - detail::to_valid_segment_functor{ - this->subtracted_ring_begin(), this->subtracted_ring_end(), thrust::make_counting_iterator(0), _point_begin}); -} - -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range::segment_end() -{ - return segment_begin() + num_segments(); -} - template :: subtracted_ring_begin() { - return detail::make_counting_transform_iterator(0, detail::to_segment_offset_iterator{ - _ring_begin, thrust::make_counting_iterator(0)}); + return detail::make_counting_transform_iterator( + 0, detail::to_segment_offset_iterator{_ring_begin, thrust::make_counting_iterator(0)}); } template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range::segment_methods(rmm::cuda_stream_view stream){ - auto multilinestring_range = this->as_multilinestring_range(); - return segment_method{multilinestring_range, stream}; +auto multipolygon_range::segment_methods( + rmm::cuda_stream_view stream) +{ + auto multilinestring_range = this->as_multilinestring_range(); + return segment_method{multilinestring_range, stream}; } template diff --git a/cpp/include/cuspatial/range/multipolygon_range.cuh b/cpp/include/cuspatial/range/multipolygon_range.cuh index 3619e8902..57e7589b4 100644 --- a/cpp/include/cuspatial/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/range/multipolygon_range.cuh @@ -16,7 +16,6 @@ #pragma once - #include #include #include @@ -27,7 +26,6 @@ #include - namespace cuspatial { /** @@ -185,18 +183,7 @@ class multipolygon_range { /// Returns the one past the iterator to the number of rings of the last multipolygon CUSPATIAL_HOST_DEVICE auto multipolygon_ring_count_end(); - /// Returns an iterator to the number of segments of the first multipolygon - CUSPATIAL_HOST_DEVICE auto multipolygon_segment_count_begin(); - /// Returns the one past the iterator to the number of segments of the last multipolygon - CUSPATIAL_HOST_DEVICE auto multipolygon_segment_count_end(); - - /// Returns an iterator to the start of the segment - CUSPATIAL_HOST_DEVICE auto segment_begin(); - - /// Returns an iterator to the end of the segment - CUSPATIAL_HOST_DEVICE auto segment_end(); - - CUSPATIAL_HOST_DEVICE auto segment_methods(rmm::cuda_stream_view); + auto segment_methods(rmm::cuda_stream_view); /// Range Casting From 789f7a23b82c5ebf215db05ff7ee45cc88b497e7 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Thu, 18 May 2023 17:44:55 +0000 Subject: [PATCH 11/27] Cleanup: Remove all debug prints --- .../detail/distance/linestring_polygon_distance.cuh | 9 --------- cpp/include/cuspatial/detail/functors.cuh | 5 ----- 2 files changed, 14 deletions(-) diff --git a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh index 2ee28c74c..eae5e0169 100644 --- a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh +++ b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh @@ -92,10 +92,6 @@ pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestring continue; } - printf("geometry_id: %d\n", static_cast(geometry_id)); - printf("segment_count_range size: %d\n", - static_cast(thrust::distance(multilinestring_segments.segment_count_begin(), - multilinestring_segments.segment_count_end()))); // Retrieve the number of segments in multilinestrings[geometry_id] auto num_segment_this_multilinestring = multilinestring_segments.segment_count_begin()[geometry_id]; @@ -106,12 +102,7 @@ pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestring auto multipolygon_segment_id = local_idx / num_segment_this_multilinestring + multipolygons_segment_offsets[geometry_id]; - printf("multilinestring_segment_id: %d\n", static_cast(multilinestring_segment_id)); - printf("num segments: %d\n", - static_cast(thrust::distance(multilinestring_segments.segment_begin(), - multilinestring_segments.segment_end()))); auto [a, b] = multilinestring_segments.segment_begin()[multilinestring_segment_id]; - printf("Here!\n"); auto [c, d] = multipolygon_segments.segment_begin()[multipolygon_segment_id]; atomicMin(&distances[geometry_id], sqrt(squared_segment_distance(a, b, c, d))); diff --git a/cpp/include/cuspatial/detail/functors.cuh b/cpp/include/cuspatial/detail/functors.cuh index 6a121a3fa..ce540dabf 100644 --- a/cpp/include/cuspatial/detail/functors.cuh +++ b/cpp/include/cuspatial/detail/functors.cuh @@ -138,11 +138,6 @@ struct to_valid_segment_functor { auto preceding_non_empty_linestrings = non_empty_partitions_begin[geometry_id]; auto pid = sid + preceding_non_empty_linestrings; - printf("sid: %d geometry_id: %d pid: %d\n", - static_cast(sid), - static_cast(geometry_id), - static_cast(pid)); - return segment{point_begin[pid], point_begin[pid + 1]}; } }; From 2390accb07921b06d5b8d16d772907ca0737e304 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Thu, 18 May 2023 17:47:58 +0000 Subject: [PATCH 12/27] Cleanup: Remove segment functions in multilinestring_range and dead codes --- .../detail/method/segment_method.cuh | 47 ++----------------- .../detail/method/segment_method_view.cuh | 26 ++-------- .../detail/range/multilinestring_range.cuh | 38 ++------------- .../cuspatial/range/multilinestring_range.cuh | 9 ---- 4 files changed, 9 insertions(+), 111 deletions(-) diff --git a/cpp/include/cuspatial/detail/method/segment_method.cuh b/cpp/include/cuspatial/detail/method/segment_method.cuh index e35c09b45..672db470c 100644 --- a/cpp/include/cuspatial/detail/method/segment_method.cuh +++ b/cpp/include/cuspatial/detail/method/segment_method.cuh @@ -39,16 +39,13 @@ namespace cuspatial { namespace detail { -template -struct equals { - __device__ bool operator()(IndexType x) const { return x == value; } -}; - template struct greater_than_zero_functor { __device__ IndexType operator()(IndexType x) const { return x > 0; } }; +// Optimization: for range that does not contain any empty linestrings, +// The _non_empty_linestring_prefix_sum can be initailized with counting_iterator. template class segment_method { using index_t = iterator_value_type; @@ -65,46 +62,11 @@ class segment_method { thrust::make_zip_iterator(offset_range.begin(), thrust::next(offset_range.begin())), offset_pair_to_count_functor{}); - // // Preemptive test: does the given range contain any empty ring/linestring? - // _contains_empty_geom = thrust::any_of( - // rmm::exec_policy(stream), - // count_begin, - // count_begin + _range.num_linestrings(), - // equals{} - // ); - - // std::cout << std::boolalpha << "contains empty geometry: " << _contains_empty_geom << - // std::endl; - auto count_greater_than_zero = thrust::make_transform_iterator(count_begin, greater_than_zero_functor{}); - { - thrust::device_vector count_greater_than_zero_truth( - count_greater_than_zero, count_greater_than_zero + _range.num_linestrings()); - - test::print_device_vector(count_greater_than_zero_truth, "count_greater_than_zero_truth: "); - } - - // Compute the number of empty linestrings - // auto key_begin = make_geometry_id_iterator(_range.geometry_offsets_begin(), - // _range.geometry_offsets_end()); - // { - // thrust::device_vector key_truth(key_begin, key_begin + _range.num_linestrings()); - - // test::print_device_vector(key_truth, "key_truth: "); - // } - zero_data_async( _non_empty_linestring_prefix_sum.begin(), _non_empty_linestring_prefix_sum.end(), stream); - // thrust::reduce_by_key(rmm::exec_policy(stream), - // key_begin, - // key_begin + _range.num_linestrings(), - // count_greater_than_zero, - // thrust::make_discard_iterator(), - // thrust::next(_non_empty_linestring_prefix_sum.begin()), - // thrust::equal_to{}, - // thrust::plus{}); thrust::inclusive_scan(rmm::exec_policy(stream), count_greater_than_zero, @@ -113,8 +75,6 @@ class segment_method { _num_segments = _range.num_points() - _non_empty_linestring_prefix_sum.element( _non_empty_linestring_prefix_sum.size() - 1, stream); - - test::print_device_vector(_non_empty_linestring_prefix_sum, "non_empty_geometry_prefix_sum: "); } auto view() @@ -122,12 +82,11 @@ class segment_method { auto index_range = range(_non_empty_linestring_prefix_sum.begin(), _non_empty_linestring_prefix_sum.end()); return segment_method_view{ - _range, index_range, _num_segments, _contains_empty_geom}; + _range, index_range, _num_segments}; } private: ParentRange _range; - bool _contains_empty_geom; index_t _num_segments; rmm::device_uvector _non_empty_linestring_prefix_sum; }; diff --git a/cpp/include/cuspatial/detail/method/segment_method_view.cuh b/cpp/include/cuspatial/detail/method/segment_method_view.cuh index 0ffeec18c..13a6655a5 100644 --- a/cpp/include/cuspatial/detail/method/segment_method_view.cuh +++ b/cpp/include/cuspatial/detail/method/segment_method_view.cuh @@ -16,7 +16,6 @@ #pragma once -#include "thrust/iterator/permutation_iterator.h" #include #include @@ -28,6 +27,7 @@ #include #include +#include #include namespace cuspatial { @@ -40,32 +40,13 @@ class segment_method_view { public: segment_method_view(ParentRange range, IndexRange non_empty_geometry_prefix_sum, - index_t num_segments, - bool contains_empty_ring) + index_t num_segments) : _range(range), _non_empty_geometry_prefix_sum(non_empty_geometry_prefix_sum), - _num_segments(num_segments), - _contains_empty_ring(contains_empty_ring) + _num_segments(num_segments) { - thrust::device_vector soffsets(segment_offset_begin(), segment_offset_end()); - test::print_device_vector(soffsets, "Segment offsets: "); } - // CUSPATIAL_HOST_DEVICE auto segment_count_begin() - // { - // auto num_points_begin = _range.multilinestring_point_count_begin(); - // auto n_point_linestring_pair_it = - // thrust::make_zip_iterator(num_points_begin, this->_non_empty_linestring_count_begin()); - - // return thrust::make_transform_iterator(n_point_linestring_pair_it, - // point_count_to_segment_count_functor{}); - // } - - // CUSPATIAL_HOST_DEVICE auto segment_count_end() - // { - // return thrust::next(this->segment_count_begin(), _range.num_multilinestrings()); - // } - CUSPATIAL_HOST_DEVICE index_t num_segments() { return _num_segments; } CUSPATIAL_HOST_DEVICE auto segment_offset_begin() @@ -113,7 +94,6 @@ class segment_method_view { ParentRange _range; IndexRange _non_empty_geometry_prefix_sum; index_t _num_segments; - bool _contains_empty_ring; CUSPATIAL_HOST_DEVICE auto _non_empty_linestring_count_begin() { diff --git a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh index 8dd834307..371e5a4e8 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh @@ -25,13 +25,13 @@ #include #include +#include #include #include #include #include #include #include -#include #include @@ -249,44 +249,12 @@ CUSPATIAL_HOST_DEVICE auto multilinestring_range -CUSPATIAL_HOST_DEVICE auto -multilinestring_range::segment_begin() -{ - return detail::make_counting_transform_iterator( - 0, - detail::to_valid_segment_functor{ - this->segment_offset_begin(), this->segment_offset_end(), thrust::make_counting_iterator(0), _point_begin}); -} - -template -CUSPATIAL_HOST_DEVICE auto -multilinestring_range::segment_end() -{ - return segment_begin() + num_segments(); -} - -template -auto -multilinestring_range::segment_methods(rmm::cuda_stream_view stream) +auto multilinestring_range::segment_methods( + rmm::cuda_stream_view stream) { return segment_method{*this, stream}; } -template -CUSPATIAL_HOST_DEVICE auto -multilinestring_range::segment_offset_begin() -{ - return detail::make_counting_transform_iterator(0, detail::to_segment_offset_iterator{ - _part_begin, thrust::make_counting_iterator(0)}); -} - -template -CUSPATIAL_HOST_DEVICE auto -multilinestring_range::segment_offset_end() -{ - return segment_offset_begin() + thrust::distance(_part_begin, _part_end); -} - template CUSPATIAL_HOST_DEVICE auto multilinestring_range::as_multipoint_range() diff --git a/cpp/include/cuspatial/range/multilinestring_range.cuh b/cpp/include/cuspatial/range/multilinestring_range.cuh index 119526947..93264df9d 100644 --- a/cpp/include/cuspatial/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/range/multilinestring_range.cuh @@ -162,12 +162,6 @@ class multilinestring_range { /// Returns an iterator to the counts of points per multilinestring CUSPATIAL_HOST_DEVICE auto multilinestring_linestring_count_end(); - /// Returns an iterator to the start of the segment - CUSPATIAL_HOST_DEVICE auto segment_begin(); - - /// Returns an iterator to the end of the segment - CUSPATIAL_HOST_DEVICE auto segment_end(); - /// Constructs a segment methods object, can only be constructed on host. /// To use segment methods on device, create a `segment_methods_view` /// See: `segment_methods::view` @@ -199,9 +193,6 @@ class multilinestring_range { VecIterator _point_begin; VecIterator _point_end; - CUSPATIAL_HOST_DEVICE auto segment_offset_begin(); - CUSPATIAL_HOST_DEVICE auto segment_offset_end(); - private: /// @internal /// Return the iterator to the part index where the point locates. From 728caa46f2d00c6f165eeaaf735842ce7630acc2 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Thu, 18 May 2023 10:58:10 -0700 Subject: [PATCH 13/27] remove num_segments --- .../cuspatial/detail/range/multilinestring_range.cuh | 7 ------- .../cuspatial/detail/range/multipolygon_range.cuh | 10 ---------- cpp/include/cuspatial/range/multilinestring_range.cuh | 3 --- cpp/include/cuspatial/range/multipolygon_range.cuh | 3 --- 4 files changed, 23 deletions(-) diff --git a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh index 371e5a4e8..e9984e803 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh @@ -117,13 +117,6 @@ multilinestring_range::num_points() return thrust::distance(_point_begin, _point_end); } -template -CUSPATIAL_HOST_DEVICE auto -multilinestring_range::num_segments() -{ - return num_points() - num_linestrings(); -} - template CUSPATIAL_HOST_DEVICE auto multilinestring_range::multilinestring_begin() diff --git a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh index 53e6c9043..ece47478d 100644 --- a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh @@ -156,16 +156,6 @@ multipolygon_range::n return thrust::distance(_point_begin, _point_end); } -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range::num_segments() -{ - return num_points() - num_rings(); -} - template Date: Thu, 18 May 2023 13:01:50 -0700 Subject: [PATCH 14/27] rename segment_methods -> multilinestring_segment --- .../distance/linestring_polygon_distance.cuh | 15 +++--- ...method.cuh => multilinestring_segment.cuh} | 31 ++++++------ .../detail/range/multilinestring_range.cuh | 6 +-- .../multilinestring_segment_range.cuh} | 50 ++++++++----------- .../detail/range/multipolygon_range.cuh | 8 +-- .../cuspatial/range/multilinestring_range.cuh | 7 ++- .../cuspatial/range/multipolygon_range.cuh | 2 +- cpp/tests/range/multilinestring_range_test.cu | 29 ++++------- cpp/tests/range/multipolygon_range_test.cu | 24 ++++----- 9 files changed, 76 insertions(+), 96 deletions(-) rename cpp/include/cuspatial/detail/{method/segment_method.cuh => multilinestring_segment.cuh} (69%) rename cpp/include/cuspatial/detail/{method/segment_method_view.cuh => range/multilinestring_segment_range.cuh} (58%) diff --git a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh index eae5e0169..027680c4f 100644 --- a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh +++ b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh @@ -93,8 +93,7 @@ pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestring } // Retrieve the number of segments in multilinestrings[geometry_id] - auto num_segment_this_multilinestring = - multilinestring_segments.segment_count_begin()[geometry_id]; + auto num_segment_this_multilinestring = multilinestring_segments.count_begin()[geometry_id]; // The segment id from the multilinestring this thread is computing (local_id + global_offset) auto multilinestring_segment_id = local_idx % num_segment_this_multilinestring + multilinestrings_segment_offsets[geometry_id]; @@ -102,8 +101,8 @@ pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestring auto multipolygon_segment_id = local_idx / num_segment_this_multilinestring + multipolygons_segment_offsets[geometry_id]; - auto [a, b] = multilinestring_segments.segment_begin()[multilinestring_segment_id]; - auto [c, d] = multipolygon_segments.segment_begin()[multipolygon_segment_id]; + auto [a, b] = multilinestring_segments.begin()[multilinestring_segment_id]; + auto [c, d] = multipolygon_segments.begin()[multipolygon_segment_id]; atomicMin(&distances[geometry_id], sqrt(squared_segment_distance(a, b, c, d))); } @@ -130,14 +129,14 @@ OutputIt pairwise_linestring_polygon_distance(MultiLinestringRange multilinestri auto multipoint_intersects = point_polygon_intersects(multipoints, multipolygons, stream); // Make views to the segments in the multilinestring - auto multilinestring_segments = multilinestrings.segment_methods(stream); + auto multilinestring_segments = multilinestrings._segments(stream); auto multilinestring_segments_range = multilinestring_segments.view(); - auto multilinestring_segment_count_begin = multilinestring_segments_range.segment_count_begin(); + auto multilinestring_segment_count_begin = multilinestring_segments_range.count_begin(); // Make views to the segments in the multilinestring - auto multipolygon_segments = multipolygons.segment_methods(stream); + auto multipolygon_segments = multipolygons._segments(stream); auto multipolygon_segments_range = multipolygon_segments.view(); - auto multipolygon_segment_count_begin = multipolygon_segments_range.segment_count_begin(); + auto multipolygon_segment_count_begin = multipolygon_segments_range.count_begin(); // Compute the "boundary" of threads. Threads are partitioned based on the number of linestrings // times the number of polygons in a multipoint-multipolygon pair. diff --git a/cpp/include/cuspatial/detail/method/segment_method.cuh b/cpp/include/cuspatial/detail/multilinestring_segment.cuh similarity index 69% rename from cpp/include/cuspatial/detail/method/segment_method.cuh rename to cpp/include/cuspatial/detail/multilinestring_segment.cuh index 672db470c..68a29f1c5 100644 --- a/cpp/include/cuspatial/detail/method/segment_method.cuh +++ b/cpp/include/cuspatial/detail/multilinestring_segment.cuh @@ -18,7 +18,7 @@ #include -#include "segment_method_view.cuh" +#include "range/multilinestring_segment_range.cuh" #include #include @@ -46,18 +46,17 @@ struct greater_than_zero_functor { // Optimization: for range that does not contain any empty linestrings, // The _non_empty_linestring_prefix_sum can be initailized with counting_iterator. -template -class segment_method { - using index_t = iterator_value_type; +template +class multilinestring_segment { + using index_t = iterator_value_type; public: // segment_methods is always internal use, thus memory consumed is always temporary, // therefore always use default device memory resource. - segment_method(ParentRange parent_range, rmm::cuda_stream_view stream) - : _range(parent_range), - _non_empty_linestring_prefix_sum(parent_range.num_linestrings() + 1, stream) + multilinestring_segment(MultilinestringRange parent, rmm::cuda_stream_view stream) + : _parent(parent), _non_empty_linestring_prefix_sum(parent.num_linestrings() + 1, stream) { - auto offset_range = range(_range.part_offset_begin(), _range.part_offset_end()); + auto offset_range = ::cuspatial::range{_parent.part_offset_begin(), _parent.part_offset_end()}; auto count_begin = thrust::make_transform_iterator( thrust::make_zip_iterator(offset_range.begin(), thrust::next(offset_range.begin())), offset_pair_to_count_functor{}); @@ -70,23 +69,23 @@ class segment_method { thrust::inclusive_scan(rmm::exec_policy(stream), count_greater_than_zero, - count_greater_than_zero + parent_range.num_linestrings(), + count_greater_than_zero + _parent.num_linestrings(), thrust::next(_non_empty_linestring_prefix_sum.begin())); - _num_segments = _range.num_points() - _non_empty_linestring_prefix_sum.element( - _non_empty_linestring_prefix_sum.size() - 1, stream); + _num_segments = _parent.num_points() - _non_empty_linestring_prefix_sum.element( + _non_empty_linestring_prefix_sum.size() - 1, stream); } auto view() { - auto index_range = - range(_non_empty_linestring_prefix_sum.begin(), _non_empty_linestring_prefix_sum.end()); - return segment_method_view{ - _range, index_range, _num_segments}; + auto index_range = ::cuspatial::range{_non_empty_linestring_prefix_sum.begin(), + _non_empty_linestring_prefix_sum.end()}; + return multilinestring_segment_range{ + _parent, index_range, _num_segments}; } private: - ParentRange _range; + MultilinestringRange _parent; index_t _num_segments; rmm::device_uvector _non_empty_linestring_prefix_sum; }; diff --git a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh index e9984e803..f2dc48b36 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh @@ -25,7 +25,7 @@ #include #include -#include +#include #include #include #include @@ -242,10 +242,10 @@ CUSPATIAL_HOST_DEVICE auto multilinestring_range -auto multilinestring_range::segment_methods( +auto multilinestring_range::_segments( rmm::cuda_stream_view stream) { - return segment_method{*this, stream}; + return multilinestring_segment{*this, stream}; } template diff --git a/cpp/include/cuspatial/detail/method/segment_method_view.cuh b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh similarity index 58% rename from cpp/include/cuspatial/detail/method/segment_method_view.cuh rename to cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh index 13a6655a5..089b17169 100644 --- a/cpp/include/cuspatial/detail/method/segment_method_view.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh @@ -34,14 +34,14 @@ namespace cuspatial { namespace detail { template -class segment_method_view { +class multilinestring_segment_range { using index_t = typename IndexRange::value_type; public: - segment_method_view(ParentRange range, - IndexRange non_empty_geometry_prefix_sum, - index_t num_segments) - : _range(range), + multilinestring_segment_range(ParentRange parent, + IndexRange non_empty_geometry_prefix_sum, + index_t num_segments) + : _parent(parent), _non_empty_geometry_prefix_sum(non_empty_geometry_prefix_sum), _num_segments(num_segments) { @@ -49,23 +49,23 @@ class segment_method_view { CUSPATIAL_HOST_DEVICE index_t num_segments() { return _num_segments; } - CUSPATIAL_HOST_DEVICE auto segment_offset_begin() + CUSPATIAL_HOST_DEVICE auto offset_begin() { return make_counting_transform_iterator( 0, - to_segment_offset_iterator{_range.part_offset_begin(), + to_segment_offset_iterator{_parent.part_offset_begin(), _non_empty_geometry_prefix_sum.begin()}); } - CUSPATIAL_HOST_DEVICE auto segment_offset_end() + CUSPATIAL_HOST_DEVICE auto offset_end() { - return segment_offset_begin() + _non_empty_geometry_prefix_sum.size(); + return offset_begin() + _non_empty_geometry_prefix_sum.size(); } - CUSPATIAL_HOST_DEVICE auto segment_count_begin() + CUSPATIAL_HOST_DEVICE auto count_begin() { auto permuted_offsets_it = - thrust::make_permutation_iterator(segment_offset_begin(), _range.geometry_offsets_begin()); + thrust::make_permutation_iterator(offset_begin(), _parent.geometry_offsets_begin()); auto zipped_offset_it = thrust::make_zip_iterator(permuted_offsets_it, thrust::next(permuted_offsets_it)); @@ -73,39 +73,29 @@ class segment_method_view { return thrust::make_transform_iterator(zipped_offset_it, offset_pair_to_count_functor{}); } - CUSPATIAL_HOST_DEVICE auto segment_count_end() - { - return segment_count_begin() + _range.num_multilinestrings(); - } + CUSPATIAL_HOST_DEVICE auto count_end() { return count_begin() + _parent.num_multilinestrings(); } - CUSPATIAL_HOST_DEVICE auto segment_begin() + CUSPATIAL_HOST_DEVICE auto begin() { return make_counting_transform_iterator( 0, - to_valid_segment_functor{segment_offset_begin(), - segment_offset_end(), + to_valid_segment_functor{offset_begin(), + offset_end(), _non_empty_geometry_prefix_sum.begin(), - _range.point_begin()}); + _parent.point_begin()}); } - CUSPATIAL_HOST_DEVICE auto segment_end() { return segment_begin() + _num_segments; } + CUSPATIAL_HOST_DEVICE auto end() { return begin() + _num_segments; } private: - ParentRange _range; + ParentRange _parent; IndexRange _non_empty_geometry_prefix_sum; index_t _num_segments; - - CUSPATIAL_HOST_DEVICE auto _non_empty_linestring_count_begin() - { - auto begin = _non_empty_geometry_prefix_sum.begin(); - auto paired_begin = thrust::make_zip_iterator(begin, thrust::next(begin)); - return thrust::make_transform_iterator(paired_begin, offset_pair_to_count_functor{}); - } }; template -segment_method_view(ParentRange, IndexRange, typename IndexRange::value_type, bool) - -> segment_method_view; +multilinestring_segment_range(ParentRange, IndexRange, typename IndexRange::value_type, bool) + -> multilinestring_segment_range; } // namespace detail diff --git a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh index ece47478d..c6f9e75b7 100644 --- a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh @@ -16,9 +16,9 @@ #pragma once -#include "thrust/iterator/counting_iterator.h" #include #include +#include #include #include #include @@ -29,8 +29,10 @@ #include #include + #include #include +#include #include #include #include @@ -379,11 +381,11 @@ template -auto multipolygon_range::segment_methods( +auto multipolygon_range::_segments( rmm::cuda_stream_view stream) { auto multilinestring_range = this->as_multilinestring_range(); - return segment_method{multilinestring_range, stream}; + return multilinestring_segment{multilinestring_range, stream}; } template diff --git a/cpp/include/cuspatial/range/multipolygon_range.cuh b/cpp/include/cuspatial/range/multipolygon_range.cuh index 3ad310f04..d36c83633 100644 --- a/cpp/include/cuspatial/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/range/multipolygon_range.cuh @@ -180,7 +180,7 @@ class multipolygon_range { /// Returns the one past the iterator to the number of rings of the last multipolygon CUSPATIAL_HOST_DEVICE auto multipolygon_ring_count_end(); - auto segment_methods(rmm::cuda_stream_view); + auto _segments(rmm::cuda_stream_view); /// Range Casting diff --git a/cpp/tests/range/multilinestring_range_test.cu b/cpp/tests/range/multilinestring_range_test.cu index 88d5dace4..71f823979 100644 --- a/cpp/tests/range/multilinestring_range_test.cu +++ b/cpp/tests/range/multilinestring_range_test.cu @@ -56,20 +56,13 @@ struct MultilinestringRangeTest : public BaseFixture { { auto multilinestring_array = make_multilinestring_array(geometry_offset, part_offset, coordinates); - auto rng = multilinestring_array.range(); - auto segment_methods = rng.segment_methods(stream()); - auto segment_view = segment_methods.view(); + auto rng = multilinestring_array.range(); + auto segments = rng._segments(stream()); + auto segments_range = segments.view(); - auto segment_offsets = thrust::device_vector(segment_view.segment_offset_begin(), - segment_view.segment_offset_end()); - - test::print_device_vector(segment_offsets, "segment_offsets: "); - - rmm::device_uvector> got(segment_view.num_segments(), stream()); - thrust::copy(rmm::exec_policy(stream()), - segment_view.segment_begin(), - segment_view.segment_end(), - got.begin()); + rmm::device_uvector> got(segments_range.num_segments(), stream()); + thrust::copy( + rmm::exec_policy(stream()), segments_range.begin(), segments_range.end(), got.begin()); auto d_expected = thrust::device_vector>(expected.begin(), expected.end()); CUSPATIAL_EXPECT_VEC2D_PAIRS_EQUIVALENT(d_expected, got); @@ -104,16 +97,16 @@ struct MultilinestringRangeTest : public BaseFixture { { auto multilinestring_array = make_multilinestring_array(geometry_offset, part_offset, coordinates); - auto rng = multilinestring_array.range(); - auto multilinestring_with_segment_methods = rng.segment_methods(stream()); - auto methods_view = multilinestring_with_segment_methods.view(); + auto rng = multilinestring_array.range(); + auto segments = rng._segments(stream()); + auto segments_range = segments.view(); auto d_expected = thrust::device_vector(expected.begin(), expected.end()); rmm::device_uvector got(rng.num_multilinestrings(), stream()); thrust::copy(rmm::exec_policy(stream()), - methods_view.segment_count_begin(), - methods_view.segment_count_end(), + segments_range.count_begin(), + segments_range.count_end(), got.begin()); CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); diff --git a/cpp/tests/range/multipolygon_range_test.cu b/cpp/tests/range/multipolygon_range_test.cu index 1862f4179..550891115 100644 --- a/cpp/tests/range/multipolygon_range_test.cu +++ b/cpp/tests/range/multipolygon_range_test.cu @@ -42,16 +42,14 @@ struct MultipolygonRangeTest : public BaseFixture { { auto multipolygon_array = make_multipolygon_array(geometry_offset, part_offset, ring_offset, coordinates); - auto rng = multipolygon_array.range(); - auto segment_methods = rng.segment_methods(stream()); - auto segment_view = segment_methods.view(); + auto rng = multipolygon_array.range(); + auto segments = rng._segments(stream()); + auto segment_range = segments.view(); - auto got = rmm::device_uvector>(segment_view.num_segments(), stream()); + auto got = rmm::device_uvector>(segment_range.num_segments(), stream()); - thrust::copy(rmm::exec_policy(stream()), - segment_view.segment_begin(), - segment_view.segment_end(), - got.begin()); + thrust::copy( + rmm::exec_policy(stream()), segment_range.begin(), segment_range.end(), got.begin()); auto d_expected = thrust::device_vector>(expected.begin(), expected.end()); @@ -91,15 +89,15 @@ struct MultipolygonRangeTest : public BaseFixture { { auto multipolygon_array = make_multipolygon_array(geometry_offset, part_offset, ring_offset, coordinates); - auto rng = multipolygon_array.range(); - auto segment_methods = rng.segment_methods(stream()); - auto segment_view = segment_methods.view(); + auto rng = multipolygon_array.range(); + auto segments = rng._segments(stream()); + auto segment_range = segments.view(); auto got = rmm::device_uvector(rng.num_multipolygons(), stream()); thrust::copy(rmm::exec_policy(stream()), - segment_view.segment_count_begin(), - segment_view.segment_count_end(), + segment_range.count_begin(), + segment_range.count_end(), got.begin()); auto d_expected = thrust::device_vector(expected_segment_counts.begin(), From e03044971fe85426a6c8ac502373754f287a775f Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Fri, 19 May 2023 16:18:36 +0000 Subject: [PATCH 15/27] Updates segment_range and tests to make sure APIs exposed are all multi-geometry based --- .../detail/multilinestring_segment.cuh | 5 +- .../range/multilinestring_segment_range.cuh | 45 +++++--- cpp/tests/range/multilinestring_range_test.cu | 103 +++++++++++++++++- cpp/tests/range/multipolygon_range_test.cu | 87 ++++++++++++++- 4 files changed, 212 insertions(+), 28 deletions(-) diff --git a/cpp/include/cuspatial/detail/multilinestring_segment.cuh b/cpp/include/cuspatial/detail/multilinestring_segment.cuh index 68a29f1c5..3bf0b0195 100644 --- a/cpp/include/cuspatial/detail/multilinestring_segment.cuh +++ b/cpp/include/cuspatial/detail/multilinestring_segment.cuh @@ -16,11 +16,8 @@ #pragma once -#include - -#include "range/multilinestring_segment_range.cuh" - #include +#include #include #include #include diff --git a/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh index 089b17169..40b042a20 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh @@ -16,8 +16,6 @@ #pragma once -#include - #include #include #include @@ -49,38 +47,36 @@ class multilinestring_segment_range { CUSPATIAL_HOST_DEVICE index_t num_segments() { return _num_segments; } - CUSPATIAL_HOST_DEVICE auto offset_begin() + CUSPATIAL_HOST_DEVICE auto multigeometry_offset_begin() { - return make_counting_transform_iterator( - 0, - to_segment_offset_iterator{_parent.part_offset_begin(), - _non_empty_geometry_prefix_sum.begin()}); + return thrust::make_permutation_iterator(_per_linestring_offset_begin(), + _parent.geometry_offsets_begin()); } - CUSPATIAL_HOST_DEVICE auto offset_end() + CUSPATIAL_HOST_DEVICE auto multigeometry_offset_end() { - return offset_begin() + _non_empty_geometry_prefix_sum.size(); + return multigeometry_offset_begin() + _parent.num_multilinestrings() + 1; } - CUSPATIAL_HOST_DEVICE auto count_begin() + CUSPATIAL_HOST_DEVICE auto multigeometry_count_begin() { - auto permuted_offsets_it = - thrust::make_permutation_iterator(offset_begin(), _parent.geometry_offsets_begin()); - - auto zipped_offset_it = - thrust::make_zip_iterator(permuted_offsets_it, thrust::next(permuted_offsets_it)); + auto zipped_offset_it = thrust::make_zip_iterator(multigeometry_offset_begin(), + thrust::next(multigeometry_offset_begin())); return thrust::make_transform_iterator(zipped_offset_it, offset_pair_to_count_functor{}); } - CUSPATIAL_HOST_DEVICE auto count_end() { return count_begin() + _parent.num_multilinestrings(); } + CUSPATIAL_HOST_DEVICE auto multigeometry_count_end() + { + return multigeometry_count_begin() + _parent.num_multilinestrings(); + } CUSPATIAL_HOST_DEVICE auto begin() { return make_counting_transform_iterator( 0, - to_valid_segment_functor{offset_begin(), - offset_end(), + to_valid_segment_functor{_per_linestring_offset_begin(), + _per_linestring_offset_end(), _non_empty_geometry_prefix_sum.begin(), _parent.point_begin()}); } @@ -91,6 +87,19 @@ class multilinestring_segment_range { ParentRange _parent; IndexRange _non_empty_geometry_prefix_sum; index_t _num_segments; + + CUSPATIAL_HOST_DEVICE auto _per_linestring_offset_begin() + { + return make_counting_transform_iterator( + 0, + to_segment_offset_iterator{_parent.part_offset_begin(), + _non_empty_geometry_prefix_sum.begin()}); + } + + CUSPATIAL_HOST_DEVICE auto _per_linestring_offset_end() + { + return _per_linestring_offset_begin() + _non_empty_geometry_prefix_sum.size(); + } }; template diff --git a/cpp/tests/range/multilinestring_range_test.cu b/cpp/tests/range/multilinestring_range_test.cu index 71f823979..9e617dd7b 100644 --- a/cpp/tests/range/multilinestring_range_test.cu +++ b/cpp/tests/range/multilinestring_range_test.cu @@ -105,8 +105,31 @@ struct MultilinestringRangeTest : public BaseFixture { rmm::device_uvector got(rng.num_multilinestrings(), stream()); thrust::copy(rmm::exec_policy(stream()), - segments_range.count_begin(), - segments_range.count_end(), + segments_range.multigeometry_count_begin(), + segments_range.multigeometry_count_end(), + got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); + } + + void run_multilinestring_segment_offset_test(std::initializer_list geometry_offset, + std::initializer_list part_offset, + std::initializer_list> coordinates, + std::initializer_list expected) + { + auto multilinestring_array = + make_multilinestring_array(geometry_offset, part_offset, coordinates); + auto rng = multilinestring_array.range(); + auto segments = rng._segments(stream()); + auto segments_range = segments.view(); + + auto d_expected = thrust::device_vector(expected.begin(), expected.end()); + + rmm::device_uvector got(rng.num_multilinestrings() + 1, stream()); + + thrust::copy(rmm::exec_policy(stream()), + segments_range.multigeometry_offset_begin(), + segments_range.multigeometry_offset_end(), got.begin()); CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); @@ -263,6 +286,38 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorWithEmptyLineTest) {S{P{0, 0}, P{1, 1}}, S{P{1, 1}, P{2, 2}}, S{P{10, 10}, P{11, 11}}, S{P{11, 11}, P{12, 12}}}); } +TYPED_TEST(MultilinestringRangeTest, SegmentIteratorWithEmptyMultiLineStringTest) +{ + using T = TypeParam; + using P = vec_2d; + using S = segment; + + CUSPATIAL_RUN_TEST( + this->run_segment_test_with_method_single, + {0, 1, 1, 2}, + {0, 3, 6}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, + {S{P{0, 0}, P{1, 1}}, S{P{1, 1}, P{2, 2}}, S{P{10, 10}, P{11, 11}}, S{P{11, 11}, P{12, 12}}}); +} + +TYPED_TEST(MultilinestringRangeTest, SegmentIteratorWithEmptyMultiLineStringTest2) +{ + using T = TypeParam; + using P = vec_2d; + using S = segment; + + CUSPATIAL_RUN_TEST(this->run_segment_test_with_method_single, + + {0, 1, 1, 2}, + {0, 4, 7}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}, P{10, 10}, P{11, 11}, P{12, 12}}, + {S{P{0, 0}, P{1, 1}}, + S{P{1, 1}, P{2, 2}}, + S{P{2, 2}, P{3, 3}}, + S{P{10, 10}, P{11, 11}}, + S{P{11, 11}, P{12, 12}}}); +} + TYPED_TEST(MultilinestringRangeTest, PerMultilinestringCountTest) { using T = TypeParam; @@ -341,7 +396,7 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest4) {2, 2}); } -// FIXME: contains empty linestring +// contains empty linestring TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest5) { using T = TypeParam; @@ -355,6 +410,48 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest5) {2, 0, 2}); } +// contains empty multilinestring +TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest6) +{ + using T = TypeParam; + using P = vec_2d; + using S = segment; + + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, + {0, 1, 1, 2}, + {0, 3, 6}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, + {2, 0, 2}); +} + +// contains empty multilinestring +TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest7) +{ + using T = TypeParam; + using P = vec_2d; + using S = segment; + + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, + {0, 1, 1, 2}, + {0, 4, 7}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}, P{10, 10}, P{11, 11}, P{12, 12}}, + {3, 0, 2}); +} + +// contains empty multilinestring +TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentOffsetTest) +{ + using T = TypeParam; + using P = vec_2d; + using S = segment; + + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_offset_test, + {0, 1, 1, 2}, + {0, 4, 7}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}, P{10, 10}, P{11, 11}, P{12, 12}}, + {0, 3, 3, 5}); +} + TYPED_TEST(MultilinestringRangeTest, MultilinestringLinestringCountTest) { using T = TypeParam; diff --git a/cpp/tests/range/multipolygon_range_test.cu b/cpp/tests/range/multipolygon_range_test.cu index 550891115..e737f033b 100644 --- a/cpp/tests/range/multipolygon_range_test.cu +++ b/cpp/tests/range/multipolygon_range_test.cu @@ -96,8 +96,8 @@ struct MultipolygonRangeTest : public BaseFixture { auto got = rmm::device_uvector(rng.num_multipolygons(), stream()); thrust::copy(rmm::exec_policy(stream()), - segment_range.count_begin(), - segment_range.count_end(), + segment_range.multigeometry_count_begin(), + segment_range.multigeometry_count_end(), got.begin()); auto d_expected = thrust::device_vector(expected_segment_counts.begin(), @@ -169,7 +169,7 @@ TYPED_TEST(MultipolygonRangeTest, SegmentIterators) {0, 1}, {0, 4}, {{0, 0}, {1, 0}, {1, 1}, {0, 0}}, - {S{P{0, 0}, P{1, 0}}, S{P{1, 0}, P{1, 1}}, S{P{1, 1}, P{0, 0}}}); + {S{{0, 0}, P{1, 0}}, S{P{1, 0}, P{1, 1}}, S{P{1, 1}, P{0, 0}}}); } TYPED_TEST(MultipolygonRangeTest, SegmentIterators2) @@ -217,6 +217,87 @@ TYPED_TEST(MultipolygonRangeTest, SegmentIterators4) {{11, 11}, {10, 10}}}); } +TYPED_TEST(MultipolygonRangeTest, SegmentIterators5) +{ + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, + {0, 1, 2, 3}, + {0, 1, 2, 3}, + {0, 4, 9, 14}, + {{-1, -1}, + {-2, -2}, + {-2, -1}, + {-1, -1}, + + {-20, -20}, + {-20, -21}, + {-21, -21}, + {-21, -20}, + {-20, -20}, + + {-10, -10}, + {-10, -11}, + {-11, -11}, + {-11, -10}, + {-10, -10}}, + + {{{-1, -1}, {-2, -2}}, + {{-2, -2}, {-2, -1}}, + {{-2, -1}, {-1, -1}}, + {{-20, -20}, {-20, -21}}, + {{-20, -21}, {-21, -21}}, + {{-21, -21}, {-21, -20}}, + {{-21, -20}, {-20, -20}}, + {{-10, -10}, {-10, -11}}, + {{-10, -11}, {-11, -11}}, + {{-11, -11}, {-11, -10}}, + {{-11, -10}, {-10, -10}}}); +} + +TYPED_TEST(MultipolygonRangeTest, SegmentIterators5EmptyRing) +{ + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, + {0, 1, 2}, + {0, 1, 3}, + {0, 4, 4, 8}, + {{0, 0}, {1, 0}, {1, 1}, {0, 0}, {10, 10}, {11, 10}, {11, 11}, {10, 10}}, + {{{0, 0}, {1, 0}}, + {{1, 0}, {1, 1}}, + {{1, 1}, {0, 0}}, + {{10, 10}, {11, 10}}, + {{11, 10}, {11, 11}}, + {{11, 11}, {10, 10}}}); +} + +TYPED_TEST(MultipolygonRangeTest, SegmentIterators6EmptyPolygon) +{ + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, + {0, 1, 3}, + {0, 1, 1, 2}, + {0, 4, 8}, + {{0, 0}, {1, 0}, {1, 1}, {0, 0}, {10, 10}, {11, 10}, {11, 11}, {10, 10}}, + {{{0, 0}, {1, 0}}, + {{1, 0}, {1, 1}}, + {{1, 1}, {0, 0}}, + {{10, 10}, {11, 10}}, + {{11, 10}, {11, 11}}, + {{11, 11}, {10, 10}}}); +} + +TYPED_TEST(MultipolygonRangeTest, SegmentIterators7EmptyMultiPolygon) +{ + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, + {0, 1, 1, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{0, 0}, {1, 0}, {1, 1}, {0, 0}, {10, 10}, {11, 10}, {11, 11}, {10, 10}}, + {{{0, 0}, {1, 0}}, + {{1, 0}, {1, 1}}, + {{1, 1}, {0, 0}}, + {{10, 10}, {11, 10}}, + {{11, 10}, {11, 11}}, + {{11, 11}, {10, 10}}}); +} + TYPED_TEST(MultipolygonRangeTest, MultipolygonCountIterator) { CUSPATIAL_RUN_TEST(this->run_multipolygon_point_count_iterator_single, From 771f340c2dc76217c2f7da894360a1c11e8fd551 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Fri, 19 May 2023 16:19:50 +0000 Subject: [PATCH 16/27] Refactors `linestring_polygon_distance` and add tests to include empty multigeometry test case. --- .../distance/linestring_polygon_distance.cuh | 120 ++++++++---------- .../multilinestring_ref.cuh | 3 + .../geometry_collection/multipolygon_ref.cuh | 3 + .../cuspatial_test/vector_equality.hpp | 12 +- .../linestring_polygon_distance_test.cu | 105 +++++++++++++++ 5 files changed, 173 insertions(+), 70 deletions(-) diff --git a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh index 027680c4f..1c5cb26d8 100644 --- a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh +++ b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh @@ -19,7 +19,6 @@ #include "distance_utils.cuh" #include -#include #include #include #include @@ -30,10 +29,10 @@ #include #include -#include -#include +#include #include #include +#include #include #include @@ -45,40 +44,44 @@ namespace cuspatial { namespace detail { +/// Device functor that returns true if any of the geometry is empty. +struct any_input_is_empty { + template + bool __device__ operator()(LhsType lhs, RhsType rhs) + { + return lhs.is_empty() || rhs.is_empty(); + } +}; + /** * @brief Computes distances between the multilinestring and multipolygons * - * @param multilinestrings Range to the multilinestring - * @param multipolygons Range to the multipolygon + * This is a load balanced distance compute kernel. Each thread compute exactly 1 pair of segments + * between the multilinestring and multipolygon. + * + * @tparam T type of the underlying coordinates + * @tparam index_t type of underlying offsets + * + * @param multilinestring_segments Range to the segments of the multilinestring + * @param multipolygon_segments Range to the segments of the multipolygon * @param thread_bounds Range to the boundary of thread partitions - * @param multilinestrings_segment_offsets Range to the indices where the first segment of each - * multilinestring begins - * @param multipolygons_segment_offsets Range to the indices where the first segment of each - * multipolygon begins * @param intersects A uint8_t array that indicates if the corresponding pair of multipoint and * multipolygon intersects * @param distances Output range of distances, pre-filled with std::numerical_limits::max() */ -template void __global__ -pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestrings, - MultiLinestringSegmentRange multilinestring_segments, - MultiPolygonRange multipolygons, +pairwise_linestring_polygon_distance_kernel(MultiLinestringSegmentRange multilinestring_segments, MultiPolygonSegmentRange multipolygon_segments, IndexRange thread_bounds, - IndexRange multilinestrings_segment_offsets, - IndexRange multipolygons_segment_offsets, uint8_t* intersects, OutputIt* distances) { - using T = typename MultiLinestringRange::element_t; - using index_t = iterator_value_type; - auto num_threads = thread_bounds[thread_bounds.size() - 1]; for (auto idx = blockDim.x * blockIdx.x + threadIdx.x; idx < num_threads; idx += blockDim.x * gridDim.x) { @@ -93,13 +96,15 @@ pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestring } // Retrieve the number of segments in multilinestrings[geometry_id] - auto num_segment_this_multilinestring = multilinestring_segments.count_begin()[geometry_id]; + auto num_segment_this_multilinestring = + multilinestring_segments.multigeometry_count_begin()[geometry_id]; // The segment id from the multilinestring this thread is computing (local_id + global_offset) auto multilinestring_segment_id = - local_idx % num_segment_this_multilinestring + multilinestrings_segment_offsets[geometry_id]; + local_idx % num_segment_this_multilinestring + + multilinestring_segments.multigeometry_offset_begin()[geometry_id]; // The segment id from the multipolygon this thread is computing (local_id + global_offset) - auto multipolygon_segment_id = - local_idx / num_segment_this_multilinestring + multipolygons_segment_offsets[geometry_id]; + auto multipolygon_segment_id = local_idx / num_segment_this_multilinestring + + multipolygon_segments.multigeometry_offset_begin()[geometry_id]; auto [a, b] = multilinestring_segments.begin()[multilinestring_segment_id]; auto [c, d] = multipolygon_segments.begin()[multipolygon_segment_id]; @@ -122,24 +127,27 @@ OutputIt pairwise_linestring_polygon_distance(MultiLinestringRange multilinestri CUSPATIAL_EXPECTS(multilinestrings.size() == multipolygons.size(), "Must have the same number of input rows."); - if (multilinestrings.size() == 0) return distances_first; + auto size = multilinestrings.size(); + + if (size == 0) return distances_first; // Create a multipoint range from multilinestrings, computes intersection auto multipoints = multilinestrings.as_multipoint_range(); auto multipoint_intersects = point_polygon_intersects(multipoints, multipolygons, stream); // Make views to the segments in the multilinestring - auto multilinestring_segments = multilinestrings._segments(stream); - auto multilinestring_segments_range = multilinestring_segments.view(); - auto multilinestring_segment_count_begin = multilinestring_segments_range.count_begin(); + auto multilinestring_segments = multilinestrings._segments(stream); + auto multilinestring_segments_range = multilinestring_segments.view(); + auto multilinestring_segment_count_begin = + multilinestring_segments_range.multigeometry_count_begin(); // Make views to the segments in the multilinestring auto multipolygon_segments = multipolygons._segments(stream); auto multipolygon_segments_range = multipolygon_segments.view(); - auto multipolygon_segment_count_begin = multipolygon_segments_range.count_begin(); + auto multipolygon_segment_count_begin = multipolygon_segments_range.multigeometry_count_begin(); // Compute the "boundary" of threads. Threads are partitioned based on the number of linestrings - // times the number of polygons in a multipoint-multipolygon pair. + // times the number of polygons in a multilinestring-multipolygon pair. auto segment_count_product_it = thrust::make_transform_iterator(thrust::make_zip_iterator(multilinestring_segment_count_begin, multipolygon_segment_count_begin), @@ -156,49 +164,33 @@ OutputIt pairwise_linestring_polygon_distance(MultiLinestringRange multilinestri segment_count_product_it + thread_bounds.size() - 1, thrust::next(thread_bounds.begin())); - // Compute offsets to the first segment of each multilinestring and multipolygon - auto multilinestring_segment_offsets = - rmm::device_uvector(multilinestrings.num_multilinestrings() + 1, stream); - detail::zero_data_async( - multilinestring_segment_offsets.begin(), multilinestring_segment_offsets.end(), stream); - - auto multipolygon_segment_offsets = - rmm::device_uvector(multipolygons.num_multipolygons() + 1, stream); - detail::zero_data_async( - multipolygon_segment_offsets.begin(), multipolygon_segment_offsets.end(), stream); - - thrust::inclusive_scan( - rmm::exec_policy(stream), - multilinestring_segment_count_begin, - multilinestring_segment_count_begin + multilinestrings.num_multilinestrings(), - thrust::next(multilinestring_segment_offsets.begin())); - - thrust::inclusive_scan(rmm::exec_policy(stream), - multipolygon_segment_count_begin, - multipolygon_segment_count_begin + multipolygons.num_multipolygons(), - thrust::next(multipolygon_segment_offsets.begin())); - // Initialize output range thrust::fill(rmm::exec_policy(stream), distances_first, - distances_first + multilinestrings.num_multilinestrings(), + distances_first + size, std::numeric_limits::max()); + // If any input multigeometries is empty, result is nan. + auto nan_it = thrust::make_constant_iterator(std::numeric_limits::quiet_NaN()); + thrust::transform_if(rmm::exec_policy(stream), + nan_it, + nan_it + size, + thrust::make_zip_iterator(multilinestrings.begin(), multipolygons.begin()), + distances_first, + thrust::identity{}, + thrust::make_zip_function(detail::any_input_is_empty{})); + auto num_threads = thread_bounds.back_element(stream); auto [tpb, num_blocks] = grid_1d(num_threads); - detail::pairwise_linestring_polygon_distance_kernel<<>>( - multilinestrings, - multilinestring_segments_range, - multipolygons, - multipolygon_segments_range, - range{thread_bounds.begin(), thread_bounds.end()}, - range{multilinestring_segment_offsets.begin(), multilinestring_segment_offsets.end()}, - range{multipolygon_segment_offsets.begin(), multipolygon_segment_offsets.end()}, - multipoint_intersects.begin(), - distances_first); - - return distances_first + multilinestrings.num_multilinestrings(); + detail::pairwise_linestring_polygon_distance_kernel + <<>>(multilinestring_segments_range, + multipolygon_segments_range, + range{thread_bounds.begin(), thread_bounds.end()}, + multipoint_intersects.begin(), + distances_first); + + return distances_first + size; } } // namespace cuspatial diff --git a/cpp/include/cuspatial/geometry_collection/multilinestring_ref.cuh b/cpp/include/cuspatial/geometry_collection/multilinestring_ref.cuh index 872135f72..36304a7bc 100644 --- a/cpp/include/cuspatial/geometry_collection/multilinestring_ref.cuh +++ b/cpp/include/cuspatial/geometry_collection/multilinestring_ref.cuh @@ -38,6 +38,9 @@ class multilinestring_ref { /// Return the number of linestrings in the multilinestring. CUSPATIAL_HOST_DEVICE auto size() const { return num_linestrings(); } + /// Return true if this multilinestring contains 0 linestrings. + CUSPATIAL_HOST_DEVICE bool is_empty() const { return num_linestrings() == 0; } + /// Return iterator to the first linestring. CUSPATIAL_HOST_DEVICE auto part_begin() const; /// Return iterator to one past the last linestring. diff --git a/cpp/include/cuspatial/geometry_collection/multipolygon_ref.cuh b/cpp/include/cuspatial/geometry_collection/multipolygon_ref.cuh index e4ff5010f..d9972d9b1 100644 --- a/cpp/include/cuspatial/geometry_collection/multipolygon_ref.cuh +++ b/cpp/include/cuspatial/geometry_collection/multipolygon_ref.cuh @@ -41,6 +41,9 @@ class multipolygon_ref { /// Return the number of polygons in the multipolygon. CUSPATIAL_HOST_DEVICE auto size() const { return num_polygons(); } + /// Returns true if the multipolygon contains 0 geometries. + CUSPATIAL_HOST_DEVICE bool is_empty() const { return num_polygons() == 0; } + /// Return iterator to the first polygon. CUSPATIAL_HOST_DEVICE auto part_begin() const; /// Return iterator to one past the last polygon. diff --git a/cpp/include/cuspatial_test/vector_equality.hpp b/cpp/include/cuspatial_test/vector_equality.hpp index 5d72f6228..4abde035d 100644 --- a/cpp/include/cuspatial_test/vector_equality.hpp +++ b/cpp/include/cuspatial_test/vector_equality.hpp @@ -37,7 +37,7 @@ namespace cuspatial { namespace test { /** - * @brief Compare two floats are close within N ULPs + * @brief Compare two floats are close within N ULPs, nans are treated equal * * N is predefined by GoogleTest * https://google.github.io/googletest/reference/assertions.html#EXPECT_FLOAT_EQ @@ -46,22 +46,22 @@ template auto floating_eq_by_ulp(T val) { if constexpr (std::is_same_v) { - return ::testing::FloatEq(val); + return ::testing::NanSensitiveFloatEq(val); } else { - return ::testing::DoubleEq(val); + return ::testing::NanSensitiveDoubleEq(val); } } /** - * @brief Compare two floats are close within `abs_error` + * @brief Compare two floats are close within `abs_error`, nans are treated equal */ template auto floating_eq_by_abs_error(T val, T abs_error) { if constexpr (std::is_same_v) { - return ::testing::FloatNear(val, abs_error); + return ::testing::NanSensitiveFloatNear(val, abs_error); } else { - return ::testing::DoubleNear(val, abs_error); + return ::testing::NanSensitiveDoubleNear(val, abs_error); } } diff --git a/cpp/tests/distance/linestring_polygon_distance_test.cu b/cpp/tests/distance/linestring_polygon_distance_test.cu index bf74b731b..7cc3b6b55 100644 --- a/cpp/tests/distance/linestring_polygon_distance_test.cu +++ b/cpp/tests/distance/linestring_polygon_distance_test.cu @@ -413,3 +413,108 @@ TYPED_TEST(PairwiseLinestringPolygonDistanceTest, TwoPairsCrosses) P{5, 5}}, {std::sqrt(T{2}), 0.0}); } + +// Empty Geometries Tests + +/// Empty MultiLinestring vs Non-empty multipolygons +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, ThreePairEmptyMultiLinestring) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + + {0, 1, 1, 2}, + {0, 4, 7}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}, P{10, 10}, P{11, 11}, P{12, 12}}, + + {0, 1, 2, 3}, + {0, 1, 2, 3}, + {0, 4, 9, 14}, + {P{-1, -1}, + P{-2, -2}, + P{-2, -1}, + P{-1, -1}, + + P{-20, -20}, + P{-20, -21}, + P{-21, -21}, + P{-21, -20}, + P{-20, -20}, + + P{-10, -10}, + P{-10, -11}, + P{-11, -11}, + P{-11, -10}, + P{-10, -10}}, + + {std::sqrt(T{2}), std::numeric_limits::quiet_NaN(), 20 * std::sqrt(T{2})}); +} + +/// Non-empty MultiLinestring vs Empty multipolygons +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, ThreePairEmptyMultiPolygon) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + + {0, 1, 2, 3}, + {0, 4, 7, 10}, + {P{0, 0}, + P{1, 1}, + P{2, 2}, + P{3, 3}, + P{20, 20}, + P{21, 21}, + P{22, 22}, + P{10, 10}, + P{11, 11}, + P{12, 12}}, + + {0, 1, 1, 2}, + {0, 1, 2}, + {0, 4, 9}, + {P{-1, -1}, + P{-2, -2}, + P{-2, -1}, + P{-1, -1}, + + P{-10, -10}, + P{-10, -11}, + P{-11, -11}, + P{-11, -10}, + P{-10, -10}}, + {std::sqrt(T{2}), std::numeric_limits::quiet_NaN(), 20 * std::sqrt(T{2})}); +} + +/// FIXME: Empty MultiLinestring vs Empty multipolygons +/// This example fails at distance util, where point-polyogn intersection kernel doesn't handle +/// empty multipoint/multipolygons. +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, + DISABLED_ThreePairEmptyMultiLineStringEmptyMultiPolygon) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + + {0, 1, 1, 3}, + {0, 4, 7}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}, P{10, 10}, P{11, 11}, P{12, 12}}, + + {0, 1, 1, 2}, + {0, 1, 2, 3}, + {0, 4, 9, 14}, + {P{-1, -1}, + P{-2, -2}, + P{-2, -1}, + P{-1, -1}, + + P{-10, -10}, + P{-10, -11}, + P{-11, -11}, + P{-11, -10}, + P{-10, -10}}, + {std::sqrt(T{2}), std::numeric_limits::quiet_NaN(), 20 * std::sqrt(T{2})}); +} From 81e05bd1a6c0e2e7447cd7a40eac5d121fe4a4e3 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Fri, 19 May 2023 16:27:15 +0000 Subject: [PATCH 17/27] Move all segment range only functors to local file --- cpp/include/cuspatial/detail/functors.cuh | 101 ------------------ .../detail/range/multilinestring_range.cuh | 13 ++- .../range/multilinestring_segment_range.cuh | 101 ++++++++++++++++++ 3 files changed, 107 insertions(+), 108 deletions(-) diff --git a/cpp/include/cuspatial/detail/functors.cuh b/cpp/include/cuspatial/detail/functors.cuh index ce540dabf..25c79885d 100644 --- a/cpp/include/cuspatial/detail/functors.cuh +++ b/cpp/include/cuspatial/detail/functors.cuh @@ -46,106 +46,5 @@ struct offset_pair_to_count_functor { } }; -/** - * @brief Convert counts of points to counts of segments in a linestring. - * - * A Multilinestring is composed of a series of Linestrings. Each Linestring is composed of a series - * of segments. The number of segments in a multilinestring is the number of points in the - * multilinestring minus the number of linestrings. - * - * Caveats: This has a strong assumption that the Multilinestring does not contain empty - * linestrings. While each non-empty linestring in the multilinestring represents 1 extra segment, - * an empty multilinestring does not introduce any extra segments since it does not contain any - * points. - * - * Used to create segment count iterators, such as `multi*_segment_count_begin`. - * - * @tparam IndexPair Must be iterator to a pair of counts - * @param n_point_linestring_pair A pair of counts, the first is the number of points, the second is - * the number of linestrings. - */ -struct point_count_to_segment_count_functor { - template - CUSPATIAL_HOST_DEVICE auto operator()(IndexPair n_point_linestring_pair) - { - auto nPoints = thrust::get<0>(n_point_linestring_pair); - auto nLinestrings = thrust::get<1>(n_point_linestring_pair); - return nPoints - nLinestrings; - } -}; - -/** - * @brief Given an offset iterator it that partitions a point range, return an offset iterator that - * partitions the segment range made from the same point range. - * - * One partition to a point range introduces one invalid segment, except empty partitions. - * Therefore, the offsets that partitions the segment range is the offset that partitions the point - * range subtracts the number of *non-empty* point partitions that precedes the current point range. - * - * @tparam OffsetIterator Iterator type to the offset - * - * Caveats: This has a strong assumption that the Multilinestring does not contain empty - * linestrings. While each non-empty linestring in the multilinestring represents 1 extra segment, - * an empty multilinestring does not introduce any extra segments since it does not contain any - * points. - * - * Used to create iterator to segment offsets, such as `segment_offset_begin`. - */ -template -struct to_segment_offset_iterator { - OffsetIterator point_partition_begin; - CountIterator non_empty_partitions_begin; - - template - CUSPATIAL_HOST_DEVICE auto operator()(IndexType i) - { - return point_partition_begin[i] - non_empty_partitions_begin[i]; - } -}; - -/// Deduction guide for to_distance_iterator -template -to_segment_offset_iterator(OffsetIterator, CountIterator) - -> to_segment_offset_iterator; - -/** - * @brief Return a segment from the a partitioned range of points - * - * Used in a counting transform iterator. Given an index of the segment, offset it by the number of - * skipped segments preceding i in the partitioned range of points. Dereference the corresponding - * point and the point following to make a segment. - * - * Used to create iterator to segments, such as `segment_begin`. - * - * @tparam OffsetIterator the iterator type indicating partitions of the point range. - * @tparam CoordinateIterator the iterator type to the point range. - */ -template -struct to_valid_segment_functor { - using element_t = iterator_vec_base_type; - - OffsetIterator segment_offset_begin; - OffsetIterator segment_offset_end; - CountIterator non_empty_partitions_begin; - CoordinateIterator point_begin; - - template - CUSPATIAL_HOST_DEVICE segment operator()(IndexType sid) - { - auto kit = - thrust::prev(thrust::upper_bound(thrust::seq, segment_offset_begin, segment_offset_end, sid)); - auto geometry_id = thrust::distance(segment_offset_begin, kit); - auto preceding_non_empty_linestrings = non_empty_partitions_begin[geometry_id]; - auto pid = sid + preceding_non_empty_linestrings; - - return segment{point_begin[pid], point_begin[pid + 1]}; - } -}; - -/// Deduction guide for to_valid_segment_functor -template -to_valid_segment_functor(OffsetIterator, OffsetIterator, CountIterator, CoordinateIterator) - -> to_valid_segment_functor; - } // namespace detail } // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh index f2dc48b36..244e622e7 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh @@ -16,13 +16,6 @@ #pragma once -#include "thrust/iterator/counting_iterator.h" -#include -#include -#include -#include -#include - #include #include #include @@ -33,7 +26,13 @@ #include #include +#include +#include +#include #include +#include +#include +#include #include #include diff --git a/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh index 40b042a20..d71e148b1 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh @@ -31,6 +31,107 @@ namespace cuspatial { namespace detail { +/** + * @brief Convert counts of points to counts of segments in a linestring. + * + * A Multilinestring is composed of a series of Linestrings. Each Linestring is composed of a series + * of segments. The number of segments in a multilinestring is the number of points in the + * multilinestring minus the number of linestrings. + * + * Caveats: This has a strong assumption that the Multilinestring does not contain empty + * linestrings. While each non-empty linestring in the multilinestring represents 1 extra segment, + * an empty multilinestring does not introduce any extra segments since it does not contain any + * points. + * + * Used to create segment count iterators, such as `multi*_segment_count_begin`. + * + * @tparam IndexPair Must be iterator to a pair of counts + * @param n_point_linestring_pair A pair of counts, the first is the number of points, the second is + * the number of linestrings. + */ +struct point_count_to_segment_count_functor { + template + CUSPATIAL_HOST_DEVICE auto operator()(IndexPair n_point_linestring_pair) + { + auto nPoints = thrust::get<0>(n_point_linestring_pair); + auto nLinestrings = thrust::get<1>(n_point_linestring_pair); + return nPoints - nLinestrings; + } +}; + +/** + * @brief Given an offset iterator it that partitions a point range, return an offset iterator that + * partitions the segment range made from the same point range. + * + * One partition to a point range introduces one invalid segment, except empty partitions. + * Therefore, the offsets that partitions the segment range is the offset that partitions the point + * range subtracts the number of *non-empty* point partitions that precedes the current point range. + * + * @tparam OffsetIterator Iterator type to the offset + * + * Caveats: This has a strong assumption that the Multilinestring does not contain empty + * linestrings. While each non-empty linestring in the multilinestring represents 1 extra segment, + * an empty multilinestring does not introduce any extra segments since it does not contain any + * points. + * + * Used to create iterator to segment offsets, such as `segment_offset_begin`. + */ +template +struct to_segment_offset_iterator { + OffsetIterator point_partition_begin; + CountIterator non_empty_partitions_begin; + + template + CUSPATIAL_HOST_DEVICE auto operator()(IndexType i) + { + return point_partition_begin[i] - non_empty_partitions_begin[i]; + } +}; + +/// Deduction guide for to_distance_iterator +template +to_segment_offset_iterator(OffsetIterator, CountIterator) + -> to_segment_offset_iterator; + +/** + * @brief Return a segment from the a partitioned range of points + * + * Used in a counting transform iterator. Given an index of the segment, offset it by the number of + * skipped segments preceding i in the partitioned range of points. Dereference the corresponding + * point and the point following to make a segment. + * + * Used to create iterator to segments, such as `segment_begin`. + * + * @tparam OffsetIterator the iterator type indicating partitions of the point range. + * @tparam CoordinateIterator the iterator type to the point range. + */ +template +struct to_valid_segment_functor { + using element_t = iterator_vec_base_type; + + OffsetIterator segment_offset_begin; + OffsetIterator segment_offset_end; + CountIterator non_empty_partitions_begin; + CoordinateIterator point_begin; + + template + CUSPATIAL_HOST_DEVICE segment operator()(IndexType sid) + { + auto kit = + thrust::prev(thrust::upper_bound(thrust::seq, segment_offset_begin, segment_offset_end, sid)); + auto geometry_id = thrust::distance(segment_offset_begin, kit); + auto preceding_non_empty_linestrings = non_empty_partitions_begin[geometry_id]; + auto pid = sid + preceding_non_empty_linestrings; + + return segment{point_begin[pid], point_begin[pid + 1]}; + } +}; + +/// Deduction guide for to_valid_segment_functor +template +to_valid_segment_functor(OffsetIterator, OffsetIterator, CountIterator, CoordinateIterator) + -> to_valid_segment_functor; + template class multilinestring_segment_range { using index_t = typename IndexRange::value_type; From 9d00b4158267de14295f9d2cf4c8273727c184c7 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Fri, 19 May 2023 17:35:43 +0000 Subject: [PATCH 18/27] Documentation for segment_range --- .../detail/multilinestring_segment.cuh | 44 ++++- .../range/multilinestring_segment_range.cuh | 172 ++++++++++++------ 2 files changed, 155 insertions(+), 61 deletions(-) diff --git a/cpp/include/cuspatial/detail/multilinestring_segment.cuh b/cpp/include/cuspatial/detail/multilinestring_segment.cuh index 3bf0b0195..0f6215117 100644 --- a/cpp/include/cuspatial/detail/multilinestring_segment.cuh +++ b/cpp/include/cuspatial/detail/multilinestring_segment.cuh @@ -36,20 +36,51 @@ namespace cuspatial { namespace detail { +/** + * @internal + * @brief Functor that returns true if current value is greater than 0. + */ template struct greater_than_zero_functor { __device__ IndexType operator()(IndexType x) const { return x > 0; } }; -// Optimization: for range that does not contain any empty linestrings, -// The _non_empty_linestring_prefix_sum can be initailized with counting_iterator. +/** + * @internal + * @brief Owning class to provide iterators to segments in a multilinestring range + * + * The owned memory in this struct is `_non_empty_linestring_prefix_sum`, which equals the + * number of linestrings in the multilinestring range plus 1. This vector holds the number + * of non empty linestrings that precedes the current linestring. + * + * This class is only meant for tracking the life time of the owned memory. To access the + * segment iterators, call `view()` function to create a non-owning object of this class. + * + * For detailed explanation on the implementation of the segment iterators, see documentation + * of `multilinestring_segment_range`. + * + * @note To use this class with a multipolygon range, cast the multipolygon range as a + * multilinestring range. + * + * TODO: Optimization: for range that does not contain any empty linestrings, + * `_non_empty_linestring_prefix_sum` can be substituted with `counting_iterator`. + * + * @tparam MultilinestringRange The multilinestring range to initialize this class with. + */ template class multilinestring_segment { using index_t = iterator_value_type; public: - // segment_methods is always internal use, thus memory consumed is always temporary, - // therefore always use default device memory resource. + /** + * @brief Construct a new multilinestring segment object + * + * @note multilinestring_segment is always internal use, thus memory consumed is always + * temporary, therefore always use default device memory resource. + * + * @param parent The parent multilinestring object to construct from + * @param stream The stream to perform computation on + */ multilinestring_segment(MultilinestringRange parent, rmm::cuda_stream_view stream) : _parent(parent), _non_empty_linestring_prefix_sum(parent.num_linestrings() + 1, stream) { @@ -73,6 +104,11 @@ class multilinestring_segment { _non_empty_linestring_prefix_sum.size() - 1, stream); } + /** + * @brief Return a non-owning `multilinestring_segment_range` object from this class + * + * @return multilinestring_segment_range + */ auto view() { auto index_range = ::cuspatial::range{_non_empty_linestring_prefix_sum.begin(), diff --git a/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh index d71e148b1..44a87fec7 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh @@ -32,78 +32,64 @@ namespace cuspatial { namespace detail { /** - * @brief Convert counts of points to counts of segments in a linestring. + * @brief Computes the offsets to the starting segment per linestring * - * A Multilinestring is composed of a series of Linestrings. Each Linestring is composed of a series - * of segments. The number of segments in a multilinestring is the number of points in the - * multilinestring minus the number of linestrings. + * The point indices and segment indices are correlated, but in a different index space. + * For example: + * ``` + * {0, 3} + * {0, 3, 3, 6} + * {A, B, C, X, Y, Z} + * ``` * - * Caveats: This has a strong assumption that the Multilinestring does not contain empty - * linestrings. While each non-empty linestring in the multilinestring represents 1 extra segment, - * an empty multilinestring does not introduce any extra segments since it does not contain any - * points. + * ``` + * segments: AB BC XY YZ + * sid: 0 1 2 3 + * points: A B C X Y Z + * pid: 0 1 2 3 4 5 + * ``` * - * Used to create segment count iterators, such as `multi*_segment_count_begin`. + * The original {0, 3, 3, 6} offsets are in the point index space. For example: + * The first and past the last point of the first linestring is at point index 0 and 3 (A, X). + * The first and past the last point of the second linestring is at point index 3 and 3 (empty), + * and so on. * - * @tparam IndexPair Must be iterator to a pair of counts - * @param n_point_linestring_pair A pair of counts, the first is the number of points, the second is - * the number of linestrings. - */ -struct point_count_to_segment_count_functor { - template - CUSPATIAL_HOST_DEVICE auto operator()(IndexPair n_point_linestring_pair) - { - auto nPoints = thrust::get<0>(n_point_linestring_pair); - auto nLinestrings = thrust::get<1>(n_point_linestring_pair); - return nPoints - nLinestrings; - } -}; - -/** - * @brief Given an offset iterator it that partitions a point range, return an offset iterator that - * partitions the segment range made from the same point range. - * - * One partition to a point range introduces one invalid segment, except empty partitions. - * Therefore, the offsets that partitions the segment range is the offset that partitions the point - * range subtracts the number of *non-empty* point partitions that precedes the current point range. + * The transformed segment offsets {0, 2, 2, 4} are in the segment index space. For example: + * The first and past the last segment of the first linestring is at segment index 0 and 2 ((AB), + * (XY)). + * The first and past the last segment of the second linestring is at segment index 2 and 2 + * (empty), and so on. * * @tparam OffsetIterator Iterator type to the offset - * - * Caveats: This has a strong assumption that the Multilinestring does not contain empty - * linestrings. While each non-empty linestring in the multilinestring represents 1 extra segment, - * an empty multilinestring does not introduce any extra segments since it does not contain any - * points. - * - * Used to create iterator to segment offsets, such as `segment_offset_begin`. */ template -struct to_segment_offset_iterator { - OffsetIterator point_partition_begin; - CountIterator non_empty_partitions_begin; +struct point_offset_to_segment_offset { + OffsetIterator part_offset_begin; + CountIterator non_empty_linestrings_count_begin; template CUSPATIAL_HOST_DEVICE auto operator()(IndexType i) { - return point_partition_begin[i] - non_empty_partitions_begin[i]; + return part_offset_begin[i] - non_empty_linestrings_count_begin[i]; } }; -/// Deduction guide for to_distance_iterator +/// Deduction guide template -to_segment_offset_iterator(OffsetIterator, CountIterator) - -> to_segment_offset_iterator; +point_offset_to_segment_offset(OffsetIterator, CountIterator) + -> point_offset_to_segment_offset; /** - * @brief Return a segment from the a partitioned range of points - * - * Used in a counting transform iterator. Given an index of the segment, offset it by the number of - * skipped segments preceding i in the partitioned range of points. Dereference the corresponding - * point and the point following to make a segment. + * @brief Given a segment index, return the corresponding segment * - * Used to create iterator to segments, such as `segment_begin`. + * Given a segment index, first find its corresponding part index by performing a binary search in + * the segment offsets range. Then, skip the segment index by the number of non empty linestrings + * that precedes the current linestring to find point index to the first point of the segment. + * Dereference this point and the following point to construct the segment. * - * @tparam OffsetIterator the iterator type indicating partitions of the point range. - * @tparam CoordinateIterator the iterator type to the point range. + * @tparam OffsetIterator Iterator to the segment offsets + * @tparam CountIterator Iterator the the range of the prefix sum of non empty linestrings + * @tparam CoordinateIterator Iterator to the point range */ template struct to_valid_segment_functor { @@ -119,19 +105,76 @@ struct to_valid_segment_functor { { auto kit = thrust::prev(thrust::upper_bound(thrust::seq, segment_offset_begin, segment_offset_end, sid)); - auto geometry_id = thrust::distance(segment_offset_begin, kit); - auto preceding_non_empty_linestrings = non_empty_partitions_begin[geometry_id]; + auto part_id = thrust::distance(segment_offset_begin, kit); + auto preceding_non_empty_linestrings = non_empty_partitions_begin[part_id]; auto pid = sid + preceding_non_empty_linestrings; return segment{point_begin[pid], point_begin[pid + 1]}; } }; -/// Deduction guide for to_valid_segment_functor +/// Deduction guide template to_valid_segment_functor(OffsetIterator, OffsetIterator, CountIterator, CoordinateIterator) -> to_valid_segment_functor; +/** + * @brief Non-owning object to `multilinestring_segment` class + * + * A `multilinestring_segment_range` provide views into the segments of a multilinestring. + * The segments of a multilinestring has a near 1-1 mapping to the points of the multilinestring, + * except that the last point of a linestring and the first point of the next linestring does not + * form a valid segment. For example, the below multilinestring (points are denoted a letters): + * + * ``` + * {0, 2} + * {0, 3, 6} + * {A, B, C, X, Y, Z} + * ``` + * + * contains 6 points, but only 4 segments. AB, BC, XY and YZ. + * If we assign an index to all four segments, and an index to all points: + * + * ``` + * segments: AB BC XY YZ + * sid: 0 1 2 3 + * points: A B C X Y Z + * pid: 0 1 2 3 4 5 + * ``` + * + * Notice that if we "skip" the segment index by a few steps, it can correctly find the correponding + * point index of the starting point of the segment. For example: + * skipping sid==0 (AB) by 0 steps, finds the starting point of A (pid==0) + * skipping sid==2 (XY) by 1 step, finds the starting point of X (pid==3) + * + * Intuitively, the *steps to skip* equals the number of linestrings that precedes the linestring + * that the current segment is in. This is because every linestring adds an "invalid" segment to the + * preceding linestring. However, consider the following edge case that contains empty linestrings: + * + * ``` + * {0, 3} + * {0, 3, 3, 6} + * {A, B, C, X, Y, Z} + * ``` + * + * For segment XY, there are 2 linestrings that precedes its linestring ((0, 3) and (3, 3)). + * However, we cannot skip the sid of XY by 2 to get its starting point index. This is because the + * empty linestring in between does not introduce the "invalid" segment. Therefore, the correct + * steps to skip equals the number of *non-empty* linestrings that precedes the current linestring + * that the segment is in. + * + * Concretely, with the above example: + * ``` + * segments: AB BC XY YZ + * sid: 0 1 2 3 + * num_preceding_non_empty_linestrings: 0 0 1 1 + * skipped sid (pid): 0 0 3 4 + * starting point: A B X Y + * ``` + * + * @tparam ParentRange The multilinestring range to construct from + * @tparam IndexRange The range to the prefix sum of the non empty linestring counts + */ template class multilinestring_segment_range { using index_t = typename IndexRange::value_type; @@ -146,19 +189,26 @@ class multilinestring_segment_range { { } + /// Returns the number of segments in the multilinestring CUSPATIAL_HOST_DEVICE index_t num_segments() { return _num_segments; } + /// Returns starting iterator to the range of the starting segment index per + /// multilinestring or multipolygon CUSPATIAL_HOST_DEVICE auto multigeometry_offset_begin() { return thrust::make_permutation_iterator(_per_linestring_offset_begin(), _parent.geometry_offsets_begin()); } + /// Returns end iterator to the range of the starting segment index per multilinestring + /// or multipolygon CUSPATIAL_HOST_DEVICE auto multigeometry_offset_end() { return multigeometry_offset_begin() + _parent.num_multilinestrings() + 1; } + /// Returns starting iterator to the range of the number of segments per multilinestring of + /// multipolygon CUSPATIAL_HOST_DEVICE auto multigeometry_count_begin() { auto zipped_offset_it = thrust::make_zip_iterator(multigeometry_offset_begin(), @@ -167,11 +217,15 @@ class multilinestring_segment_range { return thrust::make_transform_iterator(zipped_offset_it, offset_pair_to_count_functor{}); } + /// Returns end iterator to the range of the number of segments per multilinestring of + /// multipolygon CUSPATIAL_HOST_DEVICE auto multigeometry_count_end() { return multigeometry_count_begin() + _parent.num_multilinestrings(); } + /// Returns the iterator to the first segment of the geometry range + /// See `to_valid_segment_functor` for implementation detail CUSPATIAL_HOST_DEVICE auto begin() { return make_counting_transform_iterator( @@ -182,6 +236,7 @@ class multilinestring_segment_range { _parent.point_begin()}); } + /// Returns the iterator to the past the last segment of the geometry range CUSPATIAL_HOST_DEVICE auto end() { return begin() + _num_segments; } private: @@ -189,14 +244,17 @@ class multilinestring_segment_range { IndexRange _non_empty_geometry_prefix_sum; index_t _num_segments; + /// Returns begin iterator to the index that points to the starting index for each linestring + /// See documentation of `to_segment_offset_iterator` for detail. CUSPATIAL_HOST_DEVICE auto _per_linestring_offset_begin() { return make_counting_transform_iterator( 0, - to_segment_offset_iterator{_parent.part_offset_begin(), - _non_empty_geometry_prefix_sum.begin()}); + point_offset_to_segment_offset{_parent.part_offset_begin(), + _non_empty_geometry_prefix_sum.begin()}); } + /// Returns end iterator to the index that points to the starting index for each linestring CUSPATIAL_HOST_DEVICE auto _per_linestring_offset_end() { return _per_linestring_offset_begin() + _non_empty_geometry_prefix_sum.size(); From f6255b817d667227c4dd0b3831b7ffee9b606887 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Fri, 19 May 2023 17:37:29 +0000 Subject: [PATCH 19/27] typo --- .../detail/range/multilinestring_segment_range.cuh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh index 44a87fec7..da13c8c9c 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh @@ -142,10 +142,10 @@ to_valid_segment_functor(OffsetIterator, OffsetIterator, CountIterator, Coordina * pid: 0 1 2 3 4 5 * ``` * - * Notice that if we "skip" the segment index by a few steps, it can correctly find the correponding - * point index of the starting point of the segment. For example: - * skipping sid==0 (AB) by 0 steps, finds the starting point of A (pid==0) - * skipping sid==2 (XY) by 1 step, finds the starting point of X (pid==3) + * Notice that if we "skip" the segment index by a few steps, it can correctly find the + * corresponding point index of the starting point of the segment. For example: skipping sid==0 (AB) + * by 0 steps, finds the starting point of A (pid==0) skipping sid==2 (XY) by 1 step, finds the + * starting point of X (pid==3) * * Intuitively, the *steps to skip* equals the number of linestrings that precedes the linestring * that the current segment is in. This is because every linestring adds an "invalid" segment to the From 0f151a587fd856bc6bcaacd3afe49573f94709e7 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Fri, 19 May 2023 18:00:35 +0000 Subject: [PATCH 20/27] Cleanups & doc improvement --- .../detail/range/multipolygon_range.cuh | 22 ------- .../cuspatial/range/multilinestring_range.cuh | 3 +- .../cuspatial/range/multipolygon_range.cuh | 7 +- cpp/tests/range/multilinestring_range_test.cu | 64 ++++++------------- 4 files changed, 25 insertions(+), 71 deletions(-) diff --git a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh index c6f9e75b7..8339831fb 100644 --- a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh @@ -355,28 +355,6 @@ multipolygon_range::g return segment{_point_begin[segment_idx], _point_begin[segment_idx + 1]}; } -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range:: - subtracted_ring_begin() -{ - return detail::make_counting_transform_iterator( - 0, detail::to_segment_offset_iterator{_ring_begin, thrust::make_counting_iterator(0)}); -} - -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range::subtracted_ring_end() -{ - return subtracted_ring_begin() + thrust::distance(_ring_begin, _ring_end); -} - template CUSPATIAL_HOST_DEVICE bool is_valid_segment_id(IndexType1 segment_idx, IndexType2 ring_idx); diff --git a/cpp/tests/range/multilinestring_range_test.cu b/cpp/tests/range/multilinestring_range_test.cu index 9e617dd7b..fe3d94dd0 100644 --- a/cpp/tests/range/multilinestring_range_test.cu +++ b/cpp/tests/range/multilinestring_range_test.cu @@ -36,23 +36,6 @@ struct MultilinestringRangeTest : public BaseFixture { std::initializer_list part_offset, std::initializer_list> coordinates, std::initializer_list> expected) - { - auto multilinestring_array = - make_multilinestring_array(geometry_offset, part_offset, coordinates); - - auto rng = multilinestring_array.range(); - - rmm::device_uvector> got(rng.num_segments(), stream()); - thrust::copy(rmm::exec_policy(stream()), rng.segment_begin(), rng.segment_end(), got.begin()); - - auto d_expected = thrust::device_vector>(expected.begin(), expected.end()); - CUSPATIAL_EXPECT_VEC2D_PAIRS_EQUIVALENT(d_expected, got); - } - - void run_segment_test_with_method_single(std::initializer_list geometry_offset, - std::initializer_list part_offset, - std::initializer_list> coordinates, - std::initializer_list> expected) { auto multilinestring_array = make_multilinestring_array(geometry_offset, part_offset, coordinates); @@ -89,11 +72,10 @@ struct MultilinestringRangeTest : public BaseFixture { CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); } - void run_multilinestring_segment_method_count_test( - std::initializer_list geometry_offset, - std::initializer_list part_offset, - std::initializer_list> coordinates, - std::initializer_list expected) + void run_multilinestring_segment_count_test(std::initializer_list geometry_offset, + std::initializer_list part_offset, + std::initializer_list> coordinates, + std::initializer_list expected) { auto multilinestring_array = make_multilinestring_array(geometry_offset, part_offset, coordinates); @@ -196,11 +178,8 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorOneSegmentTest) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_segment_test_with_method_single, - {0, 1}, - {0, 2}, - {P{0, 0}, P{1, 1}}, - {S{P{0, 0}, P{1, 1}}}); + CUSPATIAL_RUN_TEST( + this->run_segment_test_single, {0, 1}, {0, 2}, {P{0, 0}, P{1, 1}}, {S{P{0, 0}, P{1, 1}}}); } TYPED_TEST(MultilinestringRangeTest, SegmentIteratorTwoSegmentTest) @@ -209,7 +188,7 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorTwoSegmentTest) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_segment_test_with_method_single, + CUSPATIAL_RUN_TEST(this->run_segment_test_single, {0, 2}, {0, 2, 4}, {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}}, @@ -222,7 +201,7 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorTwoSegmentTest2) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_segment_test_with_method_single, + CUSPATIAL_RUN_TEST(this->run_segment_test_single, {0, 1, 2}, {0, 2, 4}, {P{0, 0}, P{1, 1}, P{0, 0}, P{1, 1}}, @@ -235,7 +214,7 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorManyPairTest) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_segment_test_with_method_single, + CUSPATIAL_RUN_TEST(this->run_segment_test_single, {0, 1, 2, 3}, {0, 6, 11, 14}, { @@ -279,7 +258,7 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorWithEmptyLineTest) using S = segment; CUSPATIAL_RUN_TEST( - this->run_segment_test_with_method_single, + this->run_segment_test_single, {0, 1, 2, 3}, {0, 3, 3, 6}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, @@ -293,7 +272,7 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorWithEmptyMultiLineStringTest using S = segment; CUSPATIAL_RUN_TEST( - this->run_segment_test_with_method_single, + this->run_segment_test_single, {0, 1, 1, 2}, {0, 3, 6}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, @@ -306,7 +285,7 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorWithEmptyMultiLineStringTest using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_segment_test_with_method_single, + CUSPATIAL_RUN_TEST(this->run_segment_test_single, {0, 1, 1, 2}, {0, 4, 7}, @@ -350,11 +329,8 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, - {0, 1}, - {0, 3}, - {P{0, 0}, P{1, 1}, P{2, 2}}, - {2}); + CUSPATIAL_RUN_TEST( + this->run_multilinestring_segment_count_test, {0, 1}, {0, 3}, {P{0, 0}, P{1, 1}, P{2, 2}}, {2}); } TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest2) @@ -363,7 +339,7 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest2) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, {0, 1, 3}, {0, 3, 3, 6}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, @@ -376,7 +352,7 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest3) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, {0, 1, 2}, {0, 3, 6}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, @@ -389,7 +365,7 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest4) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, {0, 1, 3}, {0, 3, 5, 7}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{20, 20}, P{21, 21}}, @@ -403,7 +379,7 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest5) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, {0, 1, 2, 3}, {0, 3, 3, 6}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, @@ -417,7 +393,7 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest6) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, {0, 1, 1, 2}, {0, 3, 6}, {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, @@ -431,7 +407,7 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest7) using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_method_count_test, + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, {0, 1, 1, 2}, {0, 4, 7}, {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}, P{10, 10}, P{11, 11}, P{12, 12}}, From a45f2c79c38a449e8268117c082db0943ed9be5f Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Fri, 19 May 2023 18:16:27 +0000 Subject: [PATCH 21/27] [skip ci] initial addition of the load balanced kernel --- .../detail/algorithm/linestring_distance.cuh | 33 ++++++++++++++++--- 1 file changed, 29 insertions(+), 4 deletions(-) diff --git a/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh b/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh index 4c988612d..0c6086e09 100644 --- a/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh +++ b/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh @@ -42,10 +42,10 @@ namespace detail { * set to nullopt, no distance computation will be bypassed. */ template -__global__ void linestring_distance(MultiLinestringRange1 multilinestrings1, - MultiLinestringRange2 multilinestrings2, - thrust::optional intersects, - OutputIt distances_first) +__global__ void linestring_distance_multilinestring_loop(MultiLinestringRange1 multilinestrings1, + MultiLinestringRange2 multilinestrings2, + thrust::optional intersects, + OutputIt distances_first) { using T = typename MultiLinestringRange1::element_t; @@ -72,5 +72,30 @@ __global__ void linestring_distance(MultiLinestringRange1 multilinestrings1, } } +/** + * @internal + * @brief The kernel to compute (multi)linestring to (multi)linestring distance + * + * Load balanced kernel to compute distances between one pair of segments from the multilinestring + * and multilinestring. + * + * `intersects` is an optional pointer to a boolean range where the `i`th element indicates the + * `i`th output should be set to 0 and bypass distance computation. This argument is optional, if + * set to nullopt, no distance computation will be bypassed. + */ +template +__global__ void linestring_distance_load_balanced(MultiLinestringRange1 multilinestrings1, + MultiLinestringRange2 multilinestrings2, + thrust::optional intersects, + OutputIt distances_first) +{ + using T = typename MultiLinestringRange1::element_t; + + for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < multilinestrings1.num_points(); + idx += gridDim.x * blockDim.x) { + atomicMin(&distances_first[geometry_idx], static_cast(sqrt(min_distance_squared))); + } +} + } // namespace detail } // namespace cuspatial From 8475fe83819b90d2cfe1fd410c3433ddde03aeb8 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Mon, 22 May 2023 21:47:33 +0000 Subject: [PATCH 22/27] Move detail iterator factory into detail folder --- .../cuspatial/detail/iterator_factory.cuh | 203 ++++++++++++++++++ cpp/include/cuspatial/iterator_factory.cuh | 166 +------------- 2 files changed, 211 insertions(+), 158 deletions(-) create mode 100644 cpp/include/cuspatial/detail/iterator_factory.cuh diff --git a/cpp/include/cuspatial/detail/iterator_factory.cuh b/cpp/include/cuspatial/detail/iterator_factory.cuh new file mode 100644 index 000000000..8f66de0ed --- /dev/null +++ b/cpp/include/cuspatial/detail/iterator_factory.cuh @@ -0,0 +1,203 @@ +/* + * 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 +#include +#include +#include +#include +#include +#include + +#include + +namespace cuspatial { +namespace detail { + +/** + * @internal + * @brief Helper to create a `transform_iterator` that transforms sequential values. + */ +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); +} + +/** + * @internal + * @brief Helper to convert a tuple of elements into a `vec_2d` + */ +template > +struct tuple_to_vec_2d { + __device__ VectorType operator()(thrust::tuple const& pos) + { + return VectorType{thrust::get<0>(pos), thrust::get<1>(pos)}; + } +}; + +/** + * @internal + * @brief Helper to convert a `vec_2d` into a tuple of elements + */ +template > +struct vec_2d_to_tuple { + __device__ thrust::tuple operator()(VectorType const& xy) + { + return thrust::make_tuple(xy.x, xy.y); + } +}; + +/** + * @internal + * @brief Helper to convert a `box` into a tuple of elements + */ +template > +struct box_to_tuple { + __device__ thrust::tuple operator()(Box const& box) + { + return thrust::make_tuple(box.v1.x, box.v1.y, box.v2.x, box.v2.y); + } +}; + +/** + * @internal + * @brief Helper to convert a tuple of vec_2d into a `box` + */ +template > +struct vec_2d_tuple_to_box { + __device__ box operator()(thrust::tuple pair) + { + auto v1 = thrust::get<0>(pair); + auto v2 = thrust::get<1>(pair); + return {v1, v2}; + } +}; + +/** + * @internal + * @brief Generic to convert any iterator pointing to interleaved xy range into + * iterator of vec_2d. + * + * This generic version does not use of vectorized load. + * + * @pre `Iter` has operator[] defined. + * @pre std::iterator_traits::value_type is convertible to `T`. + */ +template +struct interleaved_to_vec_2d { + using element_t = typename std::iterator_traits::value_type; + using value_type = vec_2d; + Iter it; + constexpr interleaved_to_vec_2d(Iter it) : it{it} {} + + CUSPATIAL_HOST_DEVICE value_type operator()(std::size_t i) + { + return vec_2d{it[2 * i], it[2 * i + 1]}; + } +}; + +/** + * @brief Specialization for thrust iterators conforming to `contiguous_iterator`. (including raw + * pointer) + * + * This iterator specific version uses vectorized load. + * + * @throw cuspatial::logic_error if `Iter` is not aligned to type `vec_2d` + * @pre `Iter` is a `contiguous_iterator` (including raw pointer). + */ +template +struct interleaved_to_vec_2d>> { + using element_t = typename std::iterator_traits::value_type; + using value_type = vec_2d; + + element_t const* ptr; + + constexpr interleaved_to_vec_2d(Iter it) : ptr{&thrust::raw_reference_cast(*it)} + { + CUSPATIAL_EXPECTS(!((intptr_t)ptr % alignof(vec_2d)), + "Misaligned interleaved data."); + } + + CUSPATIAL_HOST_DEVICE value_type operator()(std::size_t i) + { + auto const aligned = + static_cast(__builtin_assume_aligned(ptr + 2 * i, 2 * sizeof(element_t))); + return vec_2d{aligned[0], aligned[1]}; + } +}; + +/** + * @internal + * @brief Functor to transform an index to strided index. + */ +template +struct strided_functor { + auto __device__ operator()(std::size_t i) { return i * stride; } +}; + +/** + * @internal + * @brief Functor to transform an index into a geometry ID determined by a range of offsets. + */ +template +struct index_to_geometry_id { + GeometryIter geometry_begin; + GeometryIter geometry_end; + + index_to_geometry_id(GeometryIter begin, GeometryIter end) + : geometry_begin(begin), geometry_end(end) + { + } + + CUSPATIAL_HOST_DEVICE auto operator()(IndexT idx) + { + return thrust::distance(geometry_begin, + thrust::upper_bound(thrust::seq, geometry_begin, geometry_end, idx)); + } +}; + +/** + * @brief Given iterator a pair of offsets, return the number of elements between the offsets. + * + * Used to create iterator to geometry counts, such as `multi*_point_count_begin`, + * `multi*_segment_count_begin`. + * + * Example: + * pair of offsets: (0, 3), (3, 5), (5, 8) + * number of elements between offsets: 3, 2, 3 + * + * @tparam OffsetPairIterator Must be iterator type to thrust::pair of indices. + * @param p Iterator of thrust::pair of indices. + */ +struct offset_pair_to_count_functor { + template + CUSPATIAL_HOST_DEVICE auto operator()(OffsetPairIterator p) + { + return thrust::get<1>(p) - thrust::get<0>(p); + } +}; + +} // namespace detail +} // namespace cuspatial diff --git a/cpp/include/cuspatial/iterator_factory.cuh b/cpp/include/cuspatial/iterator_factory.cuh index 1f026512c..49313298b 100644 --- a/cpp/include/cuspatial/iterator_factory.cuh +++ b/cpp/include/cuspatial/iterator_factory.cuh @@ -16,171 +16,14 @@ #pragma once -#include -#include -#include -#include +#include -#include -#include -#include #include #include #include #include -#include -#include - -#include namespace cuspatial { -namespace detail { - -/** - * @internal - * @brief Helper to create a `transform_iterator` that transforms sequential values. - */ -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); -} - -/** - * @internal - * @brief Helper to convert a tuple of elements into a `vec_2d` - */ -template > -struct tuple_to_vec_2d { - __device__ VectorType operator()(thrust::tuple const& pos) - { - return VectorType{thrust::get<0>(pos), thrust::get<1>(pos)}; - } -}; - -/** - * @internal - * @brief Helper to convert a `vec_2d` into a tuple of elements - */ -template > -struct vec_2d_to_tuple { - __device__ thrust::tuple operator()(VectorType const& xy) - { - return thrust::make_tuple(xy.x, xy.y); - } -}; - -/** - * @internal - * @brief Helper to convert a `box` into a tuple of elements - */ -template > -struct box_to_tuple { - __device__ thrust::tuple operator()(Box const& box) - { - return thrust::make_tuple(box.v1.x, box.v1.y, box.v2.x, box.v2.y); - } -}; - -/** - * @internal - * @brief Helper to convert a tuple of vec_2d into a `box` - */ -template > -struct vec_2d_tuple_to_box { - __device__ box operator()(thrust::tuple pair) - { - auto v1 = thrust::get<0>(pair); - auto v2 = thrust::get<1>(pair); - return {v1, v2}; - } -}; - -/** - * @internal - * @brief Generic to convert any iterator pointing to interleaved xy range into - * iterator of vec_2d. - * - * This generic version does not use of vectorized load. - * - * @pre `Iter` has operator[] defined. - * @pre std::iterator_traits::value_type is convertible to `T`. - */ -template -struct interleaved_to_vec_2d { - using element_t = typename std::iterator_traits::value_type; - using value_type = vec_2d; - Iter it; - constexpr interleaved_to_vec_2d(Iter it) : it{it} {} - - CUSPATIAL_HOST_DEVICE value_type operator()(std::size_t i) - { - return vec_2d{it[2 * i], it[2 * i + 1]}; - } -}; - -/** - * @brief Specialization for thrust iterators conforming to `contiguous_iterator`. (including raw - * pointer) - * - * This iterator specific version uses vectorized load. - * - * @throw cuspatial::logic_error if `Iter` is not aligned to type `vec_2d` - * @pre `Iter` is a `contiguous_iterator` (including raw pointer). - */ -template -struct interleaved_to_vec_2d>> { - using element_t = typename std::iterator_traits::value_type; - using value_type = vec_2d; - - element_t const* ptr; - - constexpr interleaved_to_vec_2d(Iter it) : ptr{&thrust::raw_reference_cast(*it)} - { - CUSPATIAL_EXPECTS(!((intptr_t)ptr % alignof(vec_2d)), - "Misaligned interleaved data."); - } - - CUSPATIAL_HOST_DEVICE value_type operator()(std::size_t i) - { - auto const aligned = - static_cast(__builtin_assume_aligned(ptr + 2 * i, 2 * sizeof(element_t))); - return vec_2d{aligned[0], aligned[1]}; - } -}; - -/** - * @internal - * @brief Functor to transform an index to strided index. - */ -template -struct strided_functor { - auto __device__ operator()(std::size_t i) { return i * stride; } -}; - -/** - * @internal - * @brief Functor to transform an index into a geometry ID determined by a range of offsets. - */ -template -struct index_to_geometry_id { - GeometryIter geometry_begin; - GeometryIter geometry_end; - - index_to_geometry_id(GeometryIter begin, GeometryIter end) - : geometry_begin(begin), geometry_end(end) - { - } - - CUSPATIAL_HOST_DEVICE auto operator()(IndexT idx) - { - return thrust::distance(geometry_begin, - thrust::upper_bound(thrust::seq, geometry_begin, geometry_end, idx)); - } -}; - -} // namespace detail /** * @addtogroup type_factories @@ -424,6 +267,13 @@ auto make_geometry_id_iterator(GeometryIter geometry_offsets_begin, std::distance(geometry_offsets_begin, geometry_offsets_end))); } +template +CUSPATIAL_HOST_DEVICE auto make_element_count_iterator(OffsetIterator iter) +{ + auto zipped_it = thrust::make_zip_iterator(iter, thrust::next(iter)); + return thrust::make_transform_iterator(zipped_it, detail::offset_pair_to_count_functor{}); +} + /** * @} // end of doxygen group */ From c463978c7ec93aac65f1bc72d0c17ad9a79401f5 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Mon, 22 May 2023 21:47:47 +0000 Subject: [PATCH 23/27] Remove functors.cuh, make iterator factory for it. --- cpp/include/cuspatial/detail/functors.cuh | 50 ----------------------- 1 file changed, 50 deletions(-) delete mode 100644 cpp/include/cuspatial/detail/functors.cuh diff --git a/cpp/include/cuspatial/detail/functors.cuh b/cpp/include/cuspatial/detail/functors.cuh deleted file mode 100644 index 25c79885d..000000000 --- a/cpp/include/cuspatial/detail/functors.cuh +++ /dev/null @@ -1,50 +0,0 @@ -/* - * Copyright (c) 2023, NVIDIA CORPORATION. - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#pragma once - -#include -#include -#include - -#include - -namespace cuspatial { -namespace detail { - -/** - * @brief Given iterator a pair of offsets, return the number of elements between the offsets. - * - * Used to create iterator to geometry counts, such as `multi*_point_count_begin`, - * `multi*_segment_count_begin`. - * - * Example: - * pair of offsets: (0, 3), (3, 5), (5, 8) - * number of elements between offsets: 3, 2, 3 - * - * @tparam OffsetPairIterator Must be iterator type to thrust::pair of indices. - * @param p Iterator of thrust::pair of indices. - */ -struct offset_pair_to_count_functor { - template - CUSPATIAL_HOST_DEVICE auto operator()(OffsetPairIterator p) - { - return thrust::get<1>(p) - thrust::get<0>(p); - } -}; - -} // namespace detail -} // namespace cuspatial From b442f7c1538e29745a661fe50e84f0fbf7fbec73 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Mon, 22 May 2023 21:48:16 +0000 Subject: [PATCH 24/27] Remove every functors usage --- .../detail/distance/linestring_polygon_distance.cuh | 1 - .../cuspatial/detail/multilinestring_segment.cuh | 1 - .../detail/range/multilinestring_range.cuh | 8 ++------ .../detail/range/multilinestring_segment_range.cuh | 7 ++----- .../cuspatial/detail/range/multipolygon_range.cuh | 13 ++----------- 5 files changed, 6 insertions(+), 24 deletions(-) diff --git a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh index 1c5cb26d8..345f6d76e 100644 --- a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh +++ b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh @@ -19,7 +19,6 @@ #include "distance_utils.cuh" #include -#include #include #include #include diff --git a/cpp/include/cuspatial/detail/multilinestring_segment.cuh b/cpp/include/cuspatial/detail/multilinestring_segment.cuh index 0f6215117..d57f125dc 100644 --- a/cpp/include/cuspatial/detail/multilinestring_segment.cuh +++ b/cpp/include/cuspatial/detail/multilinestring_segment.cuh @@ -16,7 +16,6 @@ #pragma once -#include #include #include #include diff --git a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh index 244e622e7..3b0aa597f 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh @@ -17,7 +17,6 @@ #pragma once #include -#include #include #include #include @@ -213,9 +212,7 @@ CUSPATIAL_HOST_DEVICE auto multilinestring_range @@ -229,8 +226,7 @@ template :: multilinestring_linestring_count_begin() { - auto paired_it = thrust::make_zip_iterator(_geometry_begin, thrust::next(_geometry_begin)); - return thrust::make_transform_iterator(paired_it, detail::offset_pair_to_count_functor{}); + return make_element_count_iterator(_geometry_begin); } template diff --git a/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh index da13c8c9c..58817c46e 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh @@ -17,7 +17,7 @@ #pragma once #include -#include +#include #include #include @@ -211,10 +211,7 @@ class multilinestring_segment_range { /// multipolygon CUSPATIAL_HOST_DEVICE auto multigeometry_count_begin() { - auto zipped_offset_it = thrust::make_zip_iterator(multigeometry_offset_begin(), - thrust::next(multigeometry_offset_begin())); - - return thrust::make_transform_iterator(zipped_offset_it, offset_pair_to_count_functor{}); + return make_element_count_iterator(multigeometry_offset_begin()); } /// Returns end iterator to the range of the number of segments per multilinestring of diff --git a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh index 8339831fb..567105fd1 100644 --- a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh @@ -17,7 +17,6 @@ #pragma once #include -#include #include #include #include @@ -269,11 +268,7 @@ multipolygon_range:: auto multipolygon_point_offset_it = thrust::make_permutation_iterator( _ring_begin, thrust::make_permutation_iterator(_part_begin, _geometry_begin)); - auto point_offset_pair_it = thrust::make_zip_iterator(multipolygon_point_offset_it, - thrust::next(multipolygon_point_offset_it)); - - return thrust::make_transform_iterator(point_offset_pair_it, - detail::offset_pair_to_count_functor{}); + return make_element_count_iterator(multipolygon_point_offset_it); } template :: auto multipolygon_ring_offset_it = thrust::make_permutation_iterator(_part_begin, _geometry_begin); - auto ring_offset_pair_it = thrust::make_zip_iterator(multipolygon_ring_offset_it, - thrust::next(multipolygon_ring_offset_it)); - - return thrust::make_transform_iterator(ring_offset_pair_it, - detail::offset_pair_to_count_functor{}); + return make_element_count_iterator(multipolygon_ring_offset_it); } template Date: Mon, 22 May 2023 21:49:06 +0000 Subject: [PATCH 25/27] Impelement load balanced linestring distance --- .../detail/algorithm/linestring_distance.cuh | 49 ++++++++++++++----- .../detail/distance/distance_utils.cuh | 41 ++++++++++++++++ .../detail/distance/linestring_distance.cuh | 41 +++++++++++----- 3 files changed, 109 insertions(+), 22 deletions(-) diff --git a/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh b/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh index 0c6086e09..7aec99d8e 100644 --- a/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh +++ b/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh @@ -21,6 +21,7 @@ #include +#include #include namespace cuspatial { @@ -42,10 +43,10 @@ namespace detail { * set to nullopt, no distance computation will be bypassed. */ template -__global__ void linestring_distance_multilinestring_loop(MultiLinestringRange1 multilinestrings1, - MultiLinestringRange2 multilinestrings2, - thrust::optional intersects, - OutputIt distances_first) +__global__ void linestring_distance(MultiLinestringRange1 multilinestrings1, + MultiLinestringRange2 multilinestrings2, + thrust::optional intersects, + OutputIt distances_first) { using T = typename MultiLinestringRange1::element_t; @@ -83,17 +84,43 @@ __global__ void linestring_distance_multilinestring_loop(MultiLinestringRange1 m * `i`th output should be set to 0 and bypass distance computation. This argument is optional, if * set to nullopt, no distance computation will be bypassed. */ -template -__global__ void linestring_distance_load_balanced(MultiLinestringRange1 multilinestrings1, - MultiLinestringRange2 multilinestrings2, +template +__global__ void linestring_distance_load_balanced(SegmentRange1 multilinestrings1, + SegmentRange2 multilinestrings2, + IndexRange thread_bounds, thrust::optional intersects, - OutputIt distances_first) + OutputIt distances) { - using T = typename MultiLinestringRange1::element_t; + using index_t = typename IndexRange::value_type; - for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < multilinestrings1.num_points(); + auto num_threads = thread_bounds[thread_bounds.size() - 1]; + for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < num_threads; idx += gridDim.x * blockDim.x) { - atomicMin(&distances_first[geometry_idx], static_cast(sqrt(min_distance_squared))); + auto it = thrust::prev( + thrust::upper_bound(thrust::seq, thread_bounds.begin(), thread_bounds.end(), idx)); + auto const geometry_id = thrust::distance(thread_bounds.begin(), it); + + if (intersects.has_value() && intersects.value()[geometry_id]) { + distances[geometry_id] = 0.0f; + continue; + } + + auto const local_idx = idx - *it; + // Retrieve the number of segments in multilinestrings[geometry_id] + auto num_segment_this_multilinestring = + multilinestrings1.multigeometry_count_begin()[geometry_id]; + // The segment id from multilinestring1 this thread is computing (local_id + global_offset) + auto multilinestring1_segment_id = local_idx % num_segment_this_multilinestring + + multilinestrings1.multigeometry_offset_begin()[geometry_id]; + + // The segment id from multilinestring2 this thread is computing (local_id + global_offset) + auto multilinestring2_segment_id = local_idx / num_segment_this_multilinestring + + multilinestrings2.multigeometry_offset_begin()[geometry_id]; + + auto [a, b] = multilinestrings1.begin()[multilinestring1_segment_id]; + auto [c, d] = multilinestrings2.begin()[multilinestring2_segment_id]; + + atomicMin(&distances[geometry_id], sqrt(squared_segment_distance(a, b, c, d))); } } diff --git a/cpp/include/cuspatial/detail/distance/distance_utils.cuh b/cpp/include/cuspatial/detail/distance/distance_utils.cuh index 71ee9d5fd..13860f600 100644 --- a/cpp/include/cuspatial/detail/distance/distance_utils.cuh +++ b/cpp/include/cuspatial/detail/distance/distance_utils.cuh @@ -29,6 +29,7 @@ #include #include #include +#include namespace cuspatial { namespace detail { @@ -115,5 +116,45 @@ rmm::device_uvector point_polygon_intersects(MultiPointRange multipoint return multipoint_intersects; } +/** + * @brief Compute the thread bound between two ranges of partitions + * + * @tparam CountIterator1 + * @tparam CountIterator2 + * @param lhs + * @param rhs + * @param stream + * @return rmm::device_uvector + */ +template > +rmm::device_uvector compute_segment_thread_bounds(CountIterator1 lhs_begin, + CountIterator1 lhs_end, + CountIterator2 rhs_begin, + rmm::cuda_stream_view stream) +{ + auto size = thrust::distance(lhs_begin, lhs_end) + 1; + + // Compute the "boundary" of threads. Threads are partitioned based on the number of linestrings + // times the number of polygons in a multilinestring-multipolygon pair. + auto segment_count_product_it = + thrust::make_transform_iterator(thrust::make_zip_iterator(lhs_begin, rhs_begin), + thrust::make_zip_function(thrust::multiplies{})); + + // Computes the "thread boundary" of each pair. This array partitions the thread range by + // geometries. E.g. threadIdx within [thread_bounds[i], thread_bounds[i+1]) computes distances of + // the ith pair. + auto thread_bounds = rmm::device_uvector(size, stream); + detail::zero_data_async(thread_bounds.begin(), thread_bounds.end(), stream); + + thrust::inclusive_scan(rmm::exec_policy(stream), + segment_count_product_it, + segment_count_product_it + thread_bounds.size() - 1, + thrust::next(thread_bounds.begin())); + + return thread_bounds; +} + } // namespace detail } // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/distance/linestring_distance.cuh b/cpp/include/cuspatial/detail/distance/linestring_distance.cuh index 19e70b51b..20159924c 100644 --- a/cpp/include/cuspatial/detail/distance/linestring_distance.cuh +++ b/cpp/include/cuspatial/detail/distance/linestring_distance.cuh @@ -17,6 +17,7 @@ #pragma once #include +#include #include #include #include @@ -33,8 +34,8 @@ namespace cuspatial { template -OutputIt pairwise_linestring_distance(MultiLinestringRange1 multilinestrings1, - MultiLinestringRange2 multilinestrings2, +OutputIt pairwise_linestring_distance(MultiLinestringRange1 lhs, + MultiLinestringRange2 rhs, OutputIt distances_first, rmm::cuda_stream_view stream) { @@ -48,25 +49,43 @@ OutputIt pairwise_linestring_distance(MultiLinestringRange1 multilinestrings1, typename MultiLinestringRange2::point_t>(), "All input types must be cuspatial::vec_2d with the same value type"); - CUSPATIAL_EXPECTS(multilinestrings1.size() == multilinestrings2.size(), - "Inputs must have the same number of rows."); + CUSPATIAL_EXPECTS(lhs.size() == rhs.size(), "Inputs must have the same number of rows."); - if (multilinestrings1.size() == 0) return distances_first; + if (lhs.size() == 0) return distances_first; + // Make views to the segments in the multilinestring + auto lhs_segments = lhs._segments(stream); + auto lhs_segments_range = lhs_segments.view(); + + // Make views to the segments in the multilinestring + auto rhs_segments = rhs._segments(stream); + auto rhs_segments_range = rhs_segments.view(); + + auto thread_bounds = + detail::compute_segment_thread_bounds(lhs_segments_range.multigeometry_count_begin(), + lhs_segments_range.multigeometry_count_end(), + rhs_segments_range.multigeometry_count_begin(), + stream); + // Initialize the output range thrust::fill(rmm::exec_policy(stream), distances_first, - distances_first + multilinestrings1.size(), + distances_first + lhs.size(), std::numeric_limits::max()); std::size_t constexpr threads_per_block = 256; - std::size_t const num_blocks = - (multilinestrings1.num_points() + threads_per_block - 1) / threads_per_block; + std::size_t num_threads = thread_bounds.element(thread_bounds.size() - 1, stream); + std::size_t const num_blocks = (num_threads + threads_per_block - 1) / threads_per_block; - detail::linestring_distance<<>>( - multilinestrings1, multilinestrings2, thrust::nullopt, distances_first); + detail::linestring_distance_load_balanced + <<>>( + lhs_segments_range, + rhs_segments_range, + range(thread_bounds.begin(), thread_bounds.end()), + thrust::nullopt, + distances_first); CUSPATIAL_CUDA_TRY(cudaGetLastError()); - return distances_first + multilinestrings1.size(); + return distances_first + lhs.size(); } } // namespace cuspatial From 3621efd76a4adb4a6cc9361dfc37d51e9fc7ca8c Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 30 May 2023 11:06:12 -0700 Subject: [PATCH 26/27] optimize load_balanced_kernel --- .../detail/algorithm/linestring_distance.cuh | 44 ++++++++++++++++--- .../detail/distance/linestring_distance.cuh | 6 ++- 2 files changed, 43 insertions(+), 7 deletions(-) diff --git a/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh b/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh index 7aec99d8e..836cc18db 100644 --- a/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh +++ b/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh @@ -24,6 +24,10 @@ #include #include +#include + +#include + namespace cuspatial { namespace detail { @@ -73,6 +77,15 @@ __global__ void linestring_distance(MultiLinestringRange1 multilinestrings1, } } +template +auto __device__ compute_geometry_id(BoundsIterator bounds_begin, BoundsIterator bounds_end, IndexType idx) +{ + auto it = thrust::prev( + thrust::upper_bound(thrust::seq, bounds_begin, bounds_end, idx)); + auto const geometry_id = thrust::distance(bounds_begin, it); + return thrust::make_pair(it, geometry_id); +} + /** * @internal * @brief The kernel to compute (multi)linestring to (multi)linestring distance @@ -84,7 +97,7 @@ __global__ void linestring_distance(MultiLinestringRange1 multilinestrings1, * `i`th output should be set to 0 and bypass distance computation. This argument is optional, if * set to nullopt, no distance computation will be bypassed. */ -template +template __global__ void linestring_distance_load_balanced(SegmentRange1 multilinestrings1, SegmentRange2 multilinestrings2, IndexRange thread_bounds, @@ -96,9 +109,14 @@ __global__ void linestring_distance_load_balanced(SegmentRange1 multilinestrings auto num_threads = thread_bounds[thread_bounds.size() - 1]; for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < num_threads; idx += gridDim.x * blockDim.x) { - auto it = thrust::prev( - thrust::upper_bound(thrust::seq, thread_bounds.begin(), thread_bounds.end(), idx)); - auto const geometry_id = thrust::distance(thread_bounds.begin(), it); + + auto [it, geometry_id] = compute_geometry_id(thread_bounds.begin(), thread_bounds.end(), idx); + + auto first_thread_of_this_block = blockDim.x * blockIdx.x; + auto last_thread_of_block = blockDim.x * (blockIdx.x + 1) - 1; + auto first_thread_of_next_geometry = thread_bounds[geometry_id + 1]; + + bool split_block = first_thread_of_this_block < first_thread_of_next_geometry && first_thread_of_next_geometry <= last_thread_of_block; if (intersects.has_value() && intersects.value()[geometry_id]) { distances[geometry_id] = 0.0f; @@ -120,7 +138,23 @@ __global__ void linestring_distance_load_balanced(SegmentRange1 multilinestrings auto [a, b] = multilinestrings1.begin()[multilinestring1_segment_id]; auto [c, d] = multilinestrings2.begin()[multilinestring2_segment_id]; - atomicMin(&distances[geometry_id], sqrt(squared_segment_distance(a, b, c, d))); + auto partial = sqrt(squared_segment_distance(a, b, c, d)); + + if (split_block) + { + atomicMin(&distances[geometry_id], partial); + } + else + { + // block reduce + typedef cub::BlockReduce BlockReduce; + __shared__ typename BlockReduce::TempStorage temp_storage; + auto aggregate = BlockReduce(temp_storage).Reduce(partial, cub::Min()); + + // atmomic with leading thread + if (cooperative_groups::this_thread_block().thread_rank() == 0) + atomicMin(&distances[geometry_id], aggregate); + } } } diff --git a/cpp/include/cuspatial/detail/distance/linestring_distance.cuh b/cpp/include/cuspatial/detail/distance/linestring_distance.cuh index 20159924c..a31ea9dea 100644 --- a/cpp/include/cuspatial/detail/distance/linestring_distance.cuh +++ b/cpp/include/cuspatial/detail/distance/linestring_distance.cuh @@ -21,6 +21,7 @@ #include #include #include +#include #include #include @@ -73,10 +74,11 @@ OutputIt pairwise_linestring_distance(MultiLinestringRange1 lhs, std::numeric_limits::max()); std::size_t constexpr threads_per_block = 256; - std::size_t num_threads = thread_bounds.element(thread_bounds.size() - 1, stream); + // std::size_t num_threads = thread_bounds.element(thread_bounds.size() - 1, stream); + std::size_t num_threads = lhs.num_points() * rhs.num_points(); std::size_t const num_blocks = (num_threads + threads_per_block - 1) / threads_per_block; - detail::linestring_distance_load_balanced + detail::linestring_distance_load_balanced <<>>( lhs_segments_range, rhs_segments_range, From 4633224611b67accab846f19e580ced2ccb1e053 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 30 May 2023 12:13:57 -0700 Subject: [PATCH 27/27] minor updates to balanced kernel --- .../cuspatial/detail/algorithm/linestring_distance.cuh | 6 ++---- .../cuspatial/detail/distance/linestring_distance.cuh | 3 +-- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh b/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh index 836cc18db..d8a4c0357 100644 --- a/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh +++ b/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh @@ -106,8 +106,8 @@ __global__ void linestring_distance_load_balanced(SegmentRange1 multilinestrings { using index_t = typename IndexRange::value_type; - auto num_threads = thread_bounds[thread_bounds.size() - 1]; - for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < num_threads; + auto num_segment_pairs = thread_bounds[thread_bounds.size() - 1]; + for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < num_segment_pairs; idx += gridDim.x * blockDim.x) { auto [it, geometry_id] = compute_geometry_id(thread_bounds.begin(), thread_bounds.end(), idx); @@ -141,9 +141,7 @@ __global__ void linestring_distance_load_balanced(SegmentRange1 multilinestrings auto partial = sqrt(squared_segment_distance(a, b, c, d)); if (split_block) - { atomicMin(&distances[geometry_id], partial); - } else { // block reduce diff --git a/cpp/include/cuspatial/detail/distance/linestring_distance.cuh b/cpp/include/cuspatial/detail/distance/linestring_distance.cuh index a31ea9dea..dc901e78b 100644 --- a/cpp/include/cuspatial/detail/distance/linestring_distance.cuh +++ b/cpp/include/cuspatial/detail/distance/linestring_distance.cuh @@ -74,8 +74,7 @@ OutputIt pairwise_linestring_distance(MultiLinestringRange1 lhs, std::numeric_limits::max()); std::size_t constexpr threads_per_block = 256; - // std::size_t num_threads = thread_bounds.element(thread_bounds.size() - 1, stream); - std::size_t num_threads = lhs.num_points() * rhs.num_points(); + std::size_t num_threads = 1e8; std::size_t const num_blocks = (num_threads + threads_per_block - 1) / threads_per_block; detail::linestring_distance_load_balanced