From c0a1b8b470dd9c2a0232dbebb7ac4ff3c19727cb Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Tue, 20 Apr 2021 17:17:53 -0400 Subject: [PATCH 001/177] HDBSCAN --- cpp/include/cuml/cluster/hdbscan.hpp | 30 +++++ cpp/src/hdbscan/hdbscan.cu | 36 ++++++ cpp/src/hdbscan/reachability.cuh | 130 ++++++++++++++++++++ cpp/src/hdbscan/runner.h | 110 +++++++++++++++++ cpp/src/hdbscan/tree.cuh | 172 +++++++++++++++++++++++++++ cpp/test/sg/hdbscan_test.cu | 0 6 files changed, 478 insertions(+) create mode 100644 cpp/include/cuml/cluster/hdbscan.hpp create mode 100644 cpp/src/hdbscan/hdbscan.cu create mode 100644 cpp/src/hdbscan/reachability.cuh create mode 100644 cpp/src/hdbscan/runner.h create mode 100644 cpp/src/hdbscan/tree.cuh create mode 100644 cpp/test/sg/hdbscan_test.cu diff --git a/cpp/include/cuml/cluster/hdbscan.hpp b/cpp/include/cuml/cluster/hdbscan.hpp new file mode 100644 index 0000000000..3a3c61ca42 --- /dev/null +++ b/cpp/include/cuml/cluster/hdbscan.hpp @@ -0,0 +1,30 @@ +/* + * Copyright (c) 2018-2021, 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 + +namespace ML { + +template +void hdbscan(const raft::handle_t &handle, value_t *X, size_t m, size_t n, + raft::distance::DistanceType metric, int k, int min_pts, + float alpha, hdbscan_output *out); +}; // end namespace ML \ No newline at end of file diff --git a/cpp/src/hdbscan/hdbscan.cu b/cpp/src/hdbscan/hdbscan.cu new file mode 100644 index 0000000000..b7df85d8b6 --- /dev/null +++ b/cpp/src/hdbscan/hdbscan.cu @@ -0,0 +1,36 @@ +/* + * Copyright (c) 2021, NVIDIA CORPORATION. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include + +#include + +namespace ML { + +template +void hdbscan(const raft::handle_t &handle, value_t *X, size_t m, size_t n, + raft::distance::DistanceType metric, int k, int min_pts, + float alpha, hdbscan_output *out) { + HDBSCAN::_fit(handle, X, m, n, metric, k, min_pts, alpha); +} + +void hdbscan(const raft::handle_t &handle, const float *X, size_t m, size_t n, + raft::distance::DistanceType metric, int k, int min_pts, + float alpha, hdbscan_output *out); + +}; // end namespace ML \ No newline at end of file diff --git a/cpp/src/hdbscan/reachability.cuh b/cpp/src/hdbscan/reachability.cuh new file mode 100644 index 0000000000..f615c2d4e7 --- /dev/null +++ b/cpp/src/hdbscan/reachability.cuh @@ -0,0 +1,130 @@ +/* + * Copyright (c) 2018-2020, 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 + +namespace ML { +namespace HDBSCAN { +namespace Reachability { + +template +__global__ void core_distances_kernel(value_t *knn_dists, int min_pts, + value_t *out) { + int row = blockIdx.x * blockDim.x + threadIdx.x; + out[row] = knn_dists[row + min_pts]; +} + +/** + * Extract core distances from KNN graph. This is essentially + * performing a knn_dists[:,min_pts] + * @tparam value_idx data type for integrals + * @tparam value_t data type for distance + * @tparam tpb block size for kernel + * @param knn_dists knn distance array + * @param min_pts + * @param n + * @param out + * @param stream + */ +template +void core_distances(value_t *knn_dists, int min_pts, size_t n, value_t *out, + cudaStream_t stream) { + int blocks = raft::ceildiv(n * min_pts, (size_t)tpb); + core_distances_kernel + <<>>(knn_dists, min_pts, out); +} + +template +__global__ void mutual_reachability_kernel(value_t *pw_dists, + value_t *core_dists, size_t m, + size_t n, value_t *out) { + int idx = blockIdx.x * blockDim.x + threadIdx.x; + int row = idx / n; + int col = idx % n; + + value_idx a = row; + value_idx b = col; + + value_t core_a = core_dists[a]; + value_t core_b = core_dists[b]; + + value_t dist = pw_dists[row + col]; + + out[row + col] = max(core_a, core_b, dist); +} + +template +void mutual_reachability(value_t *pw_dists, value_t *core_dists, size_t m, + cudaStream_t stream) { + int blocks = raft::ceildiv(m * m, (size_t)tpb); + + mutual_reachability_kernel + <<>>(pw_dists, core_dists, m, m, stream); +} + +/** + * Constructs a mutual reachability graph, which is a k-nearest neighbors + * graph projected into mutual reachability space using the following + * function for each data point, where core_distance is the distance + * to the kth neighbor: max(core_distance(a), core_distance(b), d(a, b)) + * + * @tparam value_idx + * @tparam value_t + * @param[in] handle + * @param[in] X + * @param[in] m + * @param[in] n + * @param[in] metric + * @param[out] pw_dists + * @param[in] k + */ +template +void mutual_reachability_dists(const raft::handle_t &handle, + const value_t *X, size_t m, size_t n, + raft::distance::DistanceType metric, + int k, value_idx *inds, value_t *dists, + value_t *core_dists) { + auto stream = handle.get_stream(); + + // perform knn + brute_force_knn(handle, {X}, {m}, n, X, m, inds, dists, k, true, + true, metric); + + // Slice core distances (distances to kth nearest neighbor) + core_distances(dists, k, m, core_dists, stream); + + // Project into mutual reachability space. + // Note that it's not guaranteed the knn graph will be connected + // at this point so the core distances will need to be returned + // so additional points can be added to the graph and projected + // ito mutual reachability space later. + mutual_reachability(inds, dists, + core_dists, m, n, stream); +} + +}; // end namespace Reachability +}; // end namespace HDBSCAN +}; // end namespace ML \ No newline at end of file diff --git a/cpp/src/hdbscan/runner.h b/cpp/src/hdbscan/runner.h new file mode 100644 index 0000000000..2b1e94dc63 --- /dev/null +++ b/cpp/src/hdbscan/runner.h @@ -0,0 +1,110 @@ +/* + * Copyright (c) 2021, 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 "reachability.cuh" + +namespace ML { +namespace HDBSCAN { + +template +struct MSTEpilogueReachability { + + MSTEpilogueReachability(value_idx m_, value_t *core_distances_): + core_distances(core_distances_), m(m_){} + + void operator()(raft::handle_t &handle, value_idx *coo_rows, + value_idx *coo_cols, value_t *coo_data) { + // TODO: Schedule kernel that uses the core distances + // to perform the max(core_dist(src), core_dist(dst), d(src, dst)) operation + } + + private: + value_t *core_distances; + value_idx m; +} + + +template +void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, + raft::distance::DistanceType metric, int k, int min_pts, + float alpha) { + auto d_alloc = handle.get_device_allocator(); + auto stream = handle.get_stream(); + + /** + * Mutual reachability graph + */ + + rmm::device_uvector mutual_reachability_graph_inds(k * m, stream); + rmm::device_uvector mutual_reachability_graph_dists(k * m, stream); + rmm::device_uvector core_dists(k * m, stream); + + Reachability::mutual_reachability_dists( + handle, X, m, n, metric, min_pts, k, mutual_reachability_graph_inds.data(), + mutual_reachability_graph_dists.data(), core_dists.data()); + + /** + * Construct MST sorted by weights + */ + rmm::device_uvector mst_rows(m - 1, stream); + rmm::device_uvector mst_cols(m - 1, stream); + rmm::device_uvector mst_data(m - 1, stream); + + // during knn graph connection + raft::hierarchy::detail::build_sorted_mst(handle, X, mutual_reachability_graph_inds.data(), + mutual_reachability_graph_dists.data(), + m, n, + mst_rows, mst_cols, mst_data, + k * m, metric, 10, + MSTEpilogueReachability()); + + /** + * Perform hierarchical labeling + */ + value_idx n_edges = m - 1; + + rmm::device_uvector out_src(n_edges, stream); + rmm::device_uvector out_dst(n_edges, stream); + rmm::device_uvector out_delta(n_edges, stream); + rmm::device_uvector out_size(n_edges, stream); + + raft::hierarchy::detail::build_dendrogram_host( + mst_rows.data(), mst_cols.data(), mst_data.data(), n_edges, out_src.data(), + out_dst.data(), out_delta.data(), out_size.data()); + + /** + * Condense branches of tree according to min cluster size + */ + + + + /** + * Extract labels from stability + */ +} + +}; // end namespace HDBSCAN +}; // end namespace ML \ No newline at end of file diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh new file mode 100644 index 0000000000..a1e52d1a30 --- /dev/null +++ b/cpp/src/hdbscan/tree.cuh @@ -0,0 +1,172 @@ +/* + * Copyright (c) 2021, 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 ML { +namespace HDBSCAN { +namespace Tree { + +template +__device__ value_t get_lambda(value_idx node, value_idx num_points, value_t *deltas) { + + value_t delta = deltas[node - num_points]; + if(delta > 0.0) + return 1.0 / delta; + return std::numeric_limits::max(); +} + +/** + * + * @tparam value_idx + * @tparam value_t + * @param frontier + * @param ignore Should be initialized to -1 + * @param next_label + * @param relabel + * @param hierarchy + * @param deltas + * @param sizes + * @param n_leaves + * @param num_points + * @param min_cluster_size + */ +template +__global__ void condense_hierarchy_kernel(bool *frontier, value_idx *ignore, value_idx *next_label, + value_idx *relabel, value_idx *hierarchy, + value_t *deltas, value_idx *sizes, + int n_leaves, int num_points, int min_cluster_size) { + + int node = blockDim.x * blockIdx.x + threadIdx.x; + + // If node is in frontier, flip frontier for children + if(node <= n_leaves * 2 && frontier[node]) { + frontier[node] = false; + + // TODO: Check bounds + value_idx left_child = hierarchy[(node - num_points) * 2]; + value_idx right_child = hierarchy[((node - num_points) * 2) + 1]; + + frontier[left_child] = true; + frontier[right_child] = true; + + bool ignore_val = ignore[node]; + bool should_ignore = ignore_val > -1; + + // If the current node is being ignored (e.g. > -1) then propagate the ignore + // to children, if any + ignore[left_child] = (should_ignore * ignore_val) + (!should_ignore * -1); + ignore[right_child] = (should_ignore * ignore_val) + (!should_ignore * -1); + + if (node < num_points) { + // TODO: append relabel[should_ignore], node, get_lambda(should_ignore), 1 + + } + + // If node is not ignored and is not a leaf, condense its children + // if necessary + else if (!should_ignore and node >= num_points) { + value_idx left_child = hierarchy[(node - num_points) * 2]; + value_idx right_child = hierarchy[((node - num_points) * 2) + 1]; + + value_t lambda_value = get_lambda(node, num_points, deltas); + + // TODO: Convert to boolean arithmetic + int left_count = + left_child >= num_points ? sizes[left_child - num_points] : 1; + int right_count = + right_child >= num_points ? sizes[right_child - num_points] : 1; + + // If both children are large enough, they should be relabeled and + // included directly in the output hierarchy. + if (left_count >= min_cluster_size && right_count >= min_cluster_size) { + relabel[left_child] = atomicAdd(next_label, 1); + // TODO: Output new hierarchy entry for: relabel[node], relabel[left], lambda_value, left_count + + relabel[right_child] = atomicAdd(next_label, 1); + // TODO Output new hierarchy entry for: relabel[node], relabel[right], lambda_value, right_count + } + + // Consume left or right child as necessary + bool left_child_too_small = left_count < min_cluster_size; + bool right_child_too_small = right_count < min_cluster_size; + ignore[left_child] = + (left_child_too_small * node) + (!left_child_too_small * -1); + ignore[right_child] = + (right_child_too_small * node) + (!right_child_too_small * -1); + + // If only left or right child is too small, consume it and relabel the other + // (to it can be its own cluster) + bool only_left_child_too_small = + left_child_too_small && !right_child_too_small; + bool only_right_child_too_small = + !left_child_too_small && right_child_too_small; + + relabel[right_child] = (only_left_child_too_small * relabel[node]) + + (!only_left_child_too_small * -1); + relabel[left_child] = (only_right_child_too_small * relabel[node]) + + (!only_right_child_too_small * -1); + } + } +} + +template +void condense_hierarchy(raft::handle_t &handle, value_idx *children, value_t *delta, + value_idx *sizes, + int min_pts, int n_leaves) { + + rmm::device_uvector frontier(n_leaves*2, handle.get_stream()); + rmm::device_uvector ignore(n_leaves*2, handle.get_stream()); + + int root = 2 * n_leaves; + int num_points = floor(root / 2.0) + 1; + + rmm::device_uvector relabel(root+1, handle.get_stream()); + + // TODO: Set this properly on device + relabel[root] = num_points; + + // While frontier is not empty, perform single bfs through tree + size_t grid = raft::ceildiv(n_leaves * 2, (size_t)tpb); + + value_idx n_elements_to_traverse = thrust::reduce(thrust::cuda::par.on(handle.get_stream()), + frontier.data(), frontier.data()+(n_leaves *2), 0); + + rmm::device_uvector next_label(1, handle.get_stream()); + // TODO: Set this properly on device + next_label[0] = num_points + 1; + + while(n_elements_to_traverse > 0) { + condense_hierarchy_kernel<<>>( + frontier.data(), ignore.data(), next_label.data(), relabel.data(), children, delta, + sizes, n_leaves, num_points, min_cluster_size); + + n_elements_to_traverse = thrust::reduce(thrust::cuda::par.on(handle.get_stream()), + frontier.data(), frontier.data()+(n_leaves *2), 0);; + + CUDA_CHECK(cudaStreamSynchronize(handle.get_stream())); + } +} +}; // end namespace Tree +}; // end namespace HDBSCAN +}; // end namespace ML \ No newline at end of file diff --git a/cpp/test/sg/hdbscan_test.cu b/cpp/test/sg/hdbscan_test.cu new file mode 100644 index 0000000000..e69de29bb2 From 3800f76869bcae85cae60f5c6b2e71fc3e46fabe Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Tue, 20 Apr 2021 17:19:08 -0400 Subject: [PATCH 002/177] Fixing style --- cpp/src/hdbscan/reachability.cuh | 16 +++++----- cpp/src/hdbscan/runner.h | 20 +++++------- cpp/src/hdbscan/tree.cuh | 54 ++++++++++++++++---------------- 3 files changed, 42 insertions(+), 48 deletions(-) diff --git a/cpp/src/hdbscan/reachability.cuh b/cpp/src/hdbscan/reachability.cuh index f615c2d4e7..592686530c 100644 --- a/cpp/src/hdbscan/reachability.cuh +++ b/cpp/src/hdbscan/reachability.cuh @@ -102,16 +102,16 @@ void mutual_reachability(value_t *pw_dists, value_t *core_dists, size_t m, * @param[in] k */ template -void mutual_reachability_dists(const raft::handle_t &handle, - const value_t *X, size_t m, size_t n, - raft::distance::DistanceType metric, - int k, value_idx *inds, value_t *dists, +void mutual_reachability_dists(const raft::handle_t &handle, const value_t *X, + size_t m, size_t n, + raft::distance::DistanceType metric, int k, + value_idx *inds, value_t *dists, value_t *core_dists) { auto stream = handle.get_stream(); // perform knn - brute_force_knn(handle, {X}, {m}, n, X, m, inds, dists, k, true, - true, metric); + brute_force_knn(handle, {X}, {m}, n, X, m, inds, dists, k, true, true, + metric); // Slice core distances (distances to kth nearest neighbor) core_distances(dists, k, m, core_dists, stream); @@ -121,8 +121,8 @@ void mutual_reachability_dists(const raft::handle_t &handle, // at this point so the core distances will need to be returned // so additional points can be added to the graph and projected // ito mutual reachability space later. - mutual_reachability(inds, dists, - core_dists, m, n, stream); + mutual_reachability(inds, dists, core_dists, m, n, + stream); } }; // end namespace Reachability diff --git a/cpp/src/hdbscan/runner.h b/cpp/src/hdbscan/runner.h index 2b1e94dc63..5853d4c6a4 100644 --- a/cpp/src/hdbscan/runner.h +++ b/cpp/src/hdbscan/runner.h @@ -29,11 +29,10 @@ namespace ML { namespace HDBSCAN { -template +template struct MSTEpilogueReachability { - - MSTEpilogueReachability(value_idx m_, value_t *core_distances_): - core_distances(core_distances_), m(m_){} + MSTEpilogueReachability(value_idx m_, value_t *core_distances_) + : core_distances(core_distances_), m(m_) {} void operator()(raft::handle_t &handle, value_idx *coo_rows, value_idx *coo_cols, value_t *coo_data) { @@ -46,7 +45,6 @@ struct MSTEpilogueReachability { value_idx m; } - template void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, raft::distance::DistanceType metric, int k, int min_pts, @@ -74,12 +72,10 @@ void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, rmm::device_uvector mst_data(m - 1, stream); // during knn graph connection - raft::hierarchy::detail::build_sorted_mst(handle, X, mutual_reachability_graph_inds.data(), - mutual_reachability_graph_dists.data(), - m, n, - mst_rows, mst_cols, mst_data, - k * m, metric, 10, - MSTEpilogueReachability()); + raft::hierarchy::detail::build_sorted_mst( + handle, X, mutual_reachability_graph_inds.data(), + mutual_reachability_graph_dists.data(), m, n, mst_rows, mst_cols, mst_data, + k * m, metric, 10, MSTEpilogueReachability()); /** * Perform hierarchical labeling @@ -99,8 +95,6 @@ void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, * Condense branches of tree according to min cluster size */ - - /** * Extract labels from stability */ diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index a1e52d1a30..12e45c1937 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -27,12 +27,11 @@ namespace ML { namespace HDBSCAN { namespace Tree { -template -__device__ value_t get_lambda(value_idx node, value_idx num_points, value_t *deltas) { - +template +__device__ value_t get_lambda(value_idx node, value_idx num_points, + value_t *deltas) { value_t delta = deltas[node - num_points]; - if(delta > 0.0) - return 1.0 / delta; + if (delta > 0.0) return 1.0 / delta; return std::numeric_limits::max(); } @@ -51,16 +50,15 @@ __device__ value_t get_lambda(value_idx node, value_idx num_points, value_t *del * @param num_points * @param min_cluster_size */ -template -__global__ void condense_hierarchy_kernel(bool *frontier, value_idx *ignore, value_idx *next_label, - value_idx *relabel, value_idx *hierarchy, - value_t *deltas, value_idx *sizes, - int n_leaves, int num_points, int min_cluster_size) { - +template +__global__ void condense_hierarchy_kernel( + bool *frontier, value_idx *ignore, value_idx *next_label, value_idx *relabel, + value_idx *hierarchy, value_t *deltas, value_idx *sizes, int n_leaves, + int num_points, int min_cluster_size) { int node = blockDim.x * blockIdx.x + threadIdx.x; // If node is in frontier, flip frontier for children - if(node <= n_leaves * 2 && frontier[node]) { + if (node <= n_leaves * 2 && frontier[node]) { frontier[node] = false; // TODO: Check bounds @@ -130,18 +128,17 @@ __global__ void condense_hierarchy_kernel(bool *frontier, value_idx *ignore, val } } -template -void condense_hierarchy(raft::handle_t &handle, value_idx *children, value_t *delta, - value_idx *sizes, - int min_pts, int n_leaves) { - - rmm::device_uvector frontier(n_leaves*2, handle.get_stream()); - rmm::device_uvector ignore(n_leaves*2, handle.get_stream()); +template +void condense_hierarchy(raft::handle_t &handle, value_idx *children, + value_t *delta, value_idx *sizes, int min_pts, + int n_leaves) { + rmm::device_uvector frontier(n_leaves * 2, handle.get_stream()); + rmm::device_uvector ignore(n_leaves * 2, handle.get_stream()); int root = 2 * n_leaves; int num_points = floor(root / 2.0) + 1; - rmm::device_uvector relabel(root+1, handle.get_stream()); + rmm::device_uvector relabel(root + 1, handle.get_stream()); // TODO: Set this properly on device relabel[root] = num_points; @@ -149,20 +146,23 @@ void condense_hierarchy(raft::handle_t &handle, value_idx *children, value_t *de // While frontier is not empty, perform single bfs through tree size_t grid = raft::ceildiv(n_leaves * 2, (size_t)tpb); - value_idx n_elements_to_traverse = thrust::reduce(thrust::cuda::par.on(handle.get_stream()), - frontier.data(), frontier.data()+(n_leaves *2), 0); + value_idx n_elements_to_traverse = + thrust::reduce(thrust::cuda::par.on(handle.get_stream()), frontier.data(), + frontier.data() + (n_leaves * 2), 0); rmm::device_uvector next_label(1, handle.get_stream()); // TODO: Set this properly on device next_label[0] = num_points + 1; - while(n_elements_to_traverse > 0) { + while (n_elements_to_traverse > 0) { condense_hierarchy_kernel<<>>( - frontier.data(), ignore.data(), next_label.data(), relabel.data(), children, delta, - sizes, n_leaves, num_points, min_cluster_size); + frontier.data(), ignore.data(), next_label.data(), relabel.data(), + children, delta, sizes, n_leaves, num_points, min_cluster_size); - n_elements_to_traverse = thrust::reduce(thrust::cuda::par.on(handle.get_stream()), - frontier.data(), frontier.data()+(n_leaves *2), 0);; + n_elements_to_traverse = + thrust::reduce(thrust::cuda::par.on(handle.get_stream()), frontier.data(), + frontier.data() + (n_leaves * 2), 0); + ; CUDA_CHECK(cudaStreamSynchronize(handle.get_stream())); } From 362d429bb07e8dd13fd224af4cb933f346fcf53b Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Wed, 21 Apr 2021 09:40:57 -0400 Subject: [PATCH 003/177] Progress --- cpp/src/hdbscan/tree.cuh | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index 12e45c1937..1a46eb54cd 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -52,7 +52,7 @@ __device__ value_t get_lambda(value_idx node, value_idx num_points, */ template __global__ void condense_hierarchy_kernel( - bool *frontier, value_idx *ignore, value_idx *next_label, value_idx *relabel, + bool *frontier, value_idx *ignore, value_idx *relabel, value_idx *hierarchy, value_t *deltas, value_idx *sizes, int n_leaves, int num_points, int min_cluster_size) { int node = blockDim.x * blockIdx.x + threadIdx.x; @@ -98,10 +98,10 @@ __global__ void condense_hierarchy_kernel( // If both children are large enough, they should be relabeled and // included directly in the output hierarchy. if (left_count >= min_cluster_size && right_count >= min_cluster_size) { - relabel[left_child] = atomicAdd(next_label, 1); + relabel[left_child] = node; // TODO: Output new hierarchy entry for: relabel[node], relabel[left], lambda_value, left_count - relabel[right_child] = atomicAdd(next_label, 1); + relabel[right_child] = node; // TODO Output new hierarchy entry for: relabel[node], relabel[right], lambda_value, right_count } @@ -141,7 +141,7 @@ void condense_hierarchy(raft::handle_t &handle, value_idx *children, rmm::device_uvector relabel(root + 1, handle.get_stream()); // TODO: Set this properly on device - relabel[root] = num_points; + relabel[root] = root; // While frontier is not empty, perform single bfs through tree size_t grid = raft::ceildiv(n_leaves * 2, (size_t)tpb); @@ -150,10 +150,6 @@ void condense_hierarchy(raft::handle_t &handle, value_idx *children, thrust::reduce(thrust::cuda::par.on(handle.get_stream()), frontier.data(), frontier.data() + (n_leaves * 2), 0); - rmm::device_uvector next_label(1, handle.get_stream()); - // TODO: Set this properly on device - next_label[0] = num_points + 1; - while (n_elements_to_traverse > 0) { condense_hierarchy_kernel<<>>( frontier.data(), ignore.data(), next_label.data(), relabel.data(), @@ -162,10 +158,25 @@ void condense_hierarchy(raft::handle_t &handle, value_idx *children, n_elements_to_traverse = thrust::reduce(thrust::cuda::par.on(handle.get_stream()), frontier.data(), frontier.data() + (n_leaves * 2), 0); - ; CUDA_CHECK(cudaStreamSynchronize(handle.get_stream())); } + + // TODO: Normalize labels so they are drawn from a monotonically increasing set. +} + +template +void compute_stabilities(value_idx *condensed_hierarchy, value_idx *lambdas, value_idx *sizes, + int n_leaves) { + + // TODO: Reverse topological sort (e.g. sort hierarchy, lambdas, and sizes by lambda) + + // TODO: Perform single loop through topologically sorted condensed hierarchy children, + // marking birth lambdas + + // TODO: Perform loop through condensed hierarchy in parallel, building array of stabilities for each cluster + + // TODO: Compute sizes of each } }; // end namespace Tree }; // end namespace HDBSCAN From abf226760f07eec8cb53772fdd3779c409a63959 Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Thu, 22 Apr 2021 12:34:41 -0400 Subject: [PATCH 004/177] Checking in --- cpp/src/hdbscan/tree.cuh | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index 1a46eb54cd..bcd44c5689 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -171,12 +171,9 @@ void compute_stabilities(value_idx *condensed_hierarchy, value_idx *lambdas, val // TODO: Reverse topological sort (e.g. sort hierarchy, lambdas, and sizes by lambda) - // TODO: Perform single loop through topologically sorted condensed hierarchy children, - // marking birth lambdas + // TODO: Segmented reduction on min_lambda within each cluster - // TODO: Perform loop through condensed hierarchy in parallel, building array of stabilities for each cluster - - // TODO: Compute sizes of each + // TODO: Embarassingly parallel construction of output } }; // end namespace Tree }; // end namespace HDBSCAN From c68bd8dca0a7ccc6e214091e8ec659983b394d74 Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Thu, 22 Apr 2021 14:17:18 -0400 Subject: [PATCH 005/177] Checking in todo stubs of remaining work for initial implementation --- cpp/src/hdbscan/tree.cuh | 46 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 44 insertions(+), 2 deletions(-) diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index bcd44c5689..5b0bb12fff 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -99,10 +99,12 @@ __global__ void condense_hierarchy_kernel( // included directly in the output hierarchy. if (left_count >= min_cluster_size && right_count >= min_cluster_size) { relabel[left_child] = node; - // TODO: Output new hierarchy entry for: relabel[node], relabel[left], lambda_value, left_count + // TODO: Output new hierarchy entry for: relabel[node], + // relabel[left], lambda_value, left_count relabel[right_child] = node; - // TODO Output new hierarchy entry for: relabel[node], relabel[right], lambda_value, right_count + // TODO Output new hierarchy entry for: relabel[node], + // relabel[right], lambda_value, right_count } // Consume left or right child as necessary @@ -175,6 +177,46 @@ void compute_stabilities(value_idx *condensed_hierarchy, value_idx *lambdas, val // TODO: Embarassingly parallel construction of output } + +template +void excess_of_mass() { + + // TODO: Build CSR of cluster tree with stabilities of each child as the weights + + // TODO: Segmented reduction over CSR of cluster tree + + // TODO: Perform bfs, starting at root- + // TODO: Maintain frontier and is_cluster array. + // TODO: In each iteration, children are added to tree + // TODO: If node has is_cluster[node] = false, set children to false + // TODO: else subtree stability > stability[node] or cluster_sizes[node] > max_cluster_size + // TODO: set is_cluster[node] = false and stability[node] = subtree_stability +} + +template +void get_stability_scores() { + + // TODO: Perform segmented reduction to compute cluster_size + + // TODO: Embarassingly parallel +} + +template +void do_labelling() { + + // TODO: Similar to SLHC dendrogram construction, this one is probably best done + // on host, at least for the first iteration +} + + +template +void get_probabilities() { + + // TODO: Compute deaths array similarly to compute_stabilities + + // TODO: Embarassingly parallel +} + }; // end namespace Tree }; // end namespace HDBSCAN }; // end namespace ML \ No newline at end of file From 78a8432c9ab452103078b9a4de9b68a95052ad2b Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Thu, 22 Apr 2021 15:12:44 -0400 Subject: [PATCH 006/177] Creating output condensed hierarchy --- cpp/src/hdbscan/tree.cuh | 249 +++++++++++++++++++++++++++------------ 1 file changed, 174 insertions(+), 75 deletions(-) diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index 5b0bb12fff..8bf3f48a7e 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -53,97 +53,192 @@ __device__ value_t get_lambda(value_idx node, value_idx num_points, template __global__ void condense_hierarchy_kernel( bool *frontier, value_idx *ignore, value_idx *relabel, - value_idx *hierarchy, value_t *deltas, value_idx *sizes, int n_leaves, - int num_points, int min_cluster_size) { + const value_idx *hierarchy, const value_t *deltas, + const value_idx *sizes, int n_leaves, + int num_points, int min_cluster_size, + value_idx *out_parent, value_idx *out_child, + value_t *out_lambda, value_idx *out_count) { + int node = blockDim.x * blockIdx.x + threadIdx.x; // If node is in frontier, flip frontier for children - if (node <= n_leaves * 2 && frontier[node]) { - frontier[node] = false; + if(node > n_leaves * 2 || !frontier[node]) + return; - // TODO: Check bounds - value_idx left_child = hierarchy[(node - num_points) * 2]; - value_idx right_child = hierarchy[((node - num_points) * 2) + 1]; + frontier[node] = false; - frontier[left_child] = true; - frontier[right_child] = true; + // TODO: Check bounds + value_idx left_child = hierarchy[(node - num_points) * 2]; + value_idx right_child = hierarchy[((node - num_points) * 2) + 1]; - bool ignore_val = ignore[node]; - bool should_ignore = ignore_val > -1; + frontier[left_child] = true; + frontier[right_child] = true; - // If the current node is being ignored (e.g. > -1) then propagate the ignore - // to children, if any - ignore[left_child] = (should_ignore * ignore_val) + (!should_ignore * -1); - ignore[right_child] = (should_ignore * ignore_val) + (!should_ignore * -1); + bool ignore_val = ignore[node]; + bool should_ignore = ignore_val > -1; - if (node < num_points) { - // TODO: append relabel[should_ignore], node, get_lambda(should_ignore), 1 + // If the current node is being ignored (e.g. > -1) then propagate the ignore + // to children, if any + ignore[left_child] = (should_ignore * ignore_val) + (!should_ignore * -1); + ignore[right_child] = (should_ignore * ignore_val) + (!should_ignore * -1); - } + if (node < num_points) { + out_parent[node] = relabel[should_ignore]; + out_child[node] = node; + out_lambda[node] = get_lambda(should_ignore, num_points, deltas); + out_count[node] = 1; + } + + // If node is not ignored and is not a leaf, condense its children + // if necessary + else if (!should_ignore and node >= num_points) { + value_idx left_child = hierarchy[(node - num_points) * 2]; + value_idx right_child = hierarchy[((node - num_points) * 2) + 1]; - // If node is not ignored and is not a leaf, condense its children - // if necessary - else if (!should_ignore and node >= num_points) { - value_idx left_child = hierarchy[(node - num_points) * 2]; - value_idx right_child = hierarchy[((node - num_points) * 2) + 1]; - - value_t lambda_value = get_lambda(node, num_points, deltas); - - // TODO: Convert to boolean arithmetic - int left_count = - left_child >= num_points ? sizes[left_child - num_points] : 1; - int right_count = - right_child >= num_points ? sizes[right_child - num_points] : 1; - - // If both children are large enough, they should be relabeled and - // included directly in the output hierarchy. - if (left_count >= min_cluster_size && right_count >= min_cluster_size) { - relabel[left_child] = node; - // TODO: Output new hierarchy entry for: relabel[node], - // relabel[left], lambda_value, left_count - - relabel[right_child] = node; - // TODO Output new hierarchy entry for: relabel[node], - // relabel[right], lambda_value, right_count - } - - // Consume left or right child as necessary - bool left_child_too_small = left_count < min_cluster_size; - bool right_child_too_small = right_count < min_cluster_size; - ignore[left_child] = - (left_child_too_small * node) + (!left_child_too_small * -1); - ignore[right_child] = - (right_child_too_small * node) + (!right_child_too_small * -1); - - // If only left or right child is too small, consume it and relabel the other - // (to it can be its own cluster) - bool only_left_child_too_small = - left_child_too_small && !right_child_too_small; - bool only_right_child_too_small = - !left_child_too_small && right_child_too_small; - - relabel[right_child] = (only_left_child_too_small * relabel[node]) + - (!only_left_child_too_small * -1); - relabel[left_child] = (only_right_child_too_small * relabel[node]) + - (!only_right_child_too_small * -1); + value_t lambda_value = get_lambda(node, num_points, deltas); + + int left_count = + left_child >= num_points ? sizes[left_child - num_points] : 1; + int right_count = + right_child >= num_points ? sizes[right_child - num_points] : 1; + + // If both children are large enough, they should be relabeled and + // included directly in the output hierarchy. + if (left_count >= min_cluster_size && right_count >= min_cluster_size) { + relabel[left_child] = node; + out_parent[node] = relabel[node]; + out_child[node] = relabel[left_child]; + out_lambda[node] = lambda_value; + out_count[node] = left_count; + + relabel[right_child] = node; + out_parent[node] = relabel[node]; + out_child[node] = relabel[right_child]; + out_lambda[node] = lambda_value; + out_count[node] = left_count; } + + // Consume left or right child as necessary + bool left_child_too_small = left_count < min_cluster_size; + bool right_child_too_small = right_count < min_cluster_size; + ignore[left_child] = + (left_child_too_small * node) + (!left_child_too_small * -1); + ignore[right_child] = + (right_child_too_small * node) + (!right_child_too_small * -1); + + // If only left or right child is too small, consume it and relabel the other + // (to it can be its own cluster) + bool only_left_child_too_small = + left_child_too_small && !right_child_too_small; + bool only_right_child_too_small = + !left_child_too_small && right_child_too_small; + + relabel[right_child] = (only_left_child_too_small * relabel[node]) + + (!only_left_child_too_small * -1); + relabel[left_child] = (only_right_child_too_small * relabel[node]) + + (!only_right_child_too_small * -1); } } +struct Not_Empty { + +template +__host__ __device__ __forceinline__ value_t operator()(value_t a) { + return a != -1; +} +} + + +template +struct CondensedHierarchy { + + CondensedHierarchy(value_idx n_leaves_, cudaStream_t stream_): + n_leaves(n_leaves_), parents(0, stream_), children(0, stream_), lambdas(0, stream_), sizes(0, stream_) {} + + void condense(value_idx *full_parents, value_idx *full_children, + value_t *full_lambdas, value_idx *full_sizes) { + + n_edges = thrust::transform_reduce(thrust::cuda::par.on(stream), + full_parents, full_parents + (n_leaves * 2), + Not_Empty(), 0, thrust::plus()); + + thrust::copy_if(thrust::cuda::par.on(stream), full_parents, full_parents + (n_leaves * 2), parents.data(), Not_Empty()); + thrust::copy_if(thrust::cuda::par.on(stream), full_children, full_children + (n_leaves * 2), children.data(), Not_Empty()); + thrust::copy_if(thrust::cuda::par.on(stream), full_lambdas, full_lambdas + (n_leaves * 2), lambdas.data(), Not_Empty()); + thrust::copy_if(thrust::cuda::par.on(stream), full_sizes, full_sizes + (n_leaves * 2), sizes.data(), Not_Empty()); + } + + value_idx *get_parents() { + return parents.data(); + } + + value_idx *get_children() { + return children.data() + } + + value_t *get_lambdas() { + return lambdas.data(); + } + + value_idx *get_sizes() { + return sizes.data(); + } + + private: + rmm::device_uvector parents; + rmm::device_uvector children; + rmm::device_uvector lambdas; + rmm::device_uvector sizes; + + cudaStream_t stream; + value_idx n_edges; + value_idx n_leaves; + +}; + +/** + * Condenses a binary tree dendrogram in the Scipy format + * by merging labels that fall below a minimum cluster size. + * @tparam value_idx + * @tparam value_t + * @tparam tpb + * @param handle + * @param[in] children + * @param[in] delta + * @param[in] sizes + * @param[in] min_cluster_size + * @param[in] n_leaves + * @param[out] out_parent + * @param[out] out_child + * @param[out] out_lambda + * @param[out] out_size + */ template -void condense_hierarchy(raft::handle_t &handle, value_idx *children, - value_t *delta, value_idx *sizes, int min_pts, - int n_leaves) { - rmm::device_uvector frontier(n_leaves * 2, handle.get_stream()); - rmm::device_uvector ignore(n_leaves * 2, handle.get_stream()); +void condense_hierarchy(raft::handle_t &handle, const value_idx *children, + const value_t *delta, const value_idx *sizes, + int min_cluster_size, int n_leaves, + CondensedHierarchy &condensed_tree) { + + cudaStream_t stream = handle.get_stream(); + + rmm::device_uvector frontier(n_leaves * 2, stream); + rmm::device_uvector ignore(n_leaves * 2, stream); + + rmm::device_uvector out_parent(n_leaves * 2, stream); + rmm::device_uvector out_child(n_leaves * 2, stream); + rmm::device_uvector out_lambda(n_leaves * 2, stream); + rmm::device_uvector out_size(n_leaves * 2, stream); int root = 2 * n_leaves; int num_points = floor(root / 2.0) + 1; - rmm::device_uvector relabel(root + 1, handle.get_stream()); + thrust::fill(thrust::cuda::par.on(stream), out_parent.data(), out_parent.data()+(n_leaves*2), -1); + thrust::fill(thrust::cuda::par.on(stream), out_child.data(), out_parent.data()+(n_leaves*2), -1); + thrust::fill(thrust::cuda::par.on(stream), out_lambda.data(), out_parent.data()+(n_leaves*2), -1); + thrust::fill(thrust::cuda::par.on(stream), out_size.data(), out_parent.data()+(n_leaves*2), -1); - // TODO: Set this properly on device - relabel[root] = root; + rmm::device_uvector relabel(root + 1, handle.get_stream()); + raft::update_device(relabel.data()+root, root, 1, handle.get_stream()); // While frontier is not empty, perform single bfs through tree size_t grid = raft::ceildiv(n_leaves * 2, (size_t)tpb); @@ -160,17 +255,21 @@ void condense_hierarchy(raft::handle_t &handle, value_idx *children, n_elements_to_traverse = thrust::reduce(thrust::cuda::par.on(handle.get_stream()), frontier.data(), frontier.data() + (n_leaves * 2), 0); - - CUDA_CHECK(cudaStreamSynchronize(handle.get_stream())); } // TODO: Normalize labels so they are drawn from a monotonically increasing set. + + condensed_tree.condense(out_parent.data(), out_child.data(), out_lambda.data(), out_size.data()); } template -void compute_stabilities(value_idx *condensed_hierarchy, value_idx *lambdas, value_idx *sizes, +void compute_stabilities(value_idx *condensed_parent, + value_idx *condensed_child, + value_t *lambdas, + value_idx *sizes, int n_leaves) { + // TODO: Reverse topological sort (e.g. sort hierarchy, lambdas, and sizes by lambda) // TODO: Segmented reduction on min_lambda within each cluster @@ -214,7 +313,7 @@ void get_probabilities() { // TODO: Compute deaths array similarly to compute_stabilities - // TODO: Embarassingly parallel + // TODO: Embarassingly parallel } }; // end namespace Tree From b502e1df395592aaaf6777d7e5e510dbea190faf Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Thu, 22 Apr 2021 15:14:29 -0400 Subject: [PATCH 007/177] Resizing rmm arrays --- cpp/src/hdbscan/tree.cuh | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index 8bf3f48a7e..3e8853ae08 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -153,7 +153,8 @@ template struct CondensedHierarchy { CondensedHierarchy(value_idx n_leaves_, cudaStream_t stream_): - n_leaves(n_leaves_), parents(0, stream_), children(0, stream_), lambdas(0, stream_), sizes(0, stream_) {} + n_leaves(n_leaves_), parents(0, stream_), children(0, stream_), + lambdas(0, stream_), sizes(0, stream_) {} void condense(value_idx *full_parents, value_idx *full_children, value_t *full_lambdas, value_idx *full_sizes) { @@ -162,10 +163,19 @@ struct CondensedHierarchy { full_parents, full_parents + (n_leaves * 2), Not_Empty(), 0, thrust::plus()); - thrust::copy_if(thrust::cuda::par.on(stream), full_parents, full_parents + (n_leaves * 2), parents.data(), Not_Empty()); - thrust::copy_if(thrust::cuda::par.on(stream), full_children, full_children + (n_leaves * 2), children.data(), Not_Empty()); - thrust::copy_if(thrust::cuda::par.on(stream), full_lambdas, full_lambdas + (n_leaves * 2), lambdas.data(), Not_Empty()); - thrust::copy_if(thrust::cuda::par.on(stream), full_sizes, full_sizes + (n_leaves * 2), sizes.data(), Not_Empty()); + parents.resize(n_edges, stream); + children.resize(n_edges, stream); + lambdas.resize(n_edges, stream); + sizes.resize(n_edges, stream); + + thrust::copy_if(thrust::cuda::par.on(stream), full_parents, + full_parents + (n_leaves * 2), parents.data(), Not_Empty()); + thrust::copy_if(thrust::cuda::par.on(stream), + full_children, full_children + (n_leaves * 2), children.data(), Not_Empty()); + thrust::copy_if(thrust::cuda::par.on(stream), + full_lambdas, full_lambdas + (n_leaves * 2), lambdas.data(), Not_Empty()); + thrust::copy_if(thrust::cuda::par.on(stream), + full_sizes, full_sizes + (n_leaves * 2), sizes.data(), Not_Empty()); } value_idx *get_parents() { From 146ac005da745e0930b07a2ece733f8cd0ac21d1 Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Thu, 22 Apr 2021 15:15:05 -0400 Subject: [PATCH 008/177] Adding missing semicolon --- cpp/src/hdbscan/tree.cuh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index 3e8853ae08..a5a435f112 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -146,7 +146,7 @@ template __host__ __device__ __forceinline__ value_t operator()(value_t a) { return a != -1; } -} +}; template From 3cee98b47aa7338101988dc35add3a50e1f23ff8 Mon Sep 17 00:00:00 2001 From: divyegala Date: Thu, 22 Apr 2021 17:38:08 -0700 Subject: [PATCH 009/177] initial ideas to compute stabilities --- cpp/src/hdbscan/tree.cuh | 67 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 66 insertions(+), 1 deletion(-) diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index a5a435f112..d60bf31dda 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -272,19 +272,84 @@ void condense_hierarchy(raft::handle_t &handle, const value_idx *children, condensed_tree.condense(out_parent.data(), out_child.data(), out_lambda.data(), out_size.data()); } +template +struct transform_functor { + +public: + transform_op(value_t *stabilities_, value_t *births_) : + stabilities(stabilities_), + births(births_) { + + } + + __device__ value_t operator()(const &idx) { + return stabilities[idx] - births[idx]; + } + +private: + value_t *stabilities, *births; +}; + template -void compute_stabilities(value_idx *condensed_parent, +rmm::device_uvector compute_stabilities(value_idx *condensed_parent, value_idx *condensed_child, value_t *lambdas, value_idx *sizes, + int n_points, int n_leaves) { + auto thrust_policy = rmm::exec_policy(stream); // TODO: Reverse topological sort (e.g. sort hierarchy, lambdas, and sizes by lambda) + rmm::device_uvector sorted_child(condensed_child, n_points, stream); + rmm::device_uvector sorted_lambdas(lambdas, n_points, stream); + + auto children_lambda_zip = thrust::make_zip_iterator(thrust::make_tuple(sorted_child.begin(), sorted_lambdas.begin())); + thrust::sort_by_key(policy, condensed_parent, condensed_parent + n_points, children_lambda_zip); // TODO: Segmented reduction on min_lambda within each cluster + // TODO: Converting child array to CSR offset and using CUB Segmented Reduce + // Investigate use of a kernel like coo_spmv + auto n_clusters = // max label in parent - n_leaves, which make_monotonic will provide with + rmm::device_uvector birth(n_clusters, stream); + thrust::fill(thrust_policy, birth.begin(), birth.end(), 0); + + rmm::device_uvector sorted_child_offsets(n_points + 1, stream); + auto start_offset = 0; + sorted_child_offsets.set_element_async(0, start_offset, stream); + thrust::inclusive_scan(thrust_policy, sorted_child.begin(), sorted_child.end(), sorted_child_offsets.begin() + 1); + + void *d_temp_storage = NULL; + size_t temp_storage_bytes = 0; + cub::DeviceSegmentedReduce::Min(d_temp_storage, temp_storage_bytes, lambdas.begin(), birth.begin(), + n_clusters, sorted_child_offsets.begin(), sorted_child_offsets.begin() + 1); + CUDA_CHECK(cudaMalloc(&d_temp_storage, temp_storage_bytes)); + + cub::DeviceSegmentedReduce::Min(d_temp_storage, temp_storage_bytes, lambdas.begin(), birth.begin(), + n_clusters, sorted_child_offsets.begin(), sorted_child_offsets.begin() + 1); + CUDA_CHECK(cudaFree(d_temp_storage)); // TODO: Embarassingly parallel construction of output + // TODO: It can be done with same coo_spmv kernel + // Or naive kernel, atomically write to cluster stability + rmm::device_uvector stabilities(n_clusters, stream); + thrust::fill(thrust_policy, stabilities.begin(), stabilities.end(), 0); + + *d_temp_storage = NULL; + temp_storage_bytes = 0; + cub::DeviceSegmentedReduce::Sum(d_temp_storage, temp_storage_bytes, lambdas.begin(), stabilities.begin(), + n_clusters, sorted_child_offsets.begin(), sorted_child_offsets.begin() + 1); + CUDA_CHECK(cudaMalloc(&d_temp_storage, temp_storage_bytes)); + + cub::DeviceSegmentedReduce::Sum(d_temp_storage, temp_storage_bytes, lambdas.begin(), stabilities.begin(), + n_clusters, sorted_child_offsets.begin(), sorted_child_offsets.begin() + 1); + CUDA_CHECK(cudaFree(d_temp_storage)); + + // now transform + auto transform_op = transform_functor(stabilities.data(), birth.data()); + thrust::transform(policy, thrust::make_counting_iterator(0), thrust::make_counting_iterator(n_clusters), stabilities.begin(), transform_op); + + return stabilities; } template From 60253cec12f40cf8dda37b3b14c75d59a36c2eeb Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Thu, 22 Apr 2021 20:47:14 -0400 Subject: [PATCH 010/177] Checking in --- cpp/src/hdbscan/runner.h | 11 +++- cpp/src/hdbscan/tree.cuh | 133 ++++++++++++++++++++++++++++++--------- 2 files changed, 111 insertions(+), 33 deletions(-) diff --git a/cpp/src/hdbscan/runner.h b/cpp/src/hdbscan/runner.h index 5853d4c6a4..6f879b29a3 100644 --- a/cpp/src/hdbscan/runner.h +++ b/cpp/src/hdbscan/runner.h @@ -16,14 +16,15 @@ #pragma once -#include #include -#include +#include #include #include #include + +#include "tree.cuh" #include "reachability.cuh" namespace ML { @@ -48,7 +49,7 @@ struct MSTEpilogueReachability { template void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, raft::distance::DistanceType metric, int k, int min_pts, - float alpha) { + float alpha, int min_cluster_size) { auto d_alloc = handle.get_device_allocator(); auto stream = handle.get_stream(); @@ -94,6 +95,10 @@ void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, /** * Condense branches of tree according to min cluster size */ + Tree::CondensedHierarchy condensed_tree(m, stream); + condense_hierarchy(handle, out_src.data(), out_dst.data(), + out_delta.data(), out_size.data(), + min_cluster_size, m, condensed_tree); /** * Extract labels from stability diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index a5a435f112..dc4cb16566 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -20,6 +20,9 @@ #include +#include +#include + #include #include @@ -53,7 +56,7 @@ __device__ value_t get_lambda(value_idx node, value_idx num_points, template __global__ void condense_hierarchy_kernel( bool *frontier, value_idx *ignore, value_idx *relabel, - const value_idx *hierarchy, const value_t *deltas, + const value_idx *src, const value_idx *dst, const value_t *deltas, const value_idx *sizes, int n_leaves, int num_points, int min_cluster_size, value_idx *out_parent, value_idx *out_child, @@ -68,8 +71,8 @@ __global__ void condense_hierarchy_kernel( frontier[node] = false; // TODO: Check bounds - value_idx left_child = hierarchy[(node - num_points) * 2]; - value_idx right_child = hierarchy[((node - num_points) * 2) + 1]; + value_idx left_child = src[(node - num_points) * 2]; + value_idx right_child = dst[((node - num_points) * 2)]; frontier[left_child] = true; frontier[right_child] = true; @@ -92,8 +95,8 @@ __global__ void condense_hierarchy_kernel( // If node is not ignored and is not a leaf, condense its children // if necessary else if (!should_ignore and node >= num_points) { - value_idx left_child = hierarchy[(node - num_points) * 2]; - value_idx right_child = hierarchy[((node - num_points) * 2) + 1]; + value_idx left_child = src[(node - num_points) * 2]; + value_idx right_child = dst[((node - num_points) * 2)]; value_t lambda_value = get_lambda(node, num_points, deltas); @@ -107,13 +110,13 @@ __global__ void condense_hierarchy_kernel( if (left_count >= min_cluster_size && right_count >= min_cluster_size) { relabel[left_child] = node; out_parent[node] = relabel[node]; - out_child[node] = relabel[left_child]; + out_child[node] = node; out_lambda[node] = lambda_value; out_count[node] = left_count; relabel[right_child] = node; out_parent[node] = relabel[node]; - out_child[node] = relabel[right_child]; + out_child[node] = node; out_lambda[node] = lambda_value; out_count[node] = left_count; } @@ -194,6 +197,10 @@ struct CondensedHierarchy { return sizes.data(); } + value_idx get_n_edges() { + return n_edges; + } + private: rmm::device_uvector parents; rmm::device_uvector children; @@ -224,7 +231,8 @@ struct CondensedHierarchy { * @param[out] out_size */ template -void condense_hierarchy(raft::handle_t &handle, const value_idx *children, +void condense_hierarchy(const raft::handle_t &handle, const value_idx *src, + const value_idx *dst, const value_t *delta, const value_idx *sizes, int min_cluster_size, int n_leaves, CondensedHierarchy &condensed_tree) { @@ -242,10 +250,14 @@ void condense_hierarchy(raft::handle_t &handle, const value_idx *children, int root = 2 * n_leaves; int num_points = floor(root / 2.0) + 1; - thrust::fill(thrust::cuda::par.on(stream), out_parent.data(), out_parent.data()+(n_leaves*2), -1); - thrust::fill(thrust::cuda::par.on(stream), out_child.data(), out_parent.data()+(n_leaves*2), -1); - thrust::fill(thrust::cuda::par.on(stream), out_lambda.data(), out_parent.data()+(n_leaves*2), -1); - thrust::fill(thrust::cuda::par.on(stream), out_size.data(), out_parent.data()+(n_leaves*2), -1); + thrust::fill(thrust::cuda::par.on(stream), out_parent.data(), + out_parent.data()+(n_leaves*2), -1); + thrust::fill(thrust::cuda::par.on(stream), out_child.data(), + out_child.data()+(n_leaves*2), -1); + thrust::fill(thrust::cuda::par.on(stream), out_lambda.data(), + out_lambda.data()+(n_leaves*2), -1); + thrust::fill(thrust::cuda::par.on(stream), out_size.data(), + out_size.data()+(n_leaves*2), -1); rmm::device_uvector relabel(root + 1, handle.get_stream()); raft::update_device(relabel.data()+root, root, 1, handle.get_stream()); @@ -260,7 +272,7 @@ void condense_hierarchy(raft::handle_t &handle, const value_idx *children, while (n_elements_to_traverse > 0) { condense_hierarchy_kernel<<>>( frontier.data(), ignore.data(), next_label.data(), relabel.data(), - children, delta, sizes, n_leaves, num_points, min_cluster_size); + src, dst, delta, sizes, n_leaves, num_points, min_cluster_size); n_elements_to_traverse = thrust::reduce(thrust::cuda::par.on(handle.get_stream()), frontier.data(), @@ -269,37 +281,98 @@ void condense_hierarchy(raft::handle_t &handle, const value_idx *children, // TODO: Normalize labels so they are drawn from a monotonically increasing set. - condensed_tree.condense(out_parent.data(), out_child.data(), out_lambda.data(), out_size.data()); + condensed_tree.condense(out_parent.data(), out_child.data(), + out_lambda.data(), out_size.data()); } template -void compute_stabilities(value_idx *condensed_parent, - value_idx *condensed_child, - value_t *lambdas, - value_idx *sizes, - int n_leaves) { +void compute_stabilities(const raft::handle_t &handle, + const CondensedHierarchy &condensed_tree) { - - // TODO: Reverse topological sort (e.g. sort hierarchy, lambdas, and sizes by lambda) + // TODO: sort hierarchy, lambdas, and sizes by lambda // TODO: Segmented reduction on min_lambda within each cluster // TODO: Embarassingly parallel construction of output + + } +struct Greater_Than_One { + + template + __host__ __device__ __forceinline__ value_t operator()(value_t a) { + return a > 1; + } +}; + template -void excess_of_mass() { +void excess_of_mass(const raft::handle_t &handle, + const CondensedHierarchy &condensed_tree, + value_t *stability, bool *is_cluster, value_idx n_clusters) { + + /** + * - If the sum of the stabilities of the child clusters is greater than the + * stability of the cluster, then we set the cluster stability to be the + * sum of the child stabilities. + * - If, on the other hand, the cluster’s stability is greater than the sum + * of its children then we declare the cluster to be a selected cluster + * and unselect all its descendants. + * - Once we reach the root node we call the current set of selected clusters + * our flat clustering and return that. + */ + + cudaStream_t stream = handle.get_stream(); + + /** + * 1. Build CSR of cluster tree from condensed tree by filtering condensed tree for + * only those entries w/ lambda > 1 and constructing a CSR from the result + */ + + + value_idx cluster_tree_edges = thrust::transform_reduce(thrust::cuda::par.on(stream), + condensed_tree.get_lambdas(), + condensed_tree.get_lambdas() + condensed_tree.get_n_edges(), + Greater_Than_One(), 0, thrust::plus()); + + rmm::device_uvector parents(cluster_tree_edges, stream); + rmm::device_uvector children(cluster_tree_edges, stream); + rmm::device_uvector sizes(cluster_tree_edges, stream); + rmm::device_uvector indptr(n_clusters, stream); + + thrust::copy_if(thrust::cuda::par.on(stream), condensed_tree.get_parents(), + condensed_tree.get_parents() + (condensed_tree.get_n_edges()), condensed_tree.get_lambdas(), + parents.data(), Greater_Than_One()); + + thrust::copy_if(thrust::cuda::par.on(stream), condensed_tree.get_children(), + condensed_tree.get_children() + (condensed_tree.get_n_edges()), condensed_tree.get_lambdas(), + children.data(), Greater_Than_One()); + + thrust::copy_if(thrust::cuda::par.on(stream), condensed_tree.get_sizes(), + condensed_tree.get_sizes() + (condensed_tree.get_n_edges()), condensed_tree.get_lambdas(), + sizes.data(), Greater_Than_One()); + + raft::sparse::op::coo_sort(0, 0, cluster_tree_edges, parents.data(), children.data(), sizes.data(), + handle.get_device_allocator(), handle.get_stream()); + + raft::sparse::convert::sorted_coo_to_csr(parents.data(), cluster_tree_edges, indptr.data(), n_clusters, + handle.get_device_allocator(), handle.get_stream()); - // TODO: Build CSR of cluster tree with stabilities of each child as the weights + /** + * 2. Iterate through each level from leaves back to root. Use the cluster + * tree CSR and warp-level reduction to sum stabilities and test whether + * or not current cluster should continue to be its own + */ + /** + * Copy indptr to host + * For each node in sorted stability keys, + * - transformed reduce + */ - // TODO: Segmented reduction over CSR of cluster tree + /** + * 3. Perform BFS through is_cluster, propagating cluster "deselection" to leaves + */ - // TODO: Perform bfs, starting at root- - // TODO: Maintain frontier and is_cluster array. - // TODO: In each iteration, children are added to tree - // TODO: If node has is_cluster[node] = false, set children to false - // TODO: else subtree stability > stability[node] or cluster_sizes[node] > max_cluster_size - // TODO: set is_cluster[node] = false and stability[node] = subtree_stability } template From e4668d3de3a86ab9b89a7f830d6dbcaf60421bca Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Fri, 23 Apr 2021 09:48:26 -0400 Subject: [PATCH 011/177] Updating --- cpp/src/hdbscan/tree.cuh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index dc4cb16566..2faa62763c 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -366,7 +366,7 @@ void excess_of_mass(const raft::handle_t &handle, /** * Copy indptr to host * For each node in sorted stability keys, - * - transformed reduce + * - transformed reducet */ /** From fd77b481a2aeb120679bd0d158ae2e907ccdb968 Mon Sep 17 00:00:00 2001 From: divyegala Date: Fri, 23 Apr 2021 09:41:21 -0700 Subject: [PATCH 012/177] allocating stabilities through caller --- cpp/src/hdbscan/runner.h | 6 ++++-- cpp/src/hdbscan/tree.cuh | 11 +++++++---- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/cpp/src/hdbscan/runner.h b/cpp/src/hdbscan/runner.h index 4e3458ade4..85f2890bfd 100644 --- a/cpp/src/hdbscan/runner.h +++ b/cpp/src/hdbscan/runner.h @@ -95,12 +95,14 @@ void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, /** * Condense branches of tree according to min cluster size */ + int n_condensed_clusters; Tree::CondensedHierarchy condensed_tree(m, stream); condense_hierarchy(handle, out_src.data(), out_dst.data(), out_delta.data(), out_size.data(), - min_cluster_size, m, condensed_tree); + min_cluster_size, m, condensed_tree, n_condensed_clusters); - rmm::device_uvector stabilities = compute_stabilities(handle, condensed_tree); + rmm::device_uvector stabilities(n_condensed_clusters, handle.get_stream()); + compute_stabilities(handle, condensed_tree, stabilities, n_condensed_clusters); /** * Extract labels from stability */ diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index ebe419eaac..25bf460af0 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -126,7 +126,8 @@ void condense_hierarchy(const raft::handle_t &handle, const value_idx *src, const value_idx *dst, const value_t *delta, const value_idx *sizes, int min_cluster_size, int n_leaves, - CondensedHierarchy &condensed_tree) { + CondensedHierarchy &condensed_tree + int &n_condensed_clusters) { cudaStream_t stream = handle.get_stream(); @@ -177,7 +178,7 @@ void condense_hierarchy(const raft::handle_t &handle, const value_idx *src, auto parents = condensed_tree.get_parents(); auto parents_len = condensed_tree.get_n_edges(); - MLCommon::Label::make_monotonic(parents, parents, parents_len, stream, handle.get_device_allocator()); + n_condensed_clusters = MLCommon::Label::make_monotonic(parents, parents, parents_len, stream, handle.get_device_allocator()); } template @@ -199,8 +200,10 @@ private: }; template -rmm::device_uvector compute_stabilities(const raft::handle_t &handle, - const CondensedHierarchy &condensed_tree) { +void compute_stabilities(const raft::handle_t &handle, + const CondensedHierarchy &condensed_tree, + rmm::device_uvector &stabilities, + int n_condensed_clusters) { auto parents = condensed_tree.get_parents(); auto children = condensed_tree.get_children(); From 7edb4a3205861704b9aa41c2f8f0561888525139 Mon Sep 17 00:00:00 2001 From: divyegala Date: Fri, 23 Apr 2021 10:47:44 -0700 Subject: [PATCH 013/177] cleaning stabilities a bit --- cpp/src/hdbscan/runner.h | 7 ++- cpp/src/hdbscan/tree.cuh | 82 ++++++++++++++--------------- cpp/src_prims/label/classlabels.cuh | 5 +- 3 files changed, 46 insertions(+), 48 deletions(-) diff --git a/cpp/src/hdbscan/runner.h b/cpp/src/hdbscan/runner.h index 85f2890bfd..bef9060b72 100644 --- a/cpp/src/hdbscan/runner.h +++ b/cpp/src/hdbscan/runner.h @@ -95,14 +95,13 @@ void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, /** * Condense branches of tree according to min cluster size */ - int n_condensed_clusters; Tree::CondensedHierarchy condensed_tree(m, stream); condense_hierarchy(handle, out_src.data(), out_dst.data(), out_delta.data(), out_size.data(), - min_cluster_size, m, condensed_tree, n_condensed_clusters); + min_cluster_size, m, condensed_tree); - rmm::device_uvector stabilities(n_condensed_clusters, handle.get_stream()); - compute_stabilities(handle, condensed_tree, stabilities, n_condensed_clusters); + rmm::device_uvector stabilities(condensed_tree.get_n_clusters(), handle.get_stream()); + compute_stabilities(handle, condensed_tree, stabilities); /** * Extract labels from stability */ diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index 25bf460af0..6ce5347f78 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -46,13 +46,15 @@ __host__ __device__ __forceinline__ value_t operator()(value_t a) { template struct CondensedHierarchy { - CondensedHierarchy(value_idx n_leaves_, cudaStream_t stream_): - n_leaves(n_leaves_), parents(0, stream_), children(0, stream_), - lambdas(0, stream_), sizes(0, stream_) {} + CondensedHierarchy(const raft::handle_t &handle_, value_idx n_leaves_): + handle(handle_), n_leaves(n_leaves_), parents(0, handle.get_stream()), children(0, handle.get_stream()), + lambdas(0, handle.get_stream()), sizes(0, handle.get_stream()) {} void condense(value_idx *full_parents, value_idx *full_children, value_t *full_lambdas, value_idx *full_sizes) { + auto stream = handle.get_stream(); + n_edges = thrust::transform_reduce(thrust::cuda::par.on(stream), full_parents, full_parents + (n_leaves * 2), Not_Empty(), 0, thrust::plus()); @@ -70,6 +72,8 @@ struct CondensedHierarchy { full_lambdas, full_lambdas + (n_leaves * 2), lambdas.data(), Not_Empty()); thrust::copy_if(thrust::cuda::par.on(stream), full_sizes, full_sizes + (n_leaves * 2), sizes.data(), Not_Empty()); + + n_clusters = MLCommon::Label::make_monotonic(handle, parents.data(), parents.begin(), parents.end()); } value_idx *get_parents() { @@ -92,15 +96,21 @@ struct CondensedHierarchy { return n_edges; } + int get_n_clusters() { + return n_clusters; + } + private: + const raft::handle_t &handle; + rmm::device_uvector parents; rmm::device_uvector children; rmm::device_uvector lambdas; rmm::device_uvector sizes; - cudaStream_t stream; value_idx n_edges; value_idx n_leaves; + int n_clusters; }; @@ -126,8 +136,7 @@ void condense_hierarchy(const raft::handle_t &handle, const value_idx *src, const value_idx *dst, const value_t *delta, const value_idx *sizes, int min_cluster_size, int n_leaves, - CondensedHierarchy &condensed_tree - int &n_condensed_clusters) { + CondensedHierarchy &condensed_tree) { cudaStream_t stream = handle.get_stream(); @@ -176,9 +185,6 @@ void condense_hierarchy(const raft::handle_t &handle, const value_idx *src, condensed_tree.condense(out_parent.data(), out_child.data(), out_lambda.data(), out_size.data()); - auto parents = condensed_tree.get_parents(); - auto parents_len = condensed_tree.get_n_edges(); - n_condensed_clusters = MLCommon::Label::make_monotonic(parents, parents, parents_len, stream, handle.get_device_allocator()); } template @@ -199,68 +205,62 @@ private: value_t *stabilities, *births; }; -template +template +void segmented_reduce(const value_t *in, value_t *out, const value_idx *offsets, cudaStream_t stream) { + void *d_temp_storage = NULL; + size_t temp_storage_bytes = 0; + cub_reduce_func(d_temp_storage, temp_storage_bytes, in, out, + n_clusters, offsets, offsets + 1, stream); + CUDA_CHECK(cudaMalloc(&d_temp_storage, temp_storage_bytes)); + + cub_reduce_func(d_temp_storage, temp_storage_bytes, in, out, + n_clusters, offsets, offsets + 1, stream); + CUDA_CHECK(cudaFree(d_temp_storage)); +} + +template void compute_stabilities(const raft::handle_t &handle, const CondensedHierarchy &condensed_tree, - rmm::device_uvector &stabilities, - int n_condensed_clusters) { + rmm::device_uvector &stabilities) { auto parents = condensed_tree.get_parents(); auto children = condensed_tree.get_children(); auto lambdas = condensed_tree.get_lambdas(); auto n_edges = condensed_tree.get_n_edges(); + auto n_clusters = condensed_tree.get_n_clusters(); auto stream = handle.get_stream(); auto thrust_policy = rmm::exec_policy(stream); // TODO: Reverse topological sort (e.g. sort hierarchy, lambdas, and sizes by lambda) - rmm::device_uvector sorted_child(condensed_child, n_points, stream); - rmm::device_uvector sorted_lambdas(lambdas, n_points, stream); + rmm::device_uvector sorted_child(condensed_child, n_edges, stream); + rmm::device_uvector sorted_lambdas(lambdas, n_edges, stream); auto children_lambda_zip = thrust::make_zip_iterator(thrust::make_tuple(sorted_child.begin(), sorted_lambdas.begin())); - thrust::sort_by_key(policy, condensed_parent, condensed_parent + n_points, children_lambda_zip); + thrust::sort_by_key(policy, parents, parents + n_edges, children_lambda_zip); // TODO: sort hierarchy, lambdas, and sizes by lambda // TODO: Segmented reduction on min_lambda within each cluster // TODO: Converting child array to CSR offset and using CUB Segmented Reduce // Investigate use of a kernel like coo_spmv - auto n_clusters = // max label in parent - n_leaves, which make_monotonic will provide with - rmm::device_uvector birth(n_clusters, stream); - thrust::fill(thrust_policy, birth.begin(), birth.end(), 0); + rmm::device_uvector births(n_clusters, stream); + thrust::fill(thrust_policy, births.begin(), births.end(), 0); - rmm::device_uvector sorted_child_offsets(n_points + 1, stream); - auto start_offset = 0; - sorted_child_offsets.set_element_async(0, start_offset, stream); - thrust::inclusive_scan(thrust_policy, sorted_child.begin(), sorted_child.end(), sorted_child_offsets.begin() + 1); + rmm::device_uvector sorted_child_offsets(n_edges + 1, stream); - void *d_temp_storage = NULL; - size_t temp_storage_bytes = 0; - cub::DeviceSegmentedReduce::Min(d_temp_storage, temp_storage_bytes, lambdas.begin(), birth.begin(), - n_clusters, sorted_child_offsets.begin(), sorted_child_offsets.begin() + 1); - CUDA_CHECK(cudaMalloc(&d_temp_storage, temp_storage_bytes)); + raft::sparse::convert::sorted_coo_to_csr(sorted_child.data(), n_edges, sorted_child_offsets.data(), n_clusters, handle.get_stream(), handle.get_device_allocator()); - cub::DeviceSegmentedReduce::Min(d_temp_storage, temp_storage_bytes, lambdas.begin(), birth.begin(), - n_clusters, sorted_child_offsets.begin(), sorted_child_offsets.begin() + 1); - CUDA_CHECK(cudaFree(d_temp_storage)); + segmented_reduce(lambdas, births.data(), sorted_child_offsets.data(), stream); // TODO: Embarassingly parallel construction of output // TODO: It can be done with same coo_spmv kernel // Or naive kernel, atomically write to cluster stability - rmm::device_uvector stabilities(n_clusters, stream); thrust::fill(thrust_policy, stabilities.begin(), stabilities.end(), 0); - *d_temp_storage = NULL; - temp_storage_bytes = 0; - cub::DeviceSegmentedReduce::Sum(d_temp_storage, temp_storage_bytes, lambdas.begin(), stabilities.begin(), - n_clusters, sorted_child_offsets.begin(), sorted_child_offsets.begin() + 1); - CUDA_CHECK(cudaMalloc(&d_temp_storage, temp_storage_bytes)); - - cub::DeviceSegmentedReduce::Sum(d_temp_storage, temp_storage_bytes, lambdas.begin(), stabilities.begin(), - n_clusters, sorted_child_offsets.begin(), sorted_child_offsets.begin() + 1); - CUDA_CHECK(cudaFree(d_temp_storage)); + segmented_reduce(lambdas, stabilities.data(), sorted_child_offsets.data(), stream); - // now transform + // now transform, and calculate summation lambda(point) - lambda(birth) auto transform_op = transform_functor(stabilities.data(), birth.data()); thrust::transform(policy, thrust::make_counting_iterator(0), thrust::make_counting_iterator(n_clusters), stabilities.begin(), transform_op); diff --git a/cpp/src_prims/label/classlabels.cuh b/cpp/src_prims/label/classlabels.cuh index 74ed2bc871..a8f511f5ce 100644 --- a/cpp/src_prims/label/classlabels.cuh +++ b/cpp/src_prims/label/classlabels.cuh @@ -191,10 +191,9 @@ void make_monotonic(Type *out, Type *in, size_t N, cudaStream_t stream, } template -int make_monotonic(Type *out, Type *in, size_t N, cudaStream_t stream, - std::shared_ptr allocator) { +int make_monotonic(const raft::handle_t &handle, Type *out, Type *in, size_t N) { return make_monotonic( - out, in, N, stream, [] __device__(Type val) { return false; }, allocator); + out, in, N, handle.get_stream(), [] __device__(Type val) { return false; }, handle.get_device_allocator()); } }; // namespace Label }; // end namespace MLCommon From 0efd58b748ec252c0020834dc46fd0e08418a2e9 Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Fri, 23 Apr 2021 13:58:57 -0400 Subject: [PATCH 014/177] Initial stub for eom is done. --- cpp/src/hdbscan/detail/tree_kernels.cuh | 23 ++++++ cpp/src/hdbscan/tree.cuh | 99 ++++++++++++++++++++----- 2 files changed, 103 insertions(+), 19 deletions(-) diff --git a/cpp/src/hdbscan/detail/tree_kernels.cuh b/cpp/src/hdbscan/detail/tree_kernels.cuh index d98d347ed1..c216e405ac 100644 --- a/cpp/src/hdbscan/detail/tree_kernels.cuh +++ b/cpp/src/hdbscan/detail/tree_kernels.cuh @@ -133,6 +133,29 @@ __global__ void condense_hierarchy_kernel( (!only_right_child_too_small * -1); } } + +template +__global__ void propagate_cluster_negation(const value_idx *indptr, + const value_idx *children, + bool *frontier, + bool *is_cluster, + int n_clusters) { + + int cluster = blockDim.x * blockIdx.x + threadIdx.x; + + if(cluster < n_clusters && frontier[cluster]) { + frontier[cluster] = false; + + value_idx children_start = indptr[cluster]; + value_idx children_stop = indptr[cluster]; + for(int i = 0; i < children_stop - children_start; i++) { + + value_idx child = children[i]; + frontier[child] = true; + is_cluster[child] = false; + } + } +} }; // end namespace detail diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index 6ce5347f78..bac2e5335c 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -81,7 +81,7 @@ struct CondensedHierarchy { } value_idx *get_children() { - return children.data() + return children.data(); } value_t *get_lambdas() { @@ -171,6 +171,9 @@ void condense_hierarchy(const raft::handle_t &handle, const value_idx *src, frontier.data() + (n_leaves * 2), 0); while (n_elements_to_traverse > 0) { + + // TODO: Investigate whether it would be worth performing a gather/argmatch in order + // to schedule only the number of threads needed. (it might not be worth it) detail::condense_hierarchy_kernel<<>>( frontier.data(), ignore.data(), next_label.data(), relabel.data(), src, dst, delta, sizes, n_leaves, num_points, min_cluster_size); @@ -194,7 +197,6 @@ public: transform_op(value_t *stabilities_, value_t *births_) : stabilities(stabilities_), births(births_) { - } __device__ value_t operator()(const &idx) { @@ -232,15 +234,13 @@ void compute_stabilities(const raft::handle_t &handle, auto stream = handle.get_stream(); auto thrust_policy = rmm::exec_policy(stream); - // TODO: Reverse topological sort (e.g. sort hierarchy, lambdas, and sizes by lambda) - rmm::device_uvector sorted_child(condensed_child, n_edges, stream); - rmm::device_uvector sorted_lambdas(lambdas, n_edges, stream); + // TODO: sort hierarchy, lambdas, and sizes by lambda + rmm::device_uvector sorted_child(condensed_child, n_points, stream); + rmm::device_uvector sorted_lambdas(lambdas, n_points, stream); auto children_lambda_zip = thrust::make_zip_iterator(thrust::make_tuple(sorted_child.begin(), sorted_lambdas.begin())); thrust::sort_by_key(policy, parents, parents + n_edges, children_lambda_zip); - // TODO: sort hierarchy, lambdas, and sizes by lambda - // TODO: Segmented reduction on min_lambda within each cluster // TODO: Converting child array to CSR offset and using CUB Segmented Reduce // Investigate use of a kernel like coo_spmv @@ -249,7 +249,10 @@ void compute_stabilities(const raft::handle_t &handle, rmm::device_uvector sorted_child_offsets(n_edges + 1, stream); - raft::sparse::convert::sorted_coo_to_csr(sorted_child.data(), n_edges, sorted_child_offsets.data(), n_clusters, handle.get_stream(), handle.get_device_allocator()); + raft::sparse::convert::sorted_coo_to_csr(sorted_child.data(), n_edges, + sorted_child_offsets.data(), n_clusters, + handle.get_stream(), + handle.get_device_allocator()); segmented_reduce(lambdas, births.data(), sorted_child_offsets.data(), stream); @@ -262,7 +265,8 @@ void compute_stabilities(const raft::handle_t &handle, // now transform, and calculate summation lambda(point) - lambda(birth) auto transform_op = transform_functor(stabilities.data(), birth.data()); - thrust::transform(policy, thrust::make_counting_iterator(0), thrust::make_counting_iterator(n_clusters), stabilities.begin(), transform_op); + thrust::transform(policy, thrust::make_counting_iterator(0), + thrust::make_counting_iterator(n_clusters), stabilities.begin(), transform_op); return stabilities; } @@ -275,10 +279,21 @@ struct Greater_Than_One { } }; -template + +struct Negate { + + template + __host__ +}; + + +template void excess_of_mass(const raft::handle_t &handle, const CondensedHierarchy &condensed_tree, - value_t *stability, bool *is_cluster, value_idx n_clusters) { + value_t *stability, + bool *is_cluster, + value_idx n_clusters, + value_idx max_cluster_size) { /** * - If the sum of the stabilities of the child clusters is greater than the @@ -298,6 +313,7 @@ void excess_of_mass(const raft::handle_t &handle, * only those entries w/ lambda > 1 and constructing a CSR from the result */ + std::vector cluster_sizes; value_idx cluster_tree_edges = thrust::transform_reduce(thrust::cuda::par.on(stream), condensed_tree.get_lambdas(), @@ -324,7 +340,7 @@ void excess_of_mass(const raft::handle_t &handle, raft::sparse::op::coo_sort(0, 0, cluster_tree_edges, parents.data(), children.data(), sizes.data(), handle.get_device_allocator(), handle.get_stream()); - raft::sparse::convert::sorted_coo_to_csr(parents.data(), cluster_tree_edges, indptr.data(), n_clusters, + raft::sparse::convert::sorted_coo_to_csr(parents.data(), cluster_tree_edges, indptr.data(), n_clusters+1, handle.get_device_allocator(), handle.get_stream()); /** @@ -332,16 +348,61 @@ void excess_of_mass(const raft::handle_t &handle, * tree CSR and warp-level reduction to sum stabilities and test whether * or not current cluster should continue to be its own */ - /** - * Copy indptr to host - * For each node in sorted stability keys, - * - transformed reducet - */ + std::vector indptr_h(indptr.size()); + raft::update_host(indptr_h.data(), indptr.data(), indptr.size(), stream); + CUDA_CHECK(cudaStreamSynchronize(stream)); + + std::vector is_cluster_h(n_clusters, true); + + for(value_idx node = 0; node < n_clusters; node++) { + + value_t node_stability; + raft::update_host(&node_stability, stability+node, 1, stream); + + value_t subtree_stability = thrust::transform_reduce(thrust::cuda::par.on(stream), + children.data()+indptr_h[node], children.data()+indptr_h[node]+1, + [=]__device__ (value_idx a) { return stability[a]; }, 0, thrust::plus()); + + if(subtree_stability > stability[node] || cluster_sizes[node] > max_cluster_size) { + // Deselect / merge cluster with children + raft::update_device(stability + node, subtree_stability, 1, stream); + } else { + // Deselect children + is_cluster_h[node] = false; + } + } /** - * 3. Perform BFS through is_cluster, propagating cluster "deselection" to leaves + * 3. Perform BFS through is_cluster, propagating cluster "deselection" through subtrees */ + raft::update_device(is_cluster, is_cluster_h.data(), n_clusters, stream); + rmm::device_uvector frontier(n_clusters, stream); + + thrust::transform(thrust::cuda::par.on(stream), + is_cluster, is_cluster+n_clusters, + frontier.data(), + [=]__device__ (value_t a) {return !a;}); + + + value_idx n_elements_to_traverse = + thrust::reduce(thrust::cuda::par.on(handle.get_stream()), frontier.data(), + frontier.data() + frontier.size(), 0); + + // TODO: Investigate whether it's worth gathering the sparse frontier into + // a dense form for purposes of uniform workload/thread scheduling + + // While frontier is not empty, perform single bfs through tree + size_t grid = raft::ceildiv(n_leaves * 2, (size_t)tpb); + + while (n_elements_to_traverse > 0) { + detail::propagate_cluster_negation<<>>( + indptr.data(), children.data(), frontier.data(), is_cluster, n_clusters); + + n_elements_to_traverse = + thrust::reduce(thrust::cuda::par.on(handle.get_stream()), frontier.data(), + frontier.data() + frontier.size(), 0); + } } template @@ -370,4 +431,4 @@ void get_probabilities() { }; // end namespace Tree }; // end namespace HDBSCAN -}; // end namespace ML \ No newline at end of file +}; // end namespace ML From bbb35dbbaf7502eff95064d68b02586320468d43 Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Fri, 23 Apr 2021 14:25:39 -0400 Subject: [PATCH 015/177] Adding data to epilogue --- cpp/src/hdbscan/runner.h | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/cpp/src/hdbscan/runner.h b/cpp/src/hdbscan/runner.h index bef9060b72..1f1c8f6aba 100644 --- a/cpp/src/hdbscan/runner.h +++ b/cpp/src/hdbscan/runner.h @@ -36,9 +36,13 @@ struct MSTEpilogueReachability { : core_distances(core_distances_), m(m_) {} void operator()(raft::handle_t &handle, value_idx *coo_rows, - value_idx *coo_cols, value_t *coo_data) { - // TODO: Schedule kernel that uses the core distances - // to perform the max(core_dist(src), core_dist(dst), d(src, dst)) operation + value_idx *coo_cols, value_t *coo_data, value_idx nnz) { + + auto first = thrust::make_zip_iterator(thrust::make_tuple(coo_rows, coo_cols, coo_data)); + thrust::transform(thrust::cuda::par.on(handle.get_stream()), first, first+nnz, + coo_data, [=] __device__ (thrust::tuple t) { + return max(core_distances[thrust::get<0>(t)], core_distances[thrust::get<1>(t)], thrust::get<2>(t)); + }); } private: @@ -100,11 +104,14 @@ void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, out_delta.data(), out_size.data(), min_cluster_size, m, condensed_tree); - rmm::device_uvector stabilities(condensed_tree.get_n_clusters(), handle.get_stream()); + rmm::device_uvector stabilities(condensed_tree.get_n_clusters(), + handle.get_stream()); compute_stabilities(handle, condensed_tree, stabilities); + /** * Extract labels from stability */ + } }; // end namespace HDBSCAN From 83d8ac3f627988412dcc1ea8369d01a6c474af1f Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Fri, 23 Apr 2021 15:09:00 -0400 Subject: [PATCH 016/177] COuple small updates --- cpp/src/hdbscan/runner.h | 9 +++++---- cpp/src/hdbscan/tree.cuh | 4 +--- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/cpp/src/hdbscan/runner.h b/cpp/src/hdbscan/runner.h index 1f1c8f6aba..4d734aae42 100644 --- a/cpp/src/hdbscan/runner.h +++ b/cpp/src/hdbscan/runner.h @@ -35,20 +35,22 @@ struct MSTEpilogueReachability { MSTEpilogueReachability(value_idx m_, value_t *core_distances_) : core_distances(core_distances_), m(m_) {} - void operator()(raft::handle_t &handle, value_idx *coo_rows, + void operator()(const raft::handle_t &handle, value_idx *coo_rows, value_idx *coo_cols, value_t *coo_data, value_idx nnz) { auto first = thrust::make_zip_iterator(thrust::make_tuple(coo_rows, coo_cols, coo_data)); thrust::transform(thrust::cuda::par.on(handle.get_stream()), first, first+nnz, coo_data, [=] __device__ (thrust::tuple t) { - return max(core_distances[thrust::get<0>(t)], core_distances[thrust::get<1>(t)], thrust::get<2>(t)); + return max(core_distances[thrust::get<0>(t)], + core_distances[thrust::get<1>(t)], + thrust::get<2>(t)); }); } private: value_t *core_distances; value_idx m; -} +}; template void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, @@ -60,7 +62,6 @@ void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, /** * Mutual reachability graph */ - rmm::device_uvector mutual_reachability_graph_inds(k * m, stream); rmm::device_uvector mutual_reachability_graph_dists(k * m, stream); rmm::device_uvector core_dists(k * m, stream); diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index bac2e5335c..3e93a9839c 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -416,9 +416,7 @@ void get_stability_scores() { template void do_labelling() { - // TODO: Similar to SLHC dendrogram construction, this one is probably best done - // on host, at least for the first iteration -} + // TODO: union find is constructed on host } template From f563acacb32f49f3a1b97511464ae905d450b32c Mon Sep 17 00:00:00 2001 From: "Corey J. Nolet" Date: Fri, 23 Apr 2021 15:45:08 -0400 Subject: [PATCH 017/177] Adding get_stability scores --- cpp/src/hdbscan/tree.cuh | 51 ++++++++++++++++++++++++++++++++++------ 1 file changed, 44 insertions(+), 7 deletions(-) diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index 3e93a9839c..49e830431c 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -238,7 +238,8 @@ void compute_stabilities(const raft::handle_t &handle, rmm::device_uvector sorted_child(condensed_child, n_points, stream); rmm::device_uvector sorted_lambdas(lambdas, n_points, stream); - auto children_lambda_zip = thrust::make_zip_iterator(thrust::make_tuple(sorted_child.begin(), sorted_lambdas.begin())); + auto children_lambda_zip = thrust::make_zip_iterator(thrust::make_tuple(sorted_child.begin(), + sorted_lambdas.begin())); thrust::sort_by_key(policy, parents, parents + n_edges, children_lambda_zip); // TODO: Segmented reduction on min_lambda within each cluster @@ -323,7 +324,7 @@ void excess_of_mass(const raft::handle_t &handle, rmm::device_uvector parents(cluster_tree_edges, stream); rmm::device_uvector children(cluster_tree_edges, stream); rmm::device_uvector sizes(cluster_tree_edges, stream); - rmm::device_uvector indptr(n_clusters, stream); + rmm::device_uvector indptr(n_clusters+1, stream); thrust::copy_if(thrust::cuda::par.on(stream), condensed_tree.get_parents(), condensed_tree.get_parents() + (condensed_tree.get_n_edges()), condensed_tree.get_lambdas(), @@ -405,19 +406,55 @@ void excess_of_mass(const raft::handle_t &handle, } } -template -void get_stability_scores() { +template +struct ReduceSizes { + + ReduceSizes(value_idx *cluster_sizes_): cluster_sizes(cluster_sizes_){} + __host__ __device__ value_idx operator()(value_idx v) { + atomicAdd(cluster_sizes+v, 1); + } + + private: + value_idx *cluster_sizes; +}; + +template +void get_stability_scores(const raft::handle_t &handle, + const value_idx *labels, + const value_t *stability, + const value_idx *clusters, + value_idx n_clusters, + value_t max_lambda, + value_idx n_leaves, + value_t *result) { + + rmm::device_uvector cluster_sizes(n_clusters, handle.get_stream()); // TODO: Perform segmented reduction to compute cluster_size + thrust::for_each(thrust::cuda::par.on(handle.get_stream()), labels, labels+n_leaves, + ReduceSizes()); + // TODO: Embarassingly parallel + + auto enumeration = thrust::make_zip_iterator(thrust::make_tuple(clusters, cluster_sizes.data())); + thrust::transform(thrust::cuda::par.on(handle.get_stream()), + enumeration, enumeration+n_clusters, + result, [=] __device__ (thrust::tuple tup) { + value_idx size = thrust::get<1>(tup); + value_idx c = thrust::get<0>(tup); + + bool expr = max_lambda == std::numeric_limits::max() || + max_lambda == 0.0 || + size == 0; + return (!expr * (stability[c] / size * max_lambda)) + (expr * 1.0); + }); } template void do_labelling() { - - // TODO: union find is constructed on host } - + // TODO: This can be done efficiently on host. +} template void get_probabilities() { From bbe1ffe1552a9b2814cbe89f0ef292cc388bc371 Mon Sep 17 00:00:00 2001 From: divyegala Date: Fri, 23 Apr 2021 12:45:29 -0700 Subject: [PATCH 018/177] tests building --- cpp/CMakeLists.txt | 1 + cpp/bench/sg/linkage.cu | 1 + cpp/cmake/Dependencies.cmake | 4 +- cpp/include/cuml/cluster/hdbscan.hpp | 19 +- cpp/include/cuml/cluster/linkage.hpp | 1 + cpp/src/hdbscan/detail/tree_kernels.cuh | 101 ++-- cpp/src/hdbscan/hdbscan.cu | 20 +- cpp/src/hdbscan/reachability.cuh | 2 - cpp/src/hdbscan/runner.h | 115 ++--- cpp/src/hdbscan/tree.cuh | 228 +++++---- cpp/src_prims/label/classlabels.cuh | 13 +- cpp/test/CMakeLists.txt | 1 + cpp/test/sg/hdbscan_test.cu | 603 ++++++++++++++++++++++++ 13 files changed, 857 insertions(+), 252 deletions(-) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index d8c72d8578..dd372b5421 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -422,6 +422,7 @@ if(BUILD_CUML_CPP_LIBRARY) src/glm/glm.cu src/genetic/genetic.cu src/genetic/node.cu + src/hdbscan/hdbscan.cu src/holtwinters/holtwinters.cu src/kmeans/kmeans.cu src/knn/knn.cu diff --git a/cpp/bench/sg/linkage.cu b/cpp/bench/sg/linkage.cu index 5915f33c36..943b52e017 100644 --- a/cpp/bench/sg/linkage.cu +++ b/cpp/bench/sg/linkage.cu @@ -19,6 +19,7 @@ #include #include #include +#include #include #include "benchmark.cuh" diff --git a/cpp/cmake/Dependencies.cmake b/cpp/cmake/Dependencies.cmake index d6beff799e..9f64e2364b 100644 --- a/cpp/cmake/Dependencies.cmake +++ b/cpp/cmake/Dependencies.cmake @@ -38,8 +38,8 @@ else(DEFINED ENV{RAFT_PATH}) set(RAFT_DIR ${CMAKE_CURRENT_BINARY_DIR}/raft CACHE STRING "Path to RAFT repo") ExternalProject_Add(raft - GIT_REPOSITORY https://github.com/rapidsai/raft.git - GIT_TAG f0cd81fb49638eaddc9bf18998cc894f292bc293 + GIT_REPOSITORY https://github.com/cjnolet/raft.git + GIT_TAG fea-020-hdbscan PREFIX ${RAFT_DIR} CONFIGURE_COMMAND "" BUILD_COMMAND "" diff --git a/cpp/include/cuml/cluster/hdbscan.hpp b/cpp/include/cuml/cluster/hdbscan.hpp index 3a3c61ca42..b3840f79fa 100644 --- a/cpp/include/cuml/cluster/hdbscan.hpp +++ b/cpp/include/cuml/cluster/hdbscan.hpp @@ -17,14 +17,23 @@ #pragma once #include -#include #include +#include + namespace ML { -template -void hdbscan(const raft::handle_t &handle, value_t *X, size_t m, size_t n, - raft::distance::DistanceType metric, int k, int min_pts, - float alpha, hdbscan_output *out); +template +struct hdbscan_output { + int n_clusters; + value_idx *labels; + value_t *probabilities; +}; + +// template +void hdbscan(const raft::handle_t &handle, float *X, std::size_t m, + std::size_t n, raft::distance::DistanceType metric, int k, + int min_pts, int min_cluster_size, + hdbscan_output *out); }; // end namespace ML \ No newline at end of file diff --git a/cpp/include/cuml/cluster/linkage.hpp b/cpp/include/cuml/cluster/linkage.hpp index bbb1471581..dd9c6892df 100644 --- a/cpp/include/cuml/cluster/linkage.hpp +++ b/cpp/include/cuml/cluster/linkage.hpp @@ -18,6 +18,7 @@ #include #include +#include #include diff --git a/cpp/src/hdbscan/detail/tree_kernels.cuh b/cpp/src/hdbscan/detail/tree_kernels.cuh index d98d347ed1..f048b925d8 100644 --- a/cpp/src/hdbscan/detail/tree_kernels.cuh +++ b/cpp/src/hdbscan/detail/tree_kernels.cuh @@ -23,10 +23,10 @@ namespace detail { template __device__ value_t get_lambda(value_idx node, value_idx num_points, - value_t *deltas) { - value_t delta = deltas[node - num_points]; - if (delta > 0.0) return 1.0 / delta; - return std::numeric_limits::max(); + value_t *deltas) { + value_t delta = deltas[node - num_points]; + if (delta > 0.0) return 1.0 / delta; + return std::numeric_limits::max(); } /** @@ -46,94 +46,89 @@ __device__ value_t get_lambda(value_idx node, value_idx num_points, */ template __global__ void condense_hierarchy_kernel( - bool *frontier, value_idx *ignore, value_idx *relabel, - const value_idx *src, const value_idx *dst, const value_t *deltas, - const value_idx *sizes, int n_leaves, - int num_points, int min_cluster_size, - value_idx *out_parent, value_idx *out_child, - value_t *out_lambda, value_idx *out_count) { + bool *frontier, value_idx *ignore, value_idx *relabel, const value_idx *src, + const value_idx *dst, const value_t *deltas, const value_idx *sizes, + int n_leaves, int num_points, int min_cluster_size, value_idx *out_parent, + value_idx *out_child, value_t *out_lambda, value_idx *out_count) { + int node = blockDim.x * blockIdx.x + threadIdx.x; - int node = blockDim.x * blockIdx.x + threadIdx.x; + // If node is in frontier, flip frontier for children + if (node > n_leaves * 2 || !frontier[node]) return; - // If node is in frontier, flip frontier for children - if(node > n_leaves * 2 || !frontier[node]) - return; + frontier[node] = false; - frontier[node] = false; + // TODO: Check bounds + value_idx left_child = src[(node - num_points) * 2]; + value_idx right_child = dst[((node - num_points) * 2)]; - // TODO: Check bounds - value_idx left_child = src[(node - num_points) * 2]; - value_idx right_child = dst[((node - num_points) * 2)]; - - frontier[left_child] = true; - frontier[right_child] = true; + frontier[left_child] = true; + frontier[right_child] = true; - bool ignore_val = ignore[node]; - bool should_ignore = ignore_val > -1; + bool ignore_val = ignore[node]; + bool should_ignore = ignore_val > -1; - // If the current node is being ignored (e.g. > -1) then propagate the ignore - // to children, if any - ignore[left_child] = (should_ignore * ignore_val) + (!should_ignore * -1); - ignore[right_child] = (should_ignore * ignore_val) + (!should_ignore * -1); + // If the current node is being ignored (e.g. > -1) then propagate the ignore + // to children, if any + ignore[left_child] = (should_ignore * ignore_val) + (!should_ignore * -1); + ignore[right_child] = (should_ignore * ignore_val) + (!should_ignore * -1); - if (node < num_points) { + if (node < num_points) { out_parent[node] = relabel[should_ignore]; out_child[node] = node; out_lambda[node] = get_lambda(should_ignore, num_points, deltas); out_count[node] = 1; - } + } - // If node is not ignored and is not a leaf, condense its children - // if necessary - else if (!should_ignore and node >= num_points) { + // If node is not ignored and is not a leaf, condense its children + // if necessary + else if (!should_ignore and node >= num_points) { value_idx left_child = src[(node - num_points) * 2]; value_idx right_child = dst[((node - num_points) * 2)]; value_t lambda_value = get_lambda(node, num_points, deltas); int left_count = - left_child >= num_points ? sizes[left_child - num_points] : 1; + left_child >= num_points ? sizes[left_child - num_points] : 1; int right_count = - right_child >= num_points ? sizes[right_child - num_points] : 1; + right_child >= num_points ? sizes[right_child - num_points] : 1; // If both children are large enough, they should be relabeled and // included directly in the output hierarchy. if (left_count >= min_cluster_size && right_count >= min_cluster_size) { - relabel[left_child] = node; - out_parent[node] = relabel[node]; - out_child[node] = node; - out_lambda[node] = lambda_value; - out_count[node] = left_count; - - relabel[right_child] = node; - out_parent[node] = relabel[node]; - out_child[node] = node; - out_lambda[node] = lambda_value; - out_count[node] = left_count; + relabel[left_child] = node; + out_parent[node] = relabel[node]; + out_child[node] = node; + out_lambda[node] = lambda_value; + out_count[node] = left_count; + + relabel[right_child] = node; + out_parent[node] = relabel[node]; + out_child[node] = node; + out_lambda[node] = lambda_value; + out_count[node] = left_count; } // Consume left or right child as necessary bool left_child_too_small = left_count < min_cluster_size; bool right_child_too_small = right_count < min_cluster_size; ignore[left_child] = - (left_child_too_small * node) + (!left_child_too_small * -1); + (left_child_too_small * node) + (!left_child_too_small * -1); ignore[right_child] = - (right_child_too_small * node) + (!right_child_too_small * -1); + (right_child_too_small * node) + (!right_child_too_small * -1); // If only left or right child is too small, consume it and relabel the other // (to it can be its own cluster) bool only_left_child_too_small = - left_child_too_small && !right_child_too_small; + left_child_too_small && !right_child_too_small; bool only_right_child_too_small = - !left_child_too_small && right_child_too_small; + !left_child_too_small && right_child_too_small; relabel[right_child] = (only_left_child_too_small * relabel[node]) + - (!only_left_child_too_small * -1); + (!only_left_child_too_small * -1); relabel[left_child] = (only_right_child_too_small * relabel[node]) + - (!only_right_child_too_small * -1); - } + (!only_right_child_too_small * -1); + } } - }; // end namespace detail }; // end namespace Tree diff --git a/cpp/src/hdbscan/hdbscan.cu b/cpp/src/hdbscan/hdbscan.cu index b7df85d8b6..979fa5ce81 100644 --- a/cpp/src/hdbscan/hdbscan.cu +++ b/cpp/src/hdbscan/hdbscan.cu @@ -14,23 +14,23 @@ * limitations under the License. */ -#include -#include #include #include namespace ML { -template -void hdbscan(const raft::handle_t &handle, value_t *X, size_t m, size_t n, - raft::distance::DistanceType metric, int k, int min_pts, - float alpha, hdbscan_output *out) { - HDBSCAN::_fit(handle, X, m, n, metric, k, min_pts, alpha); +// template +void hdbscan(const raft::handle_t &handle, float *X, std::size_t m, + std::size_t n, raft::distance::DistanceType metric, int k, + int min_pts, int min_cluster_size, + hdbscan_output *out) { + HDBSCAN::_fit(handle, X, m, n, metric, k, min_pts, + min_cluster_size, out); } -void hdbscan(const raft::handle_t &handle, const float *X, size_t m, size_t n, - raft::distance::DistanceType metric, int k, int min_pts, - float alpha, hdbscan_output *out); +// void hdbscan(const raft::handle_t &handle, const float *X, int m, int n, +// raft::distance::DistanceType metric, int k, int min_pts, +// hdbscan_output *out); }; // end namespace ML \ No newline at end of file diff --git a/cpp/src/hdbscan/reachability.cuh b/cpp/src/hdbscan/reachability.cuh index 592686530c..52c9e2dce0 100644 --- a/cpp/src/hdbscan/reachability.cuh +++ b/cpp/src/hdbscan/reachability.cuh @@ -16,9 +16,7 @@ #pragma once -#include #include -#include #include #include diff --git a/cpp/src/hdbscan/runner.h b/cpp/src/hdbscan/runner.h index bef9060b72..5bf65d21de 100644 --- a/cpp/src/hdbscan/runner.h +++ b/cpp/src/hdbscan/runner.h @@ -24,8 +24,8 @@ #include #include -#include "tree.cuh" #include "reachability.cuh" +#include "tree.cuh" namespace ML { namespace HDBSCAN { @@ -44,64 +44,65 @@ struct MSTEpilogueReachability { private: value_t *core_distances; value_idx m; -} +}; template -void _fit(const raft::handle_t &handle, value_t *X, value_idx m, value_idx n, - raft::distance::DistanceType metric, int k, int min_pts, - float alpha, int min_cluster_size) { - auto d_alloc = handle.get_device_allocator(); - auto stream = handle.get_stream(); - - /** - * Mutual reachability graph - */ - - rmm::device_uvector mutual_reachability_graph_inds(k * m, stream); - rmm::device_uvector mutual_reachability_graph_dists(k * m, stream); - rmm::device_uvector core_dists(k * m, stream); - - Reachability::mutual_reachability_dists( - handle, X, m, n, metric, min_pts, k, mutual_reachability_graph_inds.data(), - mutual_reachability_graph_dists.data(), core_dists.data()); - - /** - * Construct MST sorted by weights - */ - rmm::device_uvector mst_rows(m - 1, stream); - rmm::device_uvector mst_cols(m - 1, stream); - rmm::device_uvector mst_data(m - 1, stream); - - // during knn graph connection - raft::hierarchy::detail::build_sorted_mst( - handle, X, mutual_reachability_graph_inds.data(), - mutual_reachability_graph_dists.data(), m, n, mst_rows, mst_cols, mst_data, - k * m, metric, 10, MSTEpilogueReachability()); - - /** - * Perform hierarchical labeling - */ - value_idx n_edges = m - 1; - - rmm::device_uvector out_src(n_edges, stream); - rmm::device_uvector out_dst(n_edges, stream); - rmm::device_uvector out_delta(n_edges, stream); - rmm::device_uvector out_size(n_edges, stream); - - raft::hierarchy::detail::build_dendrogram_host( - mst_rows.data(), mst_cols.data(), mst_data.data(), n_edges, out_src.data(), - out_dst.data(), out_delta.data(), out_size.data()); - - /** - * Condense branches of tree according to min cluster size - */ - Tree::CondensedHierarchy condensed_tree(m, stream); - condense_hierarchy(handle, out_src.data(), out_dst.data(), - out_delta.data(), out_size.data(), - min_cluster_size, m, condensed_tree); - - rmm::device_uvector stabilities(condensed_tree.get_n_clusters(), handle.get_stream()); - compute_stabilities(handle, condensed_tree, stabilities); +void _fit(const raft::handle_t &handle, value_t *X, std::size_t m, + std::size_t n, raft::distance::DistanceType metric, int k, + int min_pts, int min_cluster_size, + hdbscan_output *out) { + // auto d_alloc = handle.get_device_allocator(); + // auto stream = handle.get_stream(); + + // /** + // * Mutual reachability graph + // */ + + // rmm::device_uvector mutual_reachability_graph_inds(k * m, stream); + // rmm::device_uvector mutual_reachability_graph_dists(k * m, stream); + // rmm::device_uvector core_dists(k * m, stream); + + // Reachability::mutual_reachability_dists( + // handle, X, m, n, metric, min_pts, k, mutual_reachability_graph_inds.data(), + // mutual_reachability_graph_dists.data(), core_dists.data()); + + // /** + // * Construct MST sorted by weights + // */ + // rmm::device_uvector mst_rows(m - 1, stream); + // rmm::device_uvector mst_cols(m - 1, stream); + // rmm::device_uvector mst_data(m - 1, stream); + + // // during knn graph connection + // raft::hierarchy::detail::build_sorted_mst( + // handle, X, mutual_reachability_graph_inds.data(), + // mutual_reachability_graph_dists.data(), m, n, mst_rows, mst_cols, mst_data, + // k * m, metric, 10, MSTEpilogueReachability()); + + // /** + // * Perform hierarchical labeling + // */ + // value_idx n_edges = m - 1; + + // rmm::device_uvector out_src(n_edges, stream); + // rmm::device_uvector out_dst(n_edges, stream); + // rmm::device_uvector out_delta(n_edges, stream); + // rmm::device_uvector out_size(n_edges, stream); + + // raft::hierarchy::detail::build_dendrogram_host( + // mst_rows.data(), mst_cols.data(), mst_data.data(), n_edges, out_src.data(), + // out_dst.data(), out_delta.data(), out_size.data()); + + // /** + // * Condense branches of tree according to min cluster size + // */ + // Tree::CondensedHierarchy condensed_tree(m, stream); + // condense_hierarchy(handle, out_src.data(), out_dst.data(), + // out_delta.data(), out_size.data(), + // min_cluster_size, m, condensed_tree); + + // rmm::device_uvector stabilities(condensed_tree.get_n_clusters(), handle.get_stream()); + // compute_stabilities(handle, condensed_tree, stabilities); /** * Extract labels from stability */ diff --git a/cpp/src/hdbscan/tree.cuh b/cpp/src/hdbscan/tree.cuh index 6ce5347f78..206e919b68 100644 --- a/cpp/src/hdbscan/tree.cuh +++ b/cpp/src/hdbscan/tree.cuh @@ -18,7 +18,7 @@ #include "detail/tree_kernels.cuh" -#include +#include