Merge pull request #1203 from aprokop/dbscan_float

benchmarks/dbscan/ArborX_DBSCANVerification.hpp

@@ -295,9 +295,10 @@ struct IndexOnlyCallback
   }
 };
 
-template <typename ExecutionSpace, typename Primitives, typename LabelsView>
+template <typename ExecutionSpace, typename Primitives, typename LabelsView,
+          typename Coordinate>
 bool verifyDBSCAN(ExecutionSpace exec_space, Primitives const &primitives,
-                  float eps, int core_min_size, LabelsView const &labels,
+                  Coordinate eps, int core_min_size, LabelsView const &labels,
                   bool verbose = false)
 {
   Kokkos::Profiling::pushRegion("ArborX::DBSCAN::verify");

src/cluster/ArborX_DBSCAN.hpp

@@ -52,9 +52,10 @@ struct DBSCANCorePoints
   }
 };
 
+template <typename Coordinate>
 struct WithinRadiusGetter
 {
-  float _r;
+  Coordinate _r;
 
   template <typename Point, typename Index>
   KOKKOS_FUNCTION auto
@@ -63,18 +64,17 @@ struct WithinRadiusGetter
     static_assert(GeometryTraits::is_point_v<Point>);
 
     constexpr int DIM = GeometryTraits::dimension_v<Point>;
-    using Coordinate = GeometryTraits::coordinate_type_t<Point>;
     using ArborX::intersects;
     return intersects(
         Sphere{convert<::ArborX::Point<DIM, Coordinate>>(pair.value), _r});
   }
 };
 
-template <typename Primitives, typename PermuteFilter>
-struct PrimitivesWithRadiusReorderedAndFiltered
+template <typename Points, typename PermuteFilter>
+struct PointsWithRadiusReorderedAndFiltered
 {
-  Primitives _primitives;
-  float _r;
+  Points _points;
+  GeometryTraits::coordinate_type_t<typename Points::value_type> _r;
   PermuteFilter _filter;
 };
 
@@ -84,8 +84,11 @@ template <typename Points, typename DenseCellOffsets, typename CellIndices,
           typename Permutation>
 struct MixedBoxPrimitives
 {
+  using Point = typename Points::value_type;
   Points _points;
-  CartesianGrid<GeometryTraits::dimension_v<typename Points::value_type>> _grid;
+  CartesianGrid<GeometryTraits::dimension_v<Point>,
+                GeometryTraits::coordinate_type_t<Point>>
+      _grid;
   DenseCellOffsets _dense_cell_offsets;
   int _num_points_in_dense_cells; // to avoid lastElement() in AccessTraits
   CellIndices _sorted_cell_indices;
@@ -95,13 +98,12 @@ struct MixedBoxPrimitives
 } // namespace Details
 
 template <typename Primitives, typename PermuteFilter>
-struct AccessTraits<Details::PrimitivesWithRadiusReorderedAndFiltered<
-    Primitives, PermuteFilter>>
+struct AccessTraits<
+    Details::PointsWithRadiusReorderedAndFiltered<Primitives, PermuteFilter>>
 {
   using memory_space = typename Primitives::memory_space;
   using Predicates =
-      Details::PrimitivesWithRadiusReorderedAndFiltered<Primitives,
-                                                        PermuteFilter>;
+      Details::PointsWithRadiusReorderedAndFiltered<Primitives, PermuteFilter>;
 
   static KOKKOS_FUNCTION size_t size(Predicates const &w)
   {
@@ -110,13 +112,14 @@ struct AccessTraits<Details::PrimitivesWithRadiusReorderedAndFiltered<
   static KOKKOS_FUNCTION auto get(Predicates const &w, size_t i)
   {
     int index = w._filter(i);
-    auto const &point = w._primitives(index);
-    constexpr int dim =
-        GeometryTraits::dimension_v<std::decay_t<decltype(point)>>;
+    auto const &point = w._points(index);
+    using Point = std::decay_t<decltype(point)>;
+    constexpr int DIM = GeometryTraits::dimension_v<Point>;
+    using Coordinate = GeometryTraits::coordinate_type_t<Point>;
     // FIXME reinterpret_cast is dangerous here if access traits return user
     // point structure (e.g., struct MyPoint { float y; float x; })
     auto const &hyper_point =
-        reinterpret_cast<::ArborX::Point<dim> const &>(point);
+        reinterpret_cast<::ArborX::Point<DIM, Coordinate> const &>(point);
     return attach(intersects(Sphere{hyper_point, w._r}), (int)index);
   }
 };
@@ -158,12 +161,16 @@ struct AccessTraits<
     i = (i - num_dense_primitives) + w._num_points_in_dense_cells;
 
     auto const &point = w._points(w._permute(i));
-    constexpr int dim =
-        GeometryTraits::dimension_v<std::decay_t<decltype(point)>>;
+
+    using Point = std::decay_t<decltype(point)>;
+    constexpr int DIM = GeometryTraits::dimension_v<Point>;
+    using Coordinate = typename GeometryTraits::coordinate_type<Point>::type;
+
     // FIXME reinterpret_cast is dangerous here if access traits return user
     // point structure (e.g., struct MyPoint { float y; float x; })
-    auto const &hyper_point = reinterpret_cast<Point<dim> const &>(point);
-    return Box<dim>{hyper_point, hyper_point};
+    auto const &hyper_point =
+        reinterpret_cast<::ArborX::Point<DIM, Coordinate> const &>(point);
+    return Box{hyper_point, hyper_point};
   }
   using memory_space = typename MixedOffsets::memory_space;
 };
