Skip to content

Commit

Permalink
Fix DBSCAN to use user-provided point precision
Browse files Browse the repository at this point in the history
  • Loading branch information
@aprokop
aprokop committed Jan 15, 2025
1 parent 591639e commit 561ad29
Show file tree
Hide file tree
Showing 5 changed files with 90 additions and 59 deletions.
5 changes: 3 additions & 2 deletions benchmarks/dbscan/ArborX_DBSCANVerification.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -286,9 +286,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)
{
Kokkos::Profiling::pushRegion("ArborX::DBSCAN::verify");

Expand Down
69 changes: 40 additions & 29 deletions src/cluster/ArborX_DBSCAN.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,9 +52,10 @@ struct DBSCANCorePoints
}
};

template <typename Coordinate>
struct WithinRadiusGetter
{
float _r;
Coordinate _r;

template <typename Point, typename Index>
KOKKOS_FUNCTION auto
Expand All @@ -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;
};

Expand All @@ -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;
Expand All @@ -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)
{
Expand All @@ -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);
}
};
Expand Down Expand Up @@ -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;
};
Expand Down Expand Up @@ -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");
Expand All @@ -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);

Expand Down Expand Up @@ -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
Expand All @@ -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();
}
}
Expand All @@ -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);

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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,
Expand Down
9 changes: 9 additions & 0 deletions src/cluster/detail/ArborX_CartesianGrid.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,6 +141,15 @@ struct CartesianGrid
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
24 changes: 15 additions & 9 deletions src/cluster/detail/ArborX_FDBSCANDenseBox.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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;

Expand Down
42 changes: 23 additions & 19 deletions test/tstDBSCAN.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,21 +44,22 @@ BOOST_AUTO_TEST_SUITE(DBSCAN)
BOOST_AUTO_TEST_CASE_TEMPLATE(dbscan_verifier, DeviceType, ARBORX_DEVICE_TYPES)
{
using ExecutionSpace = typename DeviceType::execution_space;
using Point = ArborX::Point<3>;
using Coordinate = float;
using Point = ArborX::Point<3, Coordinate>;
using ArborX::Details::verifyDBSCAN;

ExecutionSpace space;

{
auto points = toView<DeviceType, Point>({{{0, 0, 0}}, {{1, 1, 1}}});

auto r = std::sqrt(3);
Coordinate r = std::sqrt(3);

BOOST_TEST(verifyDBSCAN(space, points, r - 0.1, 2,
BOOST_TEST(verifyDBSCAN(space, points, r - (Coordinate)0.1, 2,
toView<DeviceType, int>({-1, -1})));
BOOST_TEST(!verifyDBSCAN(space, points, r - 0.1, 2,
BOOST_TEST(!verifyDBSCAN(space, points, r - (Coordinate)0.1, 2,
toView<DeviceType, int>({1, 2})));
BOOST_TEST(!verifyDBSCAN(space, points, r - 0.1, 2,
BOOST_TEST(!verifyDBSCAN(space, points, r - (Coordinate)0.1, 2,
toView<DeviceType, int>({1, 1})));
BOOST_TEST(
verifyDBSCAN(space, points, r, 2, toView<DeviceType, int>({1, 1})));
Expand All @@ -74,9 +75,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(dbscan_verifier, DeviceType, ARBORX_DEVICE_TYPES)
auto points = toView<DeviceType, Point>(
{{{0, 0, 0}}, {{1, 1, 1}}, {{3, 3, 3}}, {{6, 6, 6}}});

auto r = std::sqrt(3.f);
auto r2 = std::sqrt(12.f);
auto r3 = std::sqrt(48.f);
Coordinate r = std::sqrt(3);
Coordinate r2 = std::sqrt(12);
Coordinate r3 = std::sqrt(48);

BOOST_TEST(verifyDBSCAN(space, points, r, 2,
toView<DeviceType, int>({1, 1, -1, -1})));
Expand Down Expand Up @@ -154,9 +155,10 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(dbscan_verifier, DeviceType, ARBORX_DEVICE_TYPES)
BOOST_AUTO_TEST_CASE_TEMPLATE(dbscan, DeviceType, ARBORX_DEVICE_TYPES)
{
using ExecutionSpace = typename DeviceType::execution_space;
using Coordinate = float;
using ArborX::dbscan;
using ArborX::Details::verifyDBSCAN;
using Point = ArborX::Point<3>;
using Point = ArborX::Point<3, Coordinate>;

ExecutionSpace space;

Expand All @@ -168,20 +170,21 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(dbscan, DeviceType, ARBORX_DEVICE_TYPES)
{
auto points = toView<DeviceType, Point>({{{0, 0, 0}}, {{1, 1, 1}}});

auto r = std::sqrt(3.1f);
Coordinate r = std::sqrt(3.1);

BOOST_TEST(verifyDBSCAN(space, points, r - 0.1, 2,
dbscan(space, points, r - 0.1, 2, params)));
BOOST_TEST(
verifyDBSCAN(space, points, r - (Coordinate)0.1, 2,
dbscan(space, points, r - (Coordinate)0.1, 2, params)));
BOOST_TEST(verifyDBSCAN(space, points, r, 2,
dbscan(space, points, r, 2, params)));
BOOST_TEST(verifyDBSCAN(space, points, r, 3,
dbscan(space, points, r, 3, params)));

// Test non-View primitives
HiddenView<decltype(points)> hidden_points{points};
BOOST_TEST(
verifyDBSCAN(space, hidden_points, r - 0.1, 2,
dbscan(space, hidden_points, r - 0.1, 2, params)));
BOOST_TEST(verifyDBSCAN(
space, hidden_points, r - (Coordinate)0.1, 2,
dbscan(space, hidden_points, r - (Coordinate)0.1, 2, params)));
BOOST_TEST(verifyDBSCAN(space, hidden_points, r, 2,
dbscan(space, hidden_points, r, 2, params)));
BOOST_TEST(verifyDBSCAN(space, hidden_points, r, 3,
Expand All @@ -192,7 +195,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(dbscan, DeviceType, ARBORX_DEVICE_TYPES)
auto points = toView<DeviceType, Point>(
{{{0, 0, 0}}, {{1, 1, 1}}, {{3, 3, 3}}, {{6, 6, 6}}});

auto r = std::sqrt(3.1f);
Coordinate r = std::sqrt(3.1);

BOOST_TEST(verifyDBSCAN(space, points, r, 2,
dbscan(space, points, r, 2, params)));
Expand Down Expand Up @@ -227,9 +230,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(dbscan, DeviceType, ARBORX_DEVICE_TYPES)
{{1, -0.5, 0}}});

BOOST_TEST(verifyDBSCAN(space, points, 1.0, 3,
dbscan(space, points, 1, 3, params)));
dbscan(space, points, (Coordinate)1, 3, params)));
BOOST_TEST(verifyDBSCAN(space, points, 1.0, 4,
dbscan(space, points, 1, 4, params)));
dbscan(space, points, (Coordinate)1, 4, params)));
}
}

Expand All @@ -256,7 +259,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(dbscan, DeviceType, ARBORX_DEVICE_TYPES)

// This does *not* guarantee to trigger the issue, as it depends on the
// specific implementation and runtime. But it may.
BOOST_TEST(verifyDBSCAN(space, points, 1, 4, dbscan(space, points, 1, 4)));
BOOST_TEST(verifyDBSCAN(space, points, (Coordinate)1, 4,
dbscan(space, points, (Coordinate)1, 4)));
}
}

Expand Down

0 comments on commit 561ad29

Please sign in to comment.