From ab661c829ddb62d13b3a363f6c7d50c4f2535465 Mon Sep 17 00:00:00 2001 From: Grace Li Date: Wed, 14 Oct 2026 13:17:52 -0700 Subject: [PATCH 1/5] Add local cache map for min_incident_edge_map to facilitate look ups in updating edge chain parent, add draft function to write single-file metall datastore for cluster labels, add [[maybe_unused]] to lambdas and some minor clean up --- ...distributed_clustering_build_hierarchy.hpp | 65 ++-- .../distributed_clustering_hdbscan.hpp | 165 +++++++-- ...distributed_clustering_mst_contraction.hpp | 334 +++++++++++++----- .../run_distributed_hdbscan_clustering.cpp | 97 +++-- 4 files changed, 463 insertions(+), 198 deletions(-) diff --git a/src/clustering/distributed_clustering_build_hierarchy.hpp b/src/clustering/distributed_clustering_build_hierarchy.hpp index b6cb5a4..fc37900 100644 --- a/src/clustering/distributed_clustering_build_hierarchy.hpp +++ b/src/clustering/distributed_clustering_build_hierarchy.hpp @@ -41,12 +41,12 @@ void fill_alpha_edge_map( if (edge.chain_parent_edge_id > 0) { chain_supernode = edge.chain_supernode; } - auto set_info_lambda = [](const id_t &edge_id, alpha_edge_info &edge, - const distance_t &distance, - const supernode_t &chain_supernode) { - edge.distance = distance; - edge.chain_supernode = chain_supernode; - }; + auto set_info_lambda = + []([[maybe_unused]] const id_t &edge_id, alpha_edge_info &edge, + const distance_t &distance, const supernode_t &chain_supernode) { + edge.distance = distance; + edge.chain_supernode = chain_supernode; + }; alpha_edge_map.async_visit(edge_id, set_info_lambda, edge.distance, chain_supernode); } @@ -437,7 +437,7 @@ void fill_missing_non_root_chain_cluster_info( // Go through each non-root chain and for each cluster fill in the birth // distance auto set_birth_distance_lambda = - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain) { // The ith cluster has birth distance equal to the min edge distance of // the (i+1)th cluster @@ -463,12 +463,12 @@ void fill_missing_non_root_chain_cluster_info( &chain) { id_t parent_edge_id = chain.second.parent_edge_id; // Visit the parent edge in the alpha edge map and get its chain - auto visit_alpha_edge_lambda = [](const id_t &edge_id, + auto visit_alpha_edge_lambda = []([[maybe_unused]] const id_t &edge_id, const alpha_edge_info &edge_info, const supernode_t &chain_name, auto chain_map_ptr) { auto set_parent_chain_and_top_birth_dist = - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain, const supernode_t &parent_chain, @@ -511,7 +511,7 @@ void fill_missing_non_root_chain_cluster_info( } auto set_bottom_children_lambda = - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain, const std::array &children) { @@ -537,7 +537,7 @@ void fill_missing_non_root_chain_cluster_info( } auto set_child_lambda = - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain, const id_t &cluster_edge_id, const supernode_t &child) { @@ -588,24 +588,25 @@ void fill_missing_leaf_cluster_info( const supernode_t &cluster_name, full_leaf_cluster_info &cluster_info) { // Visit the parent edge in the alpha edge map and get its chain - auto visit_alpha_edge_lambda = - [](const id_t &edge_id, const alpha_edge_info &edge_info, - const supernode_t &cluster_name, auto leaf_cluster_map_ptr) { - auto set_parent_chain_and_top_birth_dist = - [](const supernode_t &cluster_name, - full_leaf_cluster_info &cluster_info, - const supernode_t &parent_chain, - const distance_t &parent_edge_distance) { - // Set the chain parent - cluster_info.parent_chain = parent_chain; - - // Set the birth distance for the last cluster in the chain - cluster_info.birth_distance = parent_edge_distance; - }; - leaf_cluster_map_ptr->async_visit( - cluster_name, set_parent_chain_and_top_birth_dist, - edge_info.chain_supernode, edge_info.distance); - }; + auto visit_alpha_edge_lambda = []([[maybe_unused]] const id_t &edge_id, + const alpha_edge_info &edge_info, + const supernode_t &cluster_name, + auto leaf_cluster_map_ptr) { + auto set_parent_chain_and_top_birth_dist = + []([[maybe_unused]] const supernode_t &cluster_name, + full_leaf_cluster_info &cluster_info, + const supernode_t &parent_chain, + const distance_t &parent_edge_distance) { + // Set the chain parent + cluster_info.parent_chain = parent_chain; + + // Set the birth distance for the last cluster in the chain + cluster_info.birth_distance = parent_edge_distance; + }; + leaf_cluster_map_ptr->async_visit( + cluster_name, set_parent_chain_and_top_birth_dist, + edge_info.chain_supernode, edge_info.distance); + }; alpha_edge_map.async_visit(cluster_info.parent_edge_id, visit_alpha_edge_lambda, cluster_name, leaf_cluster_map_ptr); @@ -702,7 +703,7 @@ supernode_t fill_missing_root_chain_cluster_info( [&root_chain_cluster_edges_map, full_root_chain_cluster_array_ptr]( const id_t &index, std::pair &value) { auto visit_cluster_edge_map = - [](const id_t &cluster_edge_id, + []([[maybe_unused]] const id_t &cluster_edge_id, const std::vector &edges, const id_t &array_index, ygm::ygm_ptrasync_visit( array_index, - [](const id_t &index, + []([[maybe_unused]] const id_t &index, std::pair &value, const id_t &num_points_added, const distance_t @@ -787,7 +788,7 @@ supernode_t fill_missing_root_chain_cluster_info( // Set the child in the root chain cluster array full_root_chain_cluster_array_ptr->async_visit( array_index, - [](const id_t &index, + []([[maybe_unused]] const id_t &index, std::pair &value, const supernode_t &child) { value.second.child = child; }, child); diff --git a/src/clustering/distributed_clustering_hdbscan.hpp b/src/clustering/distributed_clustering_hdbscan.hpp index 5fa6bd7..85909c7 100644 --- a/src/clustering/distributed_clustering_hdbscan.hpp +++ b/src/clustering/distributed_clustering_hdbscan.hpp @@ -18,8 +18,20 @@ #include "distributed_clustering_types.hpp" +#include +#include +#include + namespace clams::clustering { +// Cluster ID table that can be stored in a Metall datastore. +// This table also can be used as a regular in-memory table thanks to the +// scoped_fallback_allocator_type. +using pm_cluster_id_table_t = metall::container::unordered_flat_map< + id_t, id_t, std::hash, std::equal_to, + metall::manager::scoped_fallback_allocator_type< + std::pair>>; + /** * @brief Calculate cluster size, stability, and stability traversing up * from its child cluster information for a chain once all children have been @@ -48,8 +60,7 @@ std::size_t calculate_chain_cluster_size_stability_from_children( std::pair, full_chain_info>>> chain_map_ptr, ygm::ygm_ptr> - leaf_cluster_map_ptr, - supernode_t chain_name = BLANK_SUPERNODE) { + leaf_cluster_map_ptr) { std::size_t chain_num_valid_clusters = 0; // Keep track of the chain's total size and stability traversing up so far @@ -79,15 +90,15 @@ std::size_t calculate_chain_cluster_size_stability_from_children( chain_num_valid_clusters += 2; auto set_top_chain_cluster_as_valid = - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain) { auto top_cluster_it = chain.first.rbegin(); top_cluster_it->second.valid_cluster = true; }; auto set_leaf_cluster_as_valid = - [](const supernode_t &cluster_name, - full_leaf_cluster_info &cluster_info) { + []([[maybe_unused]] const supernode_t &cluster_name, + full_leaf_cluster_info &cluster_info) { cluster_info.valid_cluster = true; }; @@ -305,7 +316,8 @@ std::size_t traverse_up_cluster_hierarchy_until_root_chain( // Lambda to send child cluster information to root chain cluster static auto send_child_to_root_chain_cluster = - [](const id_t &cluster_edge_id, root_chain_cluster_info &cluster_info, + []([[maybe_unused]] const id_t &cluster_edge_id, + root_chain_cluster_info &cluster_info, const supernode_t &child_chain_name, id_t child_size, distance_t child_stability, distance_t child_stability_traversing_up) { if (cluster_info.child == child_chain_name) { @@ -380,7 +392,7 @@ std::size_t traverse_up_cluster_hierarchy_until_root_chain( id_t chain_num_valid_clusters = calculate_chain_cluster_size_stability_from_children( chain.first, chain.second, min_cluster_size, chain_map_ptr, - leaf_cluster_map_ptr, chain_name); + leaf_cluster_map_ptr); num_valid_clusters += chain_num_valid_clusters; auto top_cluster_it = chain.first.rbegin(); @@ -526,7 +538,7 @@ void make_sure_root_chain_bottom_is_valid_and_get_second_child( if (comm.rank() == 0) { chain_map.async_visit( root_chain_second_child.name, - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain) { auto top_cluster_it = chain.first.rbegin(); @@ -614,7 +626,7 @@ void make_sure_root_chain_bottom_is_valid_and_get_second_child( << root_chain_offshoot_supernode << std::endl; // debug static auto set_new_parent_lambda = - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain, const supernode_t &parent_chain) { chain.second.parent_chain = parent_chain; @@ -636,7 +648,7 @@ void make_sure_root_chain_bottom_is_valid_and_get_second_child( chain_map.async_visit( root_chain_offshoot_supernode, - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain, const full_chain_info &chain_info) { chain.second = chain_info; }, @@ -646,11 +658,11 @@ void make_sure_root_chain_bottom_is_valid_and_get_second_child( root_chain_offshoot_supernode); auto get_first_child_lambda = - [](const id_t &cluster_edge_id, + []([[maybe_unused]] const id_t &cluster_edge_id, const root_chain_cluster_info &cluster_info, auto chain_map_ptr) { chain_map_ptr->async_visit( root_chain_offshoot_supernode, - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain, const supernode_t &child) { @@ -695,7 +707,7 @@ void make_sure_root_chain_bottom_is_valid_and_get_second_child( .birth_distance = 1.0 / root_chain_cluster_info.lambda_birth, }; auto add_cluster_to_chain_map_lambda = - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain, const id_t &cluster_edge_id, @@ -727,7 +739,7 @@ void make_sure_root_chain_bottom_is_valid_and_get_second_child( std::vector &edges, auto chain_map_ptr) { chain_map_ptr->async_visit( root_chain_offshoot_supernode, - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain, const id_t &cluster_edge_id, @@ -764,7 +776,7 @@ void make_sure_root_chain_bottom_is_valid_and_get_second_child( if (comm.rank() == 0) { chain_map.async_visit( root_chain_offshoot_supernode, - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain, ygm::ygm_ptr, full_chain_info> &chain) { auto top_cluster_it = chain.first.rbegin(); @@ -825,7 +837,7 @@ void make_sure_root_chain_bottom_is_valid_and_get_second_child( if (comm.rank() == 0) { chain_map.async_visit( root_chain_second_child.name, - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain) { auto it = chain.first.rbegin(); @@ -1221,17 +1233,19 @@ void merge_invalid_root_chain_clusters( // Set the children of the post-merged root chain clusters as valid auto set_child_cluster_as_valid_lambda = [leaf_cluster_map_ptr, chain_map_ptr]( - const id_t &cluster_edge_id, root_chain_cluster_info &cluster_info) { + [[maybe_unused]] const id_t &cluster_edge_id, + root_chain_cluster_info &cluster_info) { if (cluster_info.child.second == 1) { leaf_cluster_map_ptr->async_visit( - cluster_info.child, [](const supernode_t &cluster_name, - full_leaf_cluster_info &cluster_info) { + cluster_info.child, + []([[maybe_unused]] const supernode_t &cluster_name, + full_leaf_cluster_info &cluster_info) { cluster_info.valid_cluster = true; }); } else { chain_map_ptr->async_visit( cluster_info.child, - [](const supernode_t &chain_name, + []([[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain) { auto it = chain.first.rbegin(); @@ -1268,8 +1282,7 @@ std::vector> calculate_root_chain_size_stability( ygm::container::array> &root_chain_array, - extra_child_cluster_info &root_chain_second_child, ygm::comm &comm, - const supernode_t root_chain_supernode) { + extra_child_cluster_info &root_chain_second_child, ygm::comm &comm) { // cluster size = sum of child cluster sizes + num points added // stability traversing up = max{stability, sum child stabilities // traversing up} @@ -1603,8 +1616,9 @@ std::vector traverse_down_hierarchy_and_select_clusters( // Select the child if it's a leaf cluster if (cluster_info.child.second == 1) { leaf_cluster_map_ptr->async_visit( - cluster_info.child, [](const supernode_t &cluster_name, - full_leaf_cluster_info &cluster_info) { + cluster_info.child, + []([[maybe_unused]] const supernode_t &cluster_name, + full_leaf_cluster_info &cluster_info) { cluster_info.selected = true; }); local_selected_clusters.push_back( @@ -1628,7 +1642,7 @@ std::vector traverse_down_hierarchy_and_select_clusters( // Select the child if it's a leaf cluster if (child.second == 1) { leaf_cluster_map_ptr->async_visit( - child, [](const supernode_t &cluster_name, + child, []([[maybe_unused]] const supernode_t &cluster_name, full_leaf_cluster_info &cluster_info) { cluster_info.selected = true; }); @@ -1657,8 +1671,9 @@ std::vector traverse_down_hierarchy_and_select_clusters( // away invalid root chain clusters if (cluster_info.child.second == 1) { leaf_cluster_map_ptr->async_visit( - cluster_info.child, [](const supernode_t &cluster_name, - full_leaf_cluster_info &cluster_info) { + cluster_info.child, + []([[maybe_unused]] const supernode_t &cluster_name, + full_leaf_cluster_info &cluster_info) { cluster_info.selected = true; }); local_selected_clusters.push_back( @@ -1681,8 +1696,8 @@ std::vector traverse_down_hierarchy_and_select_clusters( if (root_chain_second_child.name.second == 1) { leaf_cluster_map_ptr->async_visit( root_chain_second_child.name, - [](const supernode_t &cluster_name, - full_leaf_cluster_info &cluster_info) { + []([[maybe_unused]] const supernode_t &cluster_name, + full_leaf_cluster_info &cluster_info) { cluster_info.selected = true; }); local_selected_clusters.push_back( @@ -1783,8 +1798,8 @@ void assign_points_cluster_ids( // Goes through all edges in a leaf cluster and labels their points with the // provided cluster id static auto label_leaf_cluster_edge_endpoints_lambda = - [](const supernode_t &cluster_name, - const full_leaf_cluster_info &cluster_info, + []([[maybe_unused]] const supernode_t &cluster_name, + const full_leaf_cluster_info &cluster_info, const cluster_id_t &cluster_id, auto edge_endpoints_array_ptr, auto point_to_cluster_id_map_ptr) { for (const edge_id_with_dist_t &edge : cluster_info.edges) { @@ -1797,7 +1812,8 @@ void assign_points_cluster_ids( // Goes through all edges added to a root chain cluster and labels their // points with the provided cluster id static auto label_root_chain_cluster_edge_endpoints_lambda = - [](const id_t &cluster_edge_id, std::vector &edges, + []([[maybe_unused]] const id_t &cluster_edge_id, + std::vector &edges, const cluster_id_t &cluster_id, auto edge_endpoints_array_ptr, auto point_to_cluster_id_map_ptr) { for (const edge_id_with_dist_t &edge : edges) { @@ -1816,8 +1832,8 @@ void assign_points_cluster_ids( ygm::ygm_ptr, full_chain_info>>> - chain_map_ptr, - const supernode_t &chain_name, + chain_map_ptr, + [[maybe_unused]] const supernode_t &chain_name, std::pair, full_chain_info> &chain, const id_t &starting_cluster_edge_id, const cluster_id_t &cluster_id, ygm::ygm_ptr> @@ -1969,8 +1985,85 @@ void assign_points_cluster_ids( } } comm.barrier(); +} - return; +void write_cluster_labels_one_txt_file_per_rank( + ygm::container::map &point_to_cluster_id_map, + std::filesystem::path output_path) { + ygm::comm &comm = point_to_cluster_id_map.comm(); + + int num_files = comm.size(); + if (comm.rank() == 0) { + std::filesystem::create_directories(output_path); + } + comm.barrier(); + + // Output points and cluster labels + std::string label_file_name = std::to_string(comm.rank()) + ".txt"; + std::filesystem::path label_file_path = output_path / label_file_name; + std::ofstream label_ofs(label_file_path); + if (!label_ofs.is_open()) { + std::cout << "Error! Unable to open label file: " + << label_file_path.string() << std::endl; + } + + label_ofs << "# point_id\tcluster_id\n"; + for (const auto &[point_id, cluster_id] : point_to_cluster_id_map) { + label_ofs << point_id << "\t" << cluster_id << "\n"; + } + label_ofs.flush(); + label_ofs.close(); + comm.barrier(); +} + +void write_cluster_labels_single_pm( + ygm::container::map &point_to_cluster_id_map, + const std::string &output_filename) { + ygm::comm &comm = point_to_cluster_id_map.comm(); + + // Get local point ids and cluster labels + std::vector point_ids; + std::vector cluster_labels; + point_ids.reserve(point_to_cluster_id_map.local_size()); + cluster_labels.reserve(point_to_cluster_id_map.local_size()); + + for (const auto &[id, label] : point_to_cluster_id_map) { + point_ids.push_back(id); + cluster_labels.push_back(label); + } + + const size_t num_points = point_to_cluster_id_map.size(); + + MPI_Comm mpi_comm = comm.get_mpi_comm(); + if (comm.rank() == 0) { + metall::manager manager(metall::create_only, output_filename); + auto cluster_id_table = manager.construct( + metall::unique_instance)(manager.get_allocator()); + cluster_id_table->reserve(num_points); + + // Write own edges + for (size_t i = 0; i < point_ids.size(); ++i) { + cluster_id_table->emplace(point_ids.at(i), cluster_labels.at(i)); + } + std::cout << "Wrote labels from rank 0" << std::endl; + + // Write edges for all other ranks + for (int r = 1; r < comm.size(); ++r) { + std::vector remote_ids = + comm.mpi_recv>(r, 0, mpi_comm); + std::vector remote_labels = + comm.mpi_recv>(r, 1, mpi_comm); + for (size_t i = 0; i < remote_ids.size(); ++i) { + cluster_id_table->emplace(remote_ids.at(i), remote_labels.at(i)); + } + std::cout << "Wrote labels from rank " << r << std::endl; + } + + } else { + comm.mpi_send(point_ids, 0, 0, mpi_comm); + comm.mpi_send(cluster_labels, 0, 1, mpi_comm); + } + comm.barrier(); } /** diff --git a/src/clustering/distributed_clustering_mst_contraction.hpp b/src/clustering/distributed_clustering_mst_contraction.hpp index e20da98..f4db88b 100644 --- a/src/clustering/distributed_clustering_mst_contraction.hpp +++ b/src/clustering/distributed_clustering_mst_contraction.hpp @@ -20,6 +20,8 @@ namespace clams::clustering { +uint32_t HDBSCAN_LOCAL_CACHE_SIZE_LIMIT = 100000; + /** * @brief Sort and process MST edges. * @@ -39,6 +41,8 @@ void sort_and_process_mst_edges_into_maps( ygm::container::array &edge_contraction_array) { ygm::comm &comm = edge_endpoints_array.comm(); + ygm::utility::timer timer; // debug + // Create a YGM array with the MST edges ygm::container::array>> edge_array(comm, mst_edge_vector); @@ -48,6 +52,10 @@ void sort_and_process_mst_edges_into_maps( edge_array.sort(); comm.barrier(); + comm.cout0() << "Time to put edges into YGM array and sort (s): " + << timer.elapsed() << std::endl; + timer.reset(); + // Fill in the edge endpoints and contraction map auto visit_sorted_array_lambda = [&edge_endpoints_array, &edge_contraction_array]( @@ -64,10 +72,14 @@ void sort_and_process_mst_edges_into_maps( }; edge_array.for_all(visit_sorted_array_lambda); comm.barrier(); + + comm.cout0() << "Time to fill YGM edge maps (s): " << timer.elapsed() + << std::endl; } /** - * @brief Fill the initial min_incident_edge_map from the edge_endpoints_array + * @brief Fill the initial min_incident_edge_map from the + edge_endpoints_array * * @param edge_endpoints_array YGM array of edge id -> (edge_endpoint1, * edge_endpoint2) @@ -78,29 +90,19 @@ void sort_and_process_mst_edges_into_maps( void fill_init_min_incident_edge_map( ygm::container::map &min_incident_edge_map, ygm::container::array> &edge_endpoints_array) { - ygm::comm &comm = min_incident_edge_map.comm(); + ygm::comm &comm = edge_endpoints_array.comm(); static auto check_incident_edge_lambda = []([[maybe_unused]] const id_t &supernode_rep, id_t &min_incident_edge_id, const id_t &edge_id) { if (edge_id < min_incident_edge_id) { min_incident_edge_id = edge_id; + } else if (min_incident_edge_id == 0) { + min_incident_edge_id = edge_id; } }; - // Fill the initial incidence map from the original MST edges - auto fill_init_min_incident_edge_map_supernode_reps = - [&min_incident_edge_map]([[maybe_unused]] const id_t &edge_id, - std::pair &edge_endpoints) { - min_incident_edge_map.async_insert(edge_endpoints.first, - std::numeric_limits::max()); - min_incident_edge_map.async_insert(edge_endpoints.second, - std::numeric_limits::max()); - }; - edge_endpoints_array.for_all(fill_init_min_incident_edge_map_supernode_reps); - comm.barrier(); - - // Fill the initial incidence map from the original MST edges + // Fill the initial incidence array from the original MST edges auto fill_init_min_incident_edge_map = [&min_incident_edge_map](const id_t &edge_id, std::pair &edge_endpoints) { @@ -111,6 +113,13 @@ void fill_init_min_incident_edge_map( }; edge_endpoints_array.for_all(fill_init_min_incident_edge_map); comm.barrier(); + + // Since map values initialize to 0, set the two endpoints of edge 0 to have + // min incident edge 0 at the end + if (comm.rank() == 0) { + edge_endpoints_array.async_visit(0, fill_init_min_incident_edge_map); + } + comm.barrier(); } /** @@ -127,8 +136,9 @@ void fill_init_min_incident_edge_map( * @param tree_components_djset YGM disjoint_set for getting connected * components after edge contraction. */ +template void contract_edges( - uint32_t _round, ygm::container::map &min_incident_edge_map, + uint32_t _round, ygm_container_type &min_incident_edge_map, ygm::container::array &edge_contraction_array, ygm::ygm_ptr> alpha_edge_map_ptr, ygm::container::disjoint_set &tree_components_djset) { @@ -138,6 +148,70 @@ void contract_edges( static uint32_t round; round = _round; + [[maybe_unused]] static auto add_dendrogram_child_lambda = + []([[maybe_unused]] const id_t &edge_id, alpha_edge_info &info, + supernode_t new_child) { + if (info.dendrogram_children[0] == BLANK_SUPERNODE) { + info.dendrogram_children[0] = new_child; + } else { + info.dendrogram_children[1] = new_child; + } + }; + + auto contract_edges_lambda = + [&edge_contraction_array, alpha_edge_map_ptr, tree_components_djset_ptr]( + const id_t &supernode_rep, const id_t &min_incident_edge_id) { + // The shortest edge incident to each supernode is contracted + + // Add the contracted edge to the disjoint set and set as + // contracted For an edge between supernodes (rep1, round-1) and + // (rep2, round-1), we add the rep1-rep2 edge to the disjoint set + auto add_edge_to_djset_lambda = []([[maybe_unused]] const id_t &edge_id, + edge_contraction_info &edge_info, + auto tree_components_djset_ptr) { + edge_info.contraction_round = round; + tree_components_djset_ptr->async_union( + edge_info.endpoint_supernode_reps.first, + edge_info.endpoint_supernode_reps.second); + }; + if (min_incident_edge_id != std::numeric_limits::max()) { + edge_contraction_array.async_visit(min_incident_edge_id, + add_edge_to_djset_lambda, + tree_components_djset_ptr); + } + // Set this supernode as the dendrogram child of the min incident edge + if (round >= 2) { + alpha_edge_map_ptr->async_visit( + min_incident_edge_id, add_dendrogram_child_lambda, + std::make_pair(supernode_rep, round - 1)); + } + }; + min_incident_edge_map.for_all(contract_edges_lambda); + comm.barrier(); +} + +/** + * @brief Find all MST edges to contract this round and add to disjoint set + * + * @param _round The round of MST contraction we're on. + * @param min_incident_edge_map YGM map of supernode representative -> + * min incident edge id. + * @param edge_contraction_array YGM array of edge id -> edge_contraction_info + * @param alpha_edge_map_ptr Pointer to YGM map of alpha edge id -> + * alpha_edge_info. Alpha edges are ones that are not contracted in the first + * round. As we contract alpha edges, we update their child supernodes in the + * dendrogram. + * @param tree_components_djset YGM disjoint_set for getting connected + * components after edge contraction. + */ +void contract_edges_first_round( + ygm::container::array &min_incident_edge_array, + ygm::container::array &edge_contraction_array, + ygm::ygm_ptr> alpha_edge_map_ptr, + ygm::container::disjoint_set &tree_components_djset) { + ygm::comm &comm = edge_contraction_array.comm(); + auto tree_components_djset_ptr = tree_components_djset.get_ygm_ptr(); + [[maybe_unused]] static auto add_dendrogram_child_lambda = []([[maybe_unused]] const id_t &edge_id, alpha_edge_info &info, supernode_t new_child) { @@ -151,7 +225,7 @@ void contract_edges( auto contract_edges_lambda = [&edge_contraction_array, alpha_edge_map_ptr, tree_components_djset_ptr]( [[maybe_unused]] const id_t &supernode_rep, - const id_t min_incident_edge_id) { + const id_t &min_incident_edge_id) { // The shortest edge incident to each supernode is contracted // Add the contracted edge to the disjoint set and set as @@ -160,26 +234,19 @@ void contract_edges( auto add_edge_to_djset_lambda = []([[maybe_unused]] const id_t &edge_id, edge_contraction_info &edge_info, auto tree_components_djset_ptr) { - edge_info.contraction_round = round; + edge_info.contraction_round = 1; tree_components_djset_ptr->async_union( edge_info.endpoint_supernode_reps.first, edge_info.endpoint_supernode_reps.second); }; - edge_contraction_array.async_visit(min_incident_edge_id, - add_edge_to_djset_lambda, - tree_components_djset_ptr); - - // Set this supernode as the dendrogram child of the min incident edge - if (round >= 2) { - alpha_edge_map_ptr->async_visit(min_incident_edge_id, - add_dendrogram_child_lambda, - std::make_pair(supernode_rep, round - 1)); + if (min_incident_edge_id != std::numeric_limits::max()) { + edge_contraction_array.async_visit(min_incident_edge_id, + add_edge_to_djset_lambda, + tree_components_djset_ptr); } }; - min_incident_edge_map.for_all(contract_edges_lambda); + min_incident_edge_array.for_all(contract_edges_lambda); comm.barrier(); - - return; } /** @@ -200,46 +267,136 @@ void update_edge_chain_parent_edge_id( ygm::comm &comm = edge_contraction_array.comm(); auto edge_contraction_array_ptr = edge_contraction_array.get_ygm_ptr(); - auto update_chain_parents_lambda = - [&min_incident_edge_map, round, edge_contraction_array_ptr]( - const id_t &edge_id, edge_contraction_info &edge_info) { - // Only see if we can update an edge's chain parent edge id if not - // assigned yet (if its 0) and if the edge is already contracted - if (edge_info.chain_parent_edge_id == 0 && - 0 < edge_info.contraction_round && - edge_info.contraction_round < round) { + static std::unordered_map local_min_incident_edge_map; + local_min_incident_edge_map.clear(); + + static std::unordered_set supernode_reps_with_visit_sent; + supernode_reps_with_visit_sent.clear(); + static std::unordered_map> + map_of_edges_awaiting_visit_return; + map_of_edges_awaiting_visit_return.clear(); + + auto update_chain_parents_lambda = [&min_incident_edge_map, round, + edge_contraction_array_ptr]( + const id_t &edge_id, + edge_contraction_info &edge_info) { + // Only see if we can update an edge's chain parent edge id if not + // assigned yet (if its 0) and if the edge is already contracted + if (edge_info.chain_parent_edge_id == 0 && + 0 < edge_info.contraction_round && + edge_info.contraction_round < round) { + // If we already have the min incident edge id stored locally, then + // directly try to update the chain parent + if (local_min_incident_edge_map.find(edge_info.chain_supernode.first) != + local_min_incident_edge_map.end()) { + id_t min_incident_edge_id = + local_min_incident_edge_map.at(edge_info.chain_supernode.first); + if (min_incident_edge_id > edge_id) { + edge_info.chain_parent_edge_id = min_incident_edge_id; + } + } else { + // If we already sent out an async visit call for this chain + // supernode, then add this edge to the local map and await the + // returned info to update it later + if (supernode_reps_with_visit_sent.find( + edge_info.chain_supernode.first) != + supernode_reps_with_visit_sent.end()) { + map_of_edges_awaiting_visit_return[edge_info.chain_supernode.first] + .push_back(edge_id); + } + // Do an async visit to min incident edge map the first time we + // encounter this chain supernode locally + else { + supernode_reps_with_visit_sent.insert( + edge_info.chain_supernode.first); min_incident_edge_map.async_visit( edge_info.chain_supernode.first, []([[maybe_unused]] const id_t &supernode_rep, const id_t &min_incident_edge_id, const id_t &edge_id, auto edge_contraction_array_ptr) { - // If the supernode's min incident edge id is larger than edge - // id, we found the chain's parent edge id in the dendrogram - if (min_incident_edge_id > edge_id) { - edge_contraction_array_ptr->async_visit( - edge_id, - []([[maybe_unused]] const id_t &edge_id, - edge_contraction_info &edge_info, - id_t chain_parent_edge_id) { - edge_info.chain_parent_edge_id = chain_parent_edge_id; - }, - min_incident_edge_id); - } + // If the supernode's min incident edge id is larger than + // edge id, we found the chain's parent edge id in the + // dendrogram + edge_contraction_array_ptr->async_visit( + edge_id, + []([[maybe_unused]] const id_t &edge_id, + edge_contraction_info &edge_info, + id_t min_incident_edge_id, id_t supernode_rep) { + if (min_incident_edge_id > edge_id) { + edge_info.chain_parent_edge_id = min_incident_edge_id; + } + + if (local_min_incident_edge_map.size() < + HDBSCAN_LOCAL_CACHE_SIZE_LIMIT) { + local_min_incident_edge_map[supernode_rep] = + min_incident_edge_id; + } + }, + min_incident_edge_id, supernode_rep); }, edge_id, edge_contraction_array_ptr); } - }; + } + } + }; edge_contraction_array.for_all(update_chain_parents_lambda); comm.barrier(); - return; + // Go through the edges awaiting return messages and update them. If we + // didn't get a return message and stored value, then do the required + // async visits + for (auto &[supernode_rep, edge_id_vector] : + map_of_edges_awaiting_visit_return) { + // If we already have the min incident edge id stored locally, then + // directly try to update the chain parents + if (local_min_incident_edge_map.find(supernode_rep) != + local_min_incident_edge_map.end()) { + id_t min_incident_edge_id = local_min_incident_edge_map.at(supernode_rep); + for (id_t edge_id : edge_id_vector) { + if (min_incident_edge_id > edge_id) { + edge_contraction_array.local_visit( + edge_id, + []([[maybe_unused]] const id_t &edge_id, + edge_contraction_info &edge_info, id_t chain_parent_edge_id) { + edge_info.chain_parent_edge_id = chain_parent_edge_id; + }, + min_incident_edge_id); + } + } + } else { + for (id_t &edge_id : edge_id_vector) { + // Do an async visit to min incident edge map + min_incident_edge_map.async_visit( + supernode_rep, + []([[maybe_unused]] const id_t &supernode_rep, + const id_t &min_incident_edge_id, const id_t &edge_id, + auto edge_contraction_array_ptr) { + // If the supernode's min incident edge id is larger than + // edge id, we found the chain's parent edge id in the + // dendrogram + if (min_incident_edge_id > edge_id) { + edge_contraction_array_ptr->async_visit( + edge_id, + []([[maybe_unused]] const id_t &edge_id, + edge_contraction_info &edge_info, + id_t chain_parent_edge_id, id_t supernode_rep) { + edge_info.chain_parent_edge_id = chain_parent_edge_id; + }, + min_incident_edge_id, supernode_rep); + } + }, + edge_id, edge_contraction_array_ptr); + } + } + } + comm.barrier(); } /** - * @brief Updates the edge contraction map by updating the chain parent edge id - * of contracted edges e by looking up the current supernode they are in and - * seeing if the supernode's min incident edge has larger id than e. Makes a - * local copy of the YGM min_incident_edge_map on each rank to avoid + * @brief Updates the edge contraction map by updating the chain parent edge + * id of contracted edges e by looking up the current supernode they are in + * and seeing if the supernode's min incident edge has larger id than e. Makes + * a local copy of the YGM min_incident_edge_map on each rank to avoid * bottlenecks. This version should be used instead of * update_edge_chain_parent_edge_id when the min_incident_edge_map is small. * @@ -255,7 +412,7 @@ void update_edge_chain_parent_edge_id_from_local_map( ygm::comm &comm = edge_contraction_array.comm(); // Give all ranks a full copy of min_incident_edge_map - std::map local_min_incident_edge_map; + std::unordered_map local_min_incident_edge_map; min_incident_edge_map.for_all( [&local_min_incident_edge_map](const id_t &supernode_rep, const id_t &min_incident_edge_id) { @@ -263,10 +420,10 @@ void update_edge_chain_parent_edge_id_from_local_map( }); comm.barrier(); - static std::map min_incident_edge_std_map; + static std::unordered_map min_incident_edge_std_map; min_incident_edge_std_map.clear(); comm.async_bcast( - [](const std::map &local_map) { + [](const std::unordered_map &local_map) { min_incident_edge_std_map.insert(local_map.begin(), local_map.end()); }, local_min_incident_edge_map); @@ -291,20 +448,18 @@ void update_edge_chain_parent_edge_id_from_local_map( }; edge_contraction_array.for_all(update_chain_parents_lambda); comm.barrier(); - - return; } /** - * @brief Updates edge_contraction_array by updating the supernode endpoints of - * uncontracted edges and updates the current chain supernode of contracted + * @brief Updates edge_contraction_array by updating the supernode endpoints + * of uncontracted edges and updates the current chain supernode of contracted * edges if the chain parent edge id is not yet found. This version directly * visits the tree_conmponents_djset. * * @param _round The round of MST contraction we're on. * @param edge_contraction_array YGM array of edge id -> edge_contraction_info - * @param tree_components_djset YGM disjoint_set of connected components of the - * contracted tree. We use this to get the new supernode reps. + * @param tree_components_djset YGM disjoint_set of connected components of + * the contracted tree. We use this to get the new supernode reps. */ void update_edge_endpoints_and_chain_supernode( uint32_t _round, @@ -356,9 +511,9 @@ void update_edge_endpoints_and_chain_supernode( edge_id, edge_contraction_array_ptr); } - // If an edge is already contracted, see if it already has a known chain - // parent edge id, or if we should update the chain supernode to find - // the chain this edge belongs to + // If an edge is already contracted, see if it already has a known + // chain parent edge id, or if we should update the chain supernode to + // find the chain this edge belongs to else { // If we haven't found a chain_parent_edge_id > edge_id in the // dendrogram yet, then update the chain supernode @@ -394,23 +549,21 @@ void update_edge_endpoints_and_chain_supernode( edge_contraction_array.for_all( update_edge_endpoints_and_chain_supernode_lambda); comm.barrier(); - - return; } /** - * @brief Updates edge_contraction_array by updating the supernode endpoints of - * uncontracted edges and updates the current chain supernode of contracted + * @brief Updates edge_contraction_array by updating the supernode endpoints + * of uncontracted edges and updates the current chain supernode of contracted * edges if the chain parent edge id is not yet found. Makes a - * local copy of the tree_components_djset representatives on each rank to avoid - * bottlenecks. This version should be used instead of - * update_edge_endpoints_and_chain_supernode when the number of tree components - * is small. + * local copy of the tree_components_djset representatives on each rank to + * avoid bottlenecks. This version should be used instead of + * update_edge_endpoints_and_chain_supernode when the number of tree + * components is small. * * @param _round The round of MST contraction we're on. * @param edge_contraction_array YGM array of edge id -> edge_contraction_info - * @param tree_components_djset YGM disjoint_set of connected components of the - * contracted tree. We use this to get the new supernode reps. + * @param tree_components_djset YGM disjoint_set of connected components of + * the contracted tree. We use this to get the new supernode reps. */ void update_edge_endpoints_and_chain_supernode_from_local_map( uint32_t round, @@ -418,8 +571,8 @@ void update_edge_endpoints_and_chain_supernode_from_local_map( ygm::container::disjoint_set &tree_components_djset) { ygm::comm &comm = edge_contraction_array.comm(); - static std::map old_to_new_supernode_rep_map; - static std::map local_old_to_new_supernode_rep_map; + static std::unordered_map old_to_new_supernode_rep_map; + static std::unordered_map local_old_to_new_supernode_rep_map; old_to_new_supernode_rep_map.clear(); local_old_to_new_supernode_rep_map.clear(); @@ -429,8 +582,9 @@ void update_edge_endpoints_and_chain_supernode_from_local_map( local_old_to_new_supernode_rep_map[old_supernode_rep] = new_supernode_rep; }); comm.barrier(); + comm.async_bcast( - [](const std::map &local_map) { + [](const std::unordered_map &local_map) { old_to_new_supernode_rep_map.insert(local_map.begin(), local_map.end()); }, local_old_to_new_supernode_rep_map); @@ -453,9 +607,9 @@ void update_edge_endpoints_and_chain_supernode_from_local_map( old_to_new_supernode_rep_map.at( edge_info.endpoint_supernode_reps.second); } - // If an edge is already contracted, see if it already has a known chain - // parent edge id, or if we should update the chain supernode to find - // the chain this edge belongs to + // If an edge is already contracted, see if it already has a known + // chain parent edge id, or if we should update the chain supernode to + // find the chain this edge belongs to else { // If we haven't found a chain_parent_edge_id > edge_id in the // dendrogram yet, then update the chain supernode @@ -478,8 +632,6 @@ void update_edge_endpoints_and_chain_supernode_from_local_map( edge_contraction_array.for_all( update_edge_endpoints_and_chain_supernode_lambda); comm.barrier(); - - return; } /** @@ -488,8 +640,8 @@ void update_edge_endpoints_and_chain_supernode_from_local_map( * * @param min_incident_edge_map YGM map of supernode representative -> * min incident edge id. We will clear and refill this map. - * @param tree_components_djset YGM disjoint_set of connected components of the - * contracted tree. We use this to get the new supernode reps. + * @param tree_components_djset YGM disjoint_set of connected components of + * the contracted tree. We use this to get the new supernode reps. */ void reinit_min_incident_edge_map( ygm::container::map &min_incident_edge_map, @@ -507,8 +659,6 @@ void reinit_min_incident_edge_map( }; tree_components_djset.for_all(reinit_min_incident_edge_map_lambda); comm.barrier(); - - return; } /** @@ -519,8 +669,8 @@ void reinit_min_incident_edge_map( * @param edge_contraction_array YGM map of edge id -> edge contraction info. * Assumes that edges already have their supernode endpoints updated. * @param min_incident_edge_map YGM map of supernode representative -> - * min incident edge id. Assumes that the map already has the correct supernode - * reps as keys. + * min incident edge id. Assumes that the map already has the correct + * supernode reps as keys. */ void update_min_incident_edge_map( ygm::container::array &edge_contraction_array, diff --git a/src/clustering/run_distributed_hdbscan_clustering.cpp b/src/clustering/run_distributed_hdbscan_clustering.cpp index e407c6f..4664903 100644 --- a/src/clustering/run_distributed_hdbscan_clustering.cpp +++ b/src/clustering/run_distributed_hdbscan_clustering.cpp @@ -35,13 +35,14 @@ using namespace clams::clustering; -uint32_t LOCAL_COPY_NUM_SUPERNODES_THRESHOLD = 25000; +uint32_t HDBSCAN_LOCAL_COPY_NUM_SUPERNODES_THRESHOLD = 25000; struct option { std::filesystem::path mst_edges_path; uint32_t min_cluster_size{2}; std::filesystem::path cluster_ids_out_path; std::filesystem::path cluster_data_out_path; + bool write_metall_output{false}; bool metall_mst{false}; bool verbose{false}; bool read_mst_with_edge_ids{false}; @@ -53,14 +54,17 @@ void show_help() { << "<>\n" "Required arguments:\n" " -i Path to a file or a directory that contains input MST. \n" - " -o Path to clustering output directory. One text file with " - "lines of the form will be written per rank. " - "The filenames will be of the form rank.txt.\n" + " -o Path to clustering output directory or file. By default, " + "one text file with lines of the form will be " + "written per rank. The filenames will be of the form rank.txt. Use " + "the -s flag to write a Metall datastore instead.\n" " -m Min cluster size parameter in HDBSCAN.\n" "Optional arguments:\n" " -M If specified, input MST is Metall datastore.\n" " -e If specified, read MST file with edge ids included. This option " "only works for text file MST input, not Metall input.\n" + " -s If specified, write cluster labels to a single Metall " + "datastore located at the provided -o path.\n" " -c Path to detailed clustering data output directory. One " "csv file will be written per rank.\n" " -v Verbose printout.\n" @@ -77,7 +81,7 @@ std::pair> parse_option(int argc, char *argv[], bool show_help = false; std::vector unknown_opts; int opt_char; - while ((opt_char = getopt(argc, argv, "i:o:m:c:n:Mevh")) != -1) { + while ((opt_char = getopt(argc, argv, "i:o:m:c:Msevh")) != -1) { switch (opt_char) { case 'i': opt.mst_edges_path = std::filesystem::path(optarg); @@ -97,6 +101,9 @@ std::pair> parse_option(int argc, char *argv[], case 'e': opt.read_mst_with_edge_ids = true; break; + case 's': + opt.write_metall_output = true; + break; case 'v': opt.verbose = true; break; @@ -186,6 +193,8 @@ int main(int argc, char *argv[]) { mst_edge_vector; static id_t local_num_edges_read; local_num_edges_read = 0; + static id_t local_max_point_id; + local_max_point_id = 0; if (opt.metall_mst) { if (world.rank() == 0) { @@ -218,6 +227,13 @@ int main(int argc, char *argv[]) { std::make_pair((*input_mst_edges)[i].ids[0], (*input_mst_edges)[i].ids[1])); mst_edge_vector.push_back(array_item); + + if ((*input_mst_edges)[i].ids[0] > local_max_point_id) { + local_max_point_id = (*input_mst_edges)[i].ids[0]; + } + if ((*input_mst_edges)[i].ids[1] > local_max_point_id) { + local_max_point_id = (*input_mst_edges)[i].ids[1]; + } } local_num_edges_read = mst_edge_vector.size(); @@ -245,7 +261,12 @@ int main(int argc, char *argv[]) { distance_t dist; ss >> edge_id >> node1 >> node2 >> dist; ++local_num_edges_read; - + if (node1 > local_max_point_id) { + local_max_point_id = node1; + } + if (node2 > local_max_point_id) { + local_max_point_id = node2; + } } catch (...) { std::cout << "Error reading mst line: " << line << std::endl; } @@ -279,6 +300,12 @@ int main(int argc, char *argv[]) { std::pair> array_item = std::make_pair(dist, std::make_pair(node1, node2)); mst_edge_vector.push_back(array_item); + if (node1 > local_max_point_id) { + local_max_point_id = node1; + } + if (node2 > local_max_point_id) { + local_max_point_id = node2; + } } catch (...) { std::cout << "Error reading mst line: " << line << std::endl; } @@ -301,9 +328,10 @@ int main(int argc, char *argv[]) { sw_step.reset(); id_t total_num_edges_read = ygm::sum(local_num_edges_read, world); + id_t max_point_id = ygm::max(local_max_point_id, world); - // Initialize edge map of basic endpoints - // and edge contraction map with edge_contraction_info + // Initialize edge array of basic endpoints + // and edge contraction array with edge_contraction_info ygm::container::array> edge_endpoints_array( world, total_num_edges_read); ygm::container::array edge_contraction_array( @@ -363,6 +391,7 @@ int main(int argc, char *argv[]) { ss << " Number of MST edges: " << edge_endpoints_array.size(); if (world.rank() == 0) { spdlog::info(ss.str()); + spdlog::info(" Max point id: {}", max_point_id); spdlog::info(" Time to ingest MST edges (s): {:.3f}", sw_round); } @@ -423,8 +452,10 @@ int main(int argc, char *argv[]) { // Find all edges to contract and add to disjoint set sw_step.reset(); + id_t num_nodes_at_start_of_round; contract_edges(round, min_incident_edge_map, edge_contraction_array, alpha_edge_map_ptr, tree_components_djset); + num_nodes_at_start_of_round = min_incident_edge_map.size(); // Count how many edges were contracted this round id_t num_edges_contracted = 0; @@ -440,11 +471,10 @@ int main(int argc, char *argv[]) { contract_edges_time += sw_step.elapsed().count(); sw_step.reset(); - // Update the succ dendrogram parents for edges that were contracted // before this round if (round >= 2) { if (min_incident_edge_map.size() <= - LOCAL_COPY_NUM_SUPERNODES_THRESHOLD) { + HDBSCAN_LOCAL_COPY_NUM_SUPERNODES_THRESHOLD) { update_edge_chain_parent_edge_id_from_local_map( round, edge_contraction_array, min_incident_edge_map); } else { @@ -459,7 +489,6 @@ int main(int argc, char *argv[]) { // clear the min incident edge map since we're done with it for this // round - id_t num_nodes_at_start_of_round = min_incident_edge_map.size(); min_incident_edge_map.clear(); world.barrier(); @@ -506,7 +535,7 @@ int main(int argc, char *argv[]) { // we have few enough new supernodes, create a local copy of // tree_components_djset to update edge info if (tree_components_djset.num_sets() <= - LOCAL_COPY_NUM_SUPERNODES_THRESHOLD) { + HDBSCAN_LOCAL_COPY_NUM_SUPERNODES_THRESHOLD) { update_edge_endpoints_and_chain_supernode_from_local_map( round, edge_contraction_array, tree_components_djset); } @@ -675,7 +704,6 @@ int main(int argc, char *argv[]) { sw_step.reset(); fill_missing_leaf_cluster_info(leaf_cluster_map, alpha_edge_map); fill_missing_non_root_chain_cluster_info(chain_map, alpha_edge_map); - fill_initial_cluster_info_time += sw_step.elapsed().count(); // Erase alpha edges from the map if they are not in the root chain @@ -774,7 +802,7 @@ int main(int argc, char *argv[]) { // root chain clusters by their alpha edge id in phase 3 full_root_chain_cluster_array.for_all( [&root_chain_cluster_map]( - const id_t &index, + [[maybe_unused]] const id_t &index, const std::pair &value) { root_chain_cluster_map.async_insert(value.first, value.second); }); @@ -1014,8 +1042,7 @@ int main(int argc, char *argv[]) { { local_possible_clusters_for_selection = calculate_root_chain_size_stability(root_chain_array, - root_chain_second_child, world, - root_chain_supernode); + root_chain_second_child, world); // Update the cluster map with new root chain info update_root_chain_cluster_info(root_chain_array, @@ -1228,36 +1255,30 @@ int main(int argc, char *argv[]) { { sw_step.reset(); - id_t num_points_clustered = 0; - - int num_files = world.size(); - if (world.rank() == 0) { - spdlog::info("Writing {} files containing cluster labels to: {}", - num_files, opt.cluster_ids_out_path.string()); - std::filesystem::create_directories(opt.cluster_ids_out_path); - } - world.barrier(); - - // Output points and cluster labels - std::string label_file_name = std::to_string(world.rank()) + ".txt"; - std::filesystem::path label_file_path = - opt.cluster_ids_out_path / label_file_name; - std::ofstream label_ofs(label_file_path); - if (!label_ofs.is_open()) { - spdlog::warn("Unable to open label file: {}", label_file_path.string()); + if (opt.write_metall_output) { + if (world.rank() == 0) { + spdlog::info("Writing cluster labels to Metall datastore at: {}", + opt.cluster_ids_out_path.string()); + } + write_cluster_labels_single_pm(point_to_cluster_id_map, + opt.cluster_ids_out_path.string()); + } else { + int num_files = world.size(); + if (world.rank() == 0) { + spdlog::info("Writing {} files containing cluster labels to: {}", + num_files, opt.cluster_ids_out_path.string()); + } + write_cluster_labels_one_txt_file_per_rank(point_to_cluster_id_map, + opt.cluster_ids_out_path); } - label_ofs << "# point_id\tcluster_id\n"; + id_t num_points_clustered = 0; for (const auto &[point_id, cluster_id] : point_to_cluster_id_map) { - label_ofs << point_id << "\t" << cluster_id << "\n"; if (cluster_id >= 0) { ++num_points_clustered; } } - label_ofs.flush(); - label_ofs.close(); world.barrier(); - id_t total_num_points_clustered = ygm::sum(num_points_clustered, world); if (world.rank() == 0) { spdlog::info("Number of points clustered: {}", From 9eb32c26d08befa6aad7dc95cb21ae2097392c58 Mon Sep 17 00:00:00 2001 From: Grace Li Date: Wed, 14 Oct 2026 13:35:44 -0700 Subject: [PATCH 2/5] Remove extra debug timer printouts in ingesting edges --- .../distributed_clustering_mst_contraction.hpp | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/clustering/distributed_clustering_mst_contraction.hpp b/src/clustering/distributed_clustering_mst_contraction.hpp index f4db88b..eb5df40 100644 --- a/src/clustering/distributed_clustering_mst_contraction.hpp +++ b/src/clustering/distributed_clustering_mst_contraction.hpp @@ -41,8 +41,6 @@ void sort_and_process_mst_edges_into_maps( ygm::container::array &edge_contraction_array) { ygm::comm &comm = edge_endpoints_array.comm(); - ygm::utility::timer timer; // debug - // Create a YGM array with the MST edges ygm::container::array>> edge_array(comm, mst_edge_vector); @@ -52,10 +50,6 @@ void sort_and_process_mst_edges_into_maps( edge_array.sort(); comm.barrier(); - comm.cout0() << "Time to put edges into YGM array and sort (s): " - << timer.elapsed() << std::endl; - timer.reset(); - // Fill in the edge endpoints and contraction map auto visit_sorted_array_lambda = [&edge_endpoints_array, &edge_contraction_array]( @@ -72,9 +66,6 @@ void sort_and_process_mst_edges_into_maps( }; edge_array.for_all(visit_sorted_array_lambda); comm.barrier(); - - comm.cout0() << "Time to fill YGM edge maps (s): " << timer.elapsed() - << std::endl; } /** From 99ed2955eed664b8a76792d9e76e15fcd2bb5f9c Mon Sep 17 00:00:00 2001 From: Grace Li Date: Wed, 14 Oct 2026 13:39:33 -0700 Subject: [PATCH 3/5] Fix a typo --- src/clustering/distributed_clustering_mst_contraction.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/clustering/distributed_clustering_mst_contraction.hpp b/src/clustering/distributed_clustering_mst_contraction.hpp index eb5df40..feca716 100644 --- a/src/clustering/distributed_clustering_mst_contraction.hpp +++ b/src/clustering/distributed_clustering_mst_contraction.hpp @@ -93,7 +93,7 @@ void fill_init_min_incident_edge_map( } }; - // Fill the initial incidence array from the original MST edges + // Fill the initial incidence map from the original MST edges auto fill_init_min_incident_edge_map = [&min_incident_edge_map](const id_t &edge_id, std::pair &edge_endpoints) { From bad3170ff3a946acd17d78dc73ece5f8e14cb3a2 Mon Sep 17 00:00:00 2001 From: Grace Li Date: Wed, 14 Oct 2026 13:44:53 -0700 Subject: [PATCH 4/5] Remove extranous code from testing different ways to initialize min incident edge map --- ...distributed_clustering_mst_contraction.hpp | 70 ++----------------- .../run_distributed_hdbscan_clustering.cpp | 23 ------ 2 files changed, 4 insertions(+), 89 deletions(-) diff --git a/src/clustering/distributed_clustering_mst_contraction.hpp b/src/clustering/distributed_clustering_mst_contraction.hpp index feca716..27af64a 100644 --- a/src/clustering/distributed_clustering_mst_contraction.hpp +++ b/src/clustering/distributed_clustering_mst_contraction.hpp @@ -69,8 +69,7 @@ void sort_and_process_mst_edges_into_maps( } /** - * @brief Fill the initial min_incident_edge_map from the - edge_endpoints_array + * @brief Fill the initial min_incident_edge_map from the edge_endpoints_array * * @param edge_endpoints_array YGM array of edge id -> (edge_endpoint1, * edge_endpoint2) @@ -165,11 +164,9 @@ void contract_edges( edge_info.endpoint_supernode_reps.first, edge_info.endpoint_supernode_reps.second); }; - if (min_incident_edge_id != std::numeric_limits::max()) { - edge_contraction_array.async_visit(min_incident_edge_id, - add_edge_to_djset_lambda, - tree_components_djset_ptr); - } + edge_contraction_array.async_visit(min_incident_edge_id, + add_edge_to_djset_lambda, + tree_components_djset_ptr); // Set this supernode as the dendrogram child of the min incident edge if (round >= 2) { alpha_edge_map_ptr->async_visit( @@ -181,65 +178,6 @@ void contract_edges( comm.barrier(); } -/** - * @brief Find all MST edges to contract this round and add to disjoint set - * - * @param _round The round of MST contraction we're on. - * @param min_incident_edge_map YGM map of supernode representative -> - * min incident edge id. - * @param edge_contraction_array YGM array of edge id -> edge_contraction_info - * @param alpha_edge_map_ptr Pointer to YGM map of alpha edge id -> - * alpha_edge_info. Alpha edges are ones that are not contracted in the first - * round. As we contract alpha edges, we update their child supernodes in the - * dendrogram. - * @param tree_components_djset YGM disjoint_set for getting connected - * components after edge contraction. - */ -void contract_edges_first_round( - ygm::container::array &min_incident_edge_array, - ygm::container::array &edge_contraction_array, - ygm::ygm_ptr> alpha_edge_map_ptr, - ygm::container::disjoint_set &tree_components_djset) { - ygm::comm &comm = edge_contraction_array.comm(); - auto tree_components_djset_ptr = tree_components_djset.get_ygm_ptr(); - - [[maybe_unused]] static auto add_dendrogram_child_lambda = - []([[maybe_unused]] const id_t &edge_id, alpha_edge_info &info, - supernode_t new_child) { - if (info.dendrogram_children[0] == BLANK_SUPERNODE) { - info.dendrogram_children[0] = new_child; - } else { - info.dendrogram_children[1] = new_child; - } - }; - - auto contract_edges_lambda = [&edge_contraction_array, alpha_edge_map_ptr, - tree_components_djset_ptr]( - [[maybe_unused]] const id_t &supernode_rep, - const id_t &min_incident_edge_id) { - // The shortest edge incident to each supernode is contracted - - // Add the contracted edge to the disjoint set and set as - // contracted For an edge between supernodes (rep1, round-1) and - // (rep2, round-1), we add the rep1-rep2 edge to the disjoint set - auto add_edge_to_djset_lambda = []([[maybe_unused]] const id_t &edge_id, - edge_contraction_info &edge_info, - auto tree_components_djset_ptr) { - edge_info.contraction_round = 1; - tree_components_djset_ptr->async_union( - edge_info.endpoint_supernode_reps.first, - edge_info.endpoint_supernode_reps.second); - }; - if (min_incident_edge_id != std::numeric_limits::max()) { - edge_contraction_array.async_visit(min_incident_edge_id, - add_edge_to_djset_lambda, - tree_components_djset_ptr); - } - }; - min_incident_edge_array.for_all(contract_edges_lambda); - comm.barrier(); -} - /** * @brief Updates the edge contraction map by updating chain parent edge id of * contracted edges e by looking up the current supernode they are in and seeing diff --git a/src/clustering/run_distributed_hdbscan_clustering.cpp b/src/clustering/run_distributed_hdbscan_clustering.cpp index 4664903..5e13251 100644 --- a/src/clustering/run_distributed_hdbscan_clustering.cpp +++ b/src/clustering/run_distributed_hdbscan_clustering.cpp @@ -193,8 +193,6 @@ int main(int argc, char *argv[]) { mst_edge_vector; static id_t local_num_edges_read; local_num_edges_read = 0; - static id_t local_max_point_id; - local_max_point_id = 0; if (opt.metall_mst) { if (world.rank() == 0) { @@ -227,13 +225,6 @@ int main(int argc, char *argv[]) { std::make_pair((*input_mst_edges)[i].ids[0], (*input_mst_edges)[i].ids[1])); mst_edge_vector.push_back(array_item); - - if ((*input_mst_edges)[i].ids[0] > local_max_point_id) { - local_max_point_id = (*input_mst_edges)[i].ids[0]; - } - if ((*input_mst_edges)[i].ids[1] > local_max_point_id) { - local_max_point_id = (*input_mst_edges)[i].ids[1]; - } } local_num_edges_read = mst_edge_vector.size(); @@ -261,12 +252,6 @@ int main(int argc, char *argv[]) { distance_t dist; ss >> edge_id >> node1 >> node2 >> dist; ++local_num_edges_read; - if (node1 > local_max_point_id) { - local_max_point_id = node1; - } - if (node2 > local_max_point_id) { - local_max_point_id = node2; - } } catch (...) { std::cout << "Error reading mst line: " << line << std::endl; } @@ -300,12 +285,6 @@ int main(int argc, char *argv[]) { std::pair> array_item = std::make_pair(dist, std::make_pair(node1, node2)); mst_edge_vector.push_back(array_item); - if (node1 > local_max_point_id) { - local_max_point_id = node1; - } - if (node2 > local_max_point_id) { - local_max_point_id = node2; - } } catch (...) { std::cout << "Error reading mst line: " << line << std::endl; } @@ -328,7 +307,6 @@ int main(int argc, char *argv[]) { sw_step.reset(); id_t total_num_edges_read = ygm::sum(local_num_edges_read, world); - id_t max_point_id = ygm::max(local_max_point_id, world); // Initialize edge array of basic endpoints // and edge contraction array with edge_contraction_info @@ -391,7 +369,6 @@ int main(int argc, char *argv[]) { ss << " Number of MST edges: " << edge_endpoints_array.size(); if (world.rank() == 0) { spdlog::info(ss.str()); - spdlog::info(" Max point id: {}", max_point_id); spdlog::info(" Time to ingest MST edges (s): {:.3f}", sw_round); } From 41e2f3a0d0bbcb273dcf635aeee8c497e517fa73 Mon Sep 17 00:00:00 2001 From: Grace Li Date: Wed, 14 Oct 2026 13:57:31 -0700 Subject: [PATCH 5/5] Add timing for initializing min incident edge map, a current bottleneck for large number of comute nodes --- src/clustering/run_distributed_hdbscan_clustering.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/clustering/run_distributed_hdbscan_clustering.cpp b/src/clustering/run_distributed_hdbscan_clustering.cpp index 5e13251..3d660ac 100644 --- a/src/clustering/run_distributed_hdbscan_clustering.cpp +++ b/src/clustering/run_distributed_hdbscan_clustering.cpp @@ -412,6 +412,11 @@ int main(int argc, char *argv[]) { fill_init_min_incident_edge_map(min_incident_edge_map, edge_endpoints_array); incidence_map_time += sw_step.elapsed().count(); + if (opt.verbose && world.rank() == 0) { + spdlog::info( + " Time to initialize and fill min incident edge map (s): {:.3f}", + sw_step); + } // Max possible number of contraction rounds is log2(number of edges in // MST)