@@ -197,10 +204,10 @@ struct Parameters
 };
 } // namespace DBSCAN
 
-template <typename ExecutionSpace, typename Primitives>
+template <typename ExecutionSpace, typename Primitives, typename Coordinate>
 Kokkos::View<int *, typename AccessTraits<Primitives>::memory_space>
 dbscan(ExecutionSpace const &exec_space, Primitives const &primitives,
-       float eps, int core_min_size,
+       Coordinate eps, int core_min_size,
        DBSCAN::Parameters const &parameters = DBSCAN::Parameters())
 {
   Kokkos::Profiling::pushRegion("ArborX::DBSCAN");
@@ -227,8 +234,12 @@ dbscan(ExecutionSpace const &exec_space, Primitives const &primitives,
 
   using Point = typename Points::value_type;
   static_assert(GeometryTraits::is_point_v<Point>);
-  constexpr int dim = GeometryTraits::dimension_v<Point>;
-  using Box = Box<dim>;
+  constexpr int DIM = GeometryTraits::dimension_v<Point>;
+  static_assert(
+      std::is_same_v<typename GeometryTraits::coordinate_type<Point>::type,
+                     Coordinate>);
+
+  using Box = Box<DIM, Coordinate>;
 
   bool const is_special_case = (core_min_size == 2);
 
@@ -265,7 +276,7 @@ dbscan(ExecutionSpace const &exec_space, Primitives const &primitives,
       Details::HalfTraversal(
           exec_space, bvh,
           Details::FDBSCANCallback<UnionFind, CorePoints>{labels, CorePoints{}},
-          Details::WithinRadiusGetter{eps});
+          Details::WithinRadiusGetter<Coordinate>{eps});
       Kokkos::Profiling::popRegion();
     }
     else
@@ -287,7 +298,7 @@ dbscan(ExecutionSpace const &exec_space, Primitives const &primitives,
       Details::HalfTraversal(exec_space, bvh,
                              Details::FDBSCANCallback<UnionFind, CorePoints>{
                                  labels, CorePoints{num_neigh, core_min_size}},
-                             Details::WithinRadiusGetter{eps});
+                             Details::WithinRadiusGetter<Coordinate>{eps});
       Kokkos::Profiling::popRegion();
     }
   }
@@ -304,8 +315,8 @@ dbscan(ExecutionSpace const &exec_space, Primitives const &primitives,
 
     // The cell length is chosen to be eps/sqrt(dimension), so that any two
     // points within the same cell are within eps distance of each other.
-    float const h = eps / std::sqrt(dim);
-    Details::CartesianGrid<dim> const grid(bounds, h);
+    auto const h = eps / std::sqrt((Coordinate)DIM);
+    Details::CartesianGrid const grid(bounds, h);
 
     auto cell_indices = Details::computeCellIndices(exec_space, points, grid);
 
@@ -339,7 +350,7 @@ dbscan(ExecutionSpace const &exec_space, Primitives const &primitives,
     if (verbose)
     {
       printf("h = %e, n = [%zu", h, grid.extent(0));
-      for (int d = 1; d < decltype(grid)::dim; ++d)
+      for (int d = 1; d < DIM; ++d)
         printf(", %zu", grid.extent(d));
       printf("]\n");
       printf("#nonempty cells     : %10d\n", num_nonempty_cells);
@@ -409,7 +420,7 @@ dbscan(ExecutionSpace const &exec_space, Primitives const &primitives,
           permute, Kokkos::make_pair(num_points_in_dense_cells, n));
 
       auto const sparse_predicates =
-          Details::PrimitivesWithRadiusReorderedAndFiltered<
+          Details::PointsWithRadiusReorderedAndFiltered<
               Points, decltype(sparse_permute)>{points, eps, sparse_permute};
       bvh.query(exec_space, sparse_predicates,
                 Details::CountUpToN_DenseBox<MemorySpace, Points,

src/cluster/detail/ArborX_CartesianGrid.hpp

@@ -23,21 +23,19 @@
 namespace ArborX::Details
 {
 
-template <int DIM>
+template <int DIM, typename Coordinate = float>
 struct CartesianGrid
 {
 public:
-  static constexpr int dim = DIM;
-
-  CartesianGrid(Box<DIM> const &bounds, float h)
+  CartesianGrid(Box<DIM, Coordinate> const &bounds, Coordinate h)
       : _bounds(bounds)
   {
     ARBORX_ASSERT(h > 0);
     for (int d = 0; d < DIM; ++d)
       _h[d] = h;
     buildGrid();
   }
-  CartesianGrid(Box<DIM> const &bounds, float const h[DIM])
+  CartesianGrid(Box<DIM, Coordinate> const &bounds, Coordinate const h[DIM])
       : _bounds(bounds)
   {
     for (int d = 0; d < DIM; ++d)
@@ -80,7 +78,7 @@ struct CartesianGrid
       max[d] = min[d] + (i + 1) * _h[d];
       min[d] += i * _h[d];
     }
-    return Box<DIM>{min, max};
+    return Box{min, max};
   }
 
   KOKKOS_FUNCTION
@@ -127,7 +125,7 @@ struct CartesianGrid
     // run with a full NGSIM datasets, values below 3 could still produce wrong
     // results. This may still not be conservative enough, but all runs passed
     // verification when this warning was not triggered.
-    constexpr auto eps = 5 * std::numeric_limits<float>::epsilon();
+    constexpr auto eps = 5 * std::numeric_limits<Coordinate>::epsilon();
     for (int d = 0; d < DIM; ++d)
     {
       if (std::abs(_h[d] / min_corner[d]) < eps)
@@ -138,11 +136,20 @@ struct CartesianGrid
     }
   }
 
-  Box<DIM> _bounds;
-  float _h[DIM];
+  Box<DIM, Coordinate> _bounds;
+  Coordinate _h[DIM];
   size_t _n[DIM];
 };
 
+template <int DIM, typename Coordinate>
+#if KOKKOS_VERSION >= 40400
+KOKKOS_DEDUCTION_GUIDE
+#else
+KOKKOS_FUNCTION
+#endif
+    CartesianGrid(Box<DIM, Coordinate>, Coordinate)
+        -> CartesianGrid<DIM, Coordinate>;
+
 } // namespace ArborX::Details
 
 #endif

src/cluster/detail/ArborX_FDBSCANDenseBox.hpp

@@ -30,20 +30,23 @@ template <typename MemorySpace, typename Primitives, typename DenseCellOffsets,
           typename Permutation>
 struct CountUpToN_DenseBox
 {
+  using Coordinate =
+      GeometryTraits::coordinate_type_t<typename Primitives::value_type>;
+
   Kokkos::View<int *, MemorySpace> _counts;
   Primitives _primitives;
   DenseCellOffsets _dense_cell_offsets;
   int _num_dense_cells;
   Permutation _permute;
   int core_min_size;
-  float eps;
+  Coordinate eps;
   int _n;
 
   CountUpToN_DenseBox(Kokkos::View<int *, MemorySpace> const &counts,
                       Primitives const &primitives,
                       DenseCellOffsets const &dense_cell_offsets,
                       Permutation const &permute, int core_min_size_in,
-                      float eps_in, int n)
+                      Coordinate eps_in, int n)
       : _counts(counts)
       , _primitives(primitives)
       , _dense_cell_offsets(dense_cell_offsets)
@@ -93,22 +96,25 @@ template <typename UnionFind, typename CorePointsType, typename Primitives,
           typename DenseCellOffsets, typename Permutation>
 struct FDBSCANDenseBoxCallback
 {
+  using Coordinate =
+      GeometryTraits::coordinate_type_t<typename Primitives::value_type>;
+
   UnionFind _union_find;
   CorePointsType _is_core_point;
   Primitives _primitives;
   DenseCellOffsets _dense_cell_offsets;
   int _num_dense_cells;
   int _num_points_in_dense_cells;
   Permutation _permute;
-  float eps;
+  Coordinate eps;
 
   template <typename ExecutionSpace>
   FDBSCANDenseBoxCallback(UnionFind const &union_find,
                           CorePointsType const &is_core_point,
                           Primitives const &primitives,
                           DenseCellOffsets const &dense_cell_offsets,
                           ExecutionSpace const &exec_space,
-                          Permutation const &permute, float eps_in)
+                          Permutation const &permute, Coordinate eps_in)
       : _union_find(union_find)
       , _is_core_point(is_core_point)
       , _primitives(primitives)
@@ -182,11 +188,11 @@ struct FDBSCANDenseBoxCallback
 };
 
 template <typename ExecutionSpace, typename Primitives>
-Kokkos::View<size_t *, typename Primitives::memory_space>
-computeCellIndices(ExecutionSpace const &exec_space,
-                   Primitives const &primitives,
-                   CartesianGrid<GeometryTraits::dimension_v<
-                       typename Primitives::value_type>> const &grid)
+Kokkos::View<size_t *, typename Primitives::memory_space> computeCellIndices(
+    ExecutionSpace const &exec_space, Primitives const &primitives,
+    CartesianGrid<GeometryTraits::dimension_v<typename Primitives::value_type>,
+                  GeometryTraits::coordinate_type_t<
+                      typename Primitives::value_type>> const &grid)
 {
   using MemorySpace = typename Primitives::memory_space;