From 8b09aa49b064d81ec4ca6a3e3488acd8cf210bb6 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 7 Aug 2026 16:47:21 -0700 Subject: [PATCH 01/37] Change default AMST approx bound value to 1.1 --- script/benchmark/run_clams_bench.py | 18 +++++++--- src/mst/build_amst.cpp | 56 ++++++++++++++--------------- 2 files changed, 42 insertions(+), 32 deletions(-) diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index 90fba81..c810d4d 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -59,6 +59,11 @@ def parse_options(): default='', help='Path to an already constructed DNND PM datastore. If specified, skip the DNND step.') + # Approximate Minimum Spanning Tree (AMST) + # -a: approximate bound, default 1.1 + parser.add_argument('-a', '--amst_approx_bound', type=float, default=1.1, + help='Approximation bound for AMST. Default is 1.1.') + # For HDBSCAN # Min cluster size, conmma separated list of min cluster sizes, or range of min cluster sizes parser.add_argument('-m', '--min_cluster_size', default='5', @@ -148,7 +153,8 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, dnnd_exe, nng_k, distance_func, points_file_format, point_path, mfc_exe, - amst_exe, clustering_exe, + amst_exe, amst_approx_bound, + clustering_exe, evaluator, ygm_cluster_eval, verbose, ground_truth_path, @@ -202,7 +208,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, job_script.write("date\n") job_script.write(f"echo \"Running AMST\"\n") amst_ds_path = f"{work_dir}/amst_pm_datastore" - amst_command = f"{amst_exe} -d {dnnd_ds_path} -p {amst_ds_path}" + amst_command = f"{amst_exe} -d {dnnd_ds_path} -p {amst_ds_path} -e {amst_approx_bound}" add_srun_cmd(num_tasks_per_node, amst_command, job_script) # Run the HPC Clustering step @@ -297,6 +303,7 @@ def main(): opts.point_path, opts.mfc_exe, opts.amst_exe, + opts.amst_approx_bound, opts.clustering_exe, evaluator, opts.ygm_cluster_eval, @@ -307,17 +314,20 @@ def main(): opts.input_dnnd_ds_path) print(f"Generated batch script: {job_script}") + job_submission_cmd = f"sbatch {opts.sbatch_opts} {job_script}" # Write job execution commands log with open(f'{job_dir}/info.txt', 'w') as f: f.write(f"Command executed by user:\n") f.write(' '.join(sys.argv)) f.write('\n\n') - f.write(f"To submit the job, run:\nsbatch {opts.sbatch_opts} {job_script}\n") + f.write(f"To submit the job, run:\n{job_submission_cmd}\n") # Submit the job if opts.submit_job: print(f"Submit job: {job_script}") - os.system(f"sbatch {opts.sbatch_opts} {job_script}") + os.system(job_submission_cmd) + else: + print(f"To submit the job, run:\n{job_submission_cmd}\n") if __name__ == '__main__': diff --git a/src/mst/build_amst.cpp b/src/mst/build_amst.cpp index 95f43b0..73f6d02 100644 --- a/src/mst/build_amst.cpp +++ b/src/mst/build_amst.cpp @@ -15,15 +15,15 @@ #include #include -#include #include +#include -#include "approx_mst_builder.hpp" #include "../common.hpp" +#include "approx_mst_builder.hpp" namespace cls = clams; -using id_t = cls::id_t; +using id_t = cls::id_t; using distance_t = cls::distance_t; using edge_t = std::tuple; @@ -45,8 +45,8 @@ void usage(ygm::comm &comm) { void parse_cmd_line(int argc, char **argv, ygm::comm &comm, float &approx_bound, std::vector &txt_input_filenames, - std::string &pm_input_path, - std::string &pm_output_filename) { + std::string &pm_input_path, + std::string &pm_output_filename) { if (comm.rank0()) { std::cout << "CMD line:"; for (int i = 0; i < argc; ++i) { @@ -55,9 +55,9 @@ void parse_cmd_line(int argc, char **argv, ygm::comm &comm, float &approx_bound, std::cout << std::endl; } - int c; + int c; bool inserting_input_filenames = false; - bool prn_help = false; + bool prn_help = false; while (true) { while ((c = getopt(argc, argv, "+e:i:d:p:h ")) != -1) { inserting_input_filenames = false; @@ -113,7 +113,7 @@ void parse_cmd_line(int argc, char **argv, ygm::comm &comm, float &approx_bound, } std::vector read_dnnd_output(const std::vector &filenames, - ygm::comm &c) { + ygm::comm &c) { std::vector to_return; ygm::container::bag filenames_bag(c); @@ -138,8 +138,8 @@ std::vector read_dnnd_output(const std::vector &filenames, std::stringstream neighbors_ss(neighbors_str); std::stringstream distances_ss(distances_str); - id_t src; - id_t ngbr; + id_t src; + id_t ngbr; distance_t dist; neighbors_ss >> src; @@ -159,7 +159,7 @@ std::vector read_dnnd_output(const std::vector &filenames, } std::vector read_pm_knng_output(std::string &dstore_path, - ygm::comm &comm) { + ygm::comm &comm) { comm.cout0() << "Reading distributed PM kNNG data from: " << dstore_path << std::endl; @@ -168,7 +168,7 @@ std::vector read_pm_knng_output(std::string &dstore_path, const auto &knng = pm_knng.get_knng(); std::vector to_return; - std::size_t local_num_neighbors = 0; + std::size_t local_num_neighbors = 0; for (const auto &[source, neighbors] : knng) { (void)source; local_num_neighbors += neighbors.size(); @@ -229,7 +229,7 @@ void pivot_sort(std::vector &in_vec, ygm::comm &world) { return a.second < b.second; }; - const size_t samples_per_pivot = 40; + const size_t samples_per_pivot = 40; std::vector to_sort; to_sort.reserve(in_vec.size() * 1.1f); @@ -238,14 +238,14 @@ void pivot_sort(std::vector &in_vec, ygm::comm &world) { // of duplicate items std::vector> samples; std::vector> pivots; - static auto &s_samples = samples; - static auto &s_to_sort = to_sort; + static auto &s_samples = samples; + static auto &s_to_sort = to_sort; samples.reserve(world.size() * samples_per_pivot); // std::default_random_engine rng; - size_t my_prefix = ygm::prefix_sum(in_vec.size(), world); + size_t my_prefix = ygm::prefix_sum(in_vec.size(), world); size_t global_size = ygm::sum(in_vec.size(), world); std::uniform_int_distribution uintdist{0, global_size - 1}; @@ -275,9 +275,9 @@ void pivot_sort(std::vector &in_vec, ygm::comm &world) { // // Partition using pivots for (size_t i = 0; i < in_vec.size(); ++i) { - auto itr = std::lower_bound(pivots.begin(), pivots.end(), - std::make_pair(in_vec[i], my_prefix + i), - edge_index_pair_comp_lambda); + auto itr = std::lower_bound(pivots.begin(), pivots.end(), + std::make_pair(in_vec[i], my_prefix + i), + edge_index_pair_comp_lambda); size_t owner = std::distance(pivots.begin(), itr); world.async( @@ -319,7 +319,7 @@ void pivot_sort(std::vector &in_vec, ygm::comm &world) { in_vec.swap(to_sort); } -void write_output_single_file(const std::string &output_filename, +void write_output_single_file(const std::string &output_filename, const std::vector &edges, ygm::comm &c) { MPI_Comm mpi_comm = c.get_mpi_comm(); @@ -348,14 +348,14 @@ void write_output_single_file(const std::string &output_filename, c.cout0("Dumped AMST edges to: ", output_filename); } -void write_output_single_pm(const std::string &output_filename, +void write_output_single_pm(const std::string &output_filename, const std::vector &edges, ygm::comm &c) { const size_t num_edges = c.all_reduce_sum(edges.size()); MPI_Comm mpi_comm = c.get_mpi_comm(); if (c.rank0()) { metall::manager manager(metall::create_only, output_filename); - auto *amst_edges = manager.construct( + auto *amst_edges = manager.construct( metall::unique_instance)(manager.get_allocator()); amst_edges->reserve(num_edges); @@ -385,18 +385,18 @@ int main(int argc, char **argv) { ygm::comm world(&argc, &argv); std::vector txt_input_filenames; - std::string pm_input_path; - std::string pm_output_filename; - float amst_approx_bound{2.0}; + std::string pm_input_path; + std::string pm_output_filename; + float amst_approx_bound{1.1}; parse_cmd_line(argc, argv, world, amst_approx_bound, txt_input_filenames, pm_input_path, pm_output_filename); world.cout0("Reading edges"); ygm::utility::timer read_timer; - auto edges = (txt_input_filenames.size() > 0) - ? read_dnnd_output(txt_input_filenames, world) - : read_pm_knng_output(pm_input_path, world); + auto edges = (txt_input_filenames.size() > 0) + ? read_dnnd_output(txt_input_filenames, world) + : read_pm_knng_output(pm_input_path, world); world.cout0("Edge reading time (s): ", read_timer.elapsed()); auto amst_edges = approx_mst(edges, amst_approx_bound, world); From a8e620494c3781ca70a8629ce37708c3f9dcf186 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 7 Aug 2026 16:47:36 -0700 Subject: [PATCH 02/37] Use MFC by default --- script/benchmark/run_clams_bench.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index c810d4d..498a3c4 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -114,7 +114,7 @@ def parse_options(): default=f'{cwd}/src/knng/build_knng', help='Path to the DNND executable.') parser.add_argument('-M', '--mfc_exe', - default=f'{cwd}/src/mfc/connect_ccs_random_dist', + default=f'{cwd}/src/mfc/connect_mfc', help='Path to the MFC executable.') parser.add_argument('-A', '--amst_exe', default=f'{cwd}/src/mst/build_amst', From 805652b86cd22aa927441621c66e7cd8baac6bf4 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 7 Aug 2026 17:15:17 -0700 Subject: [PATCH 03/37] Brush up scripts --- script/benchmark/clustering_utilities.py | 2 +- script/benchmark/run_clams_bench.py | 23 ++++++++++++++++++----- 2 files changed, 19 insertions(+), 6 deletions(-) diff --git a/script/benchmark/clustering_utilities.py b/script/benchmark/clustering_utilities.py index abb873b..d2cf2de 100644 --- a/script/benchmark/clustering_utilities.py +++ b/script/benchmark/clustering_utilities.py @@ -27,7 +27,7 @@ def eval_clusters(cluster_labels, true_labels): cluster_labels[clustered_points]) # sil = silhouette_score(raw_data[clustered_points], cluster_labels[clustered_points]) pct_clustered = (np.sum(clustered_points) / cluster_labels.shape[0]) - print(f"Clustered Points: {pct_clustered * 100:.2f}%") + print(f"Cluster coverage (%): {pct_clustered * 100:.2f}") else: ari = adjusted_rand_score(true_labels, cluster_labels) ami = adjusted_mutual_info_score(true_labels, cluster_labels) diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index 498a3c4..da950c2 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -173,9 +173,12 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, add_cmd(f'mkdir -p {work_dir}', job_script, True, False) # Run the DNND step + job_script.write("echo\n") job_script.write("date\n") if len(input_dnnd_ds_path) == 0: + job_script.write(f"echo ================================\n") job_script.write(f"echo \"Building KNNG\"\n") + job_script.write(f"echo ================================\n") dnnd_ds_path = f"{work_dir}/dnnd_pm_datastore" dnnd_batch_size = 2 ** 25 verbose_flag = '-v' if verbose else '' @@ -187,33 +190,43 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, dnnd_ds_path = input_dnnd_ds_path # Connect the CCs + job_script.write("echo\n") job_script.write("date\n") - job_script.write( - f"echo \"Randomly connect components\"\n") + job_script.write(f"echo ================================\n") + job_script.write(f"echo \"Connect components\"\n") + job_script.write(f"echo ================================\n") mfc_command = f"{mfc_exe} -d {dnnd_ds_path} -f {distance_func}" add_srun_cmd(num_tasks_per_node, mfc_command, job_script) # Convert to core distance # TODO: Implement if False and min_samples > 0: + job_script.write("echo\n") job_script.write("date\n") - job_script.write( - f"echo \"Convert to core distance kNNG\"\n") + job_script.write(f"echo ================================\n") + job_script.write(f"echo Convert to core distance kNNG\n") + job_script.write(f"echo ================================\n") knng_coredist_dir = f"{work_dir}/knng_coredist/" add_cmd(f'mkdir -p {knng_coredist_dir}', job_script) conv2coredist_cmd = f"./src/conv_knng_to_core_dist -i {dnnd_ds_path} -o {knng_coredist_dir}/knng.txt -m {min_samples}" add_cmd(conv2coredist_cmd, job_script) # Run the AMST step + job_script.write("echo\n") job_script.write("date\n") + job_script.write(f"echo ================================\n") job_script.write(f"echo \"Running AMST\"\n") + job_script.write(f"echo ================================\n") amst_ds_path = f"{work_dir}/amst_pm_datastore" amst_command = f"{amst_exe} -d {dnnd_ds_path} -p {amst_ds_path} -e {amst_approx_bound}" add_srun_cmd(num_tasks_per_node, amst_command, job_script) # Run the HPC Clustering step + job_script.write("echo\n") job_script.write("date\n") - job_script.write(f"echo \"Running HPC Clustering\"\n") + job_script.write(f"echo ================================\n") + job_script.write(f"echo \"Running CLAMS-HDBSCAN\"\n") + job_script.write(f"echo ================================\n") for try_no, set_cmd in enumerate(min_cluster_size_set_cmnds): add_cmd(set_cmd, job_script, False, False) From 6932cc6dbfe3178f65b204bf7e7e97eca8a2fcf4 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 7 Aug 2026 18:40:47 -0700 Subject: [PATCH 04/37] Update bench scripts --- script/benchmark/bench_utilities.py | 13 +++- script/benchmark/run_clams_bench.py | 112 +++++++++++++++------------- 2 files changed, 72 insertions(+), 53 deletions(-) diff --git a/script/benchmark/bench_utilities.py b/script/benchmark/bench_utilities.py index 47ee707..4d16045 100644 --- a/script/benchmark/bench_utilities.py +++ b/script/benchmark/bench_utilities.py @@ -84,6 +84,9 @@ def grep_file(file, text_pattern): lines = f.readlines() return [line for line in lines if text_pattern in line] +# Parse a range string and return a list of integers +# e.g., "1-5" -> [1, 2, 3, 4, 5] +# e.g., "1,3,5" -> [1, 3, 5] def parse_range(range_str): if range_str is None: return [None] @@ -92,4 +95,12 @@ def parse_range(range_str): start, end = range_str.split('-') return list(range(int(start), int(end) + 1)) else: - return [int(x) for x in range_str.split(',')] \ No newline at end of file + return [int(x) for x in range_str.split(',')] + +# Parse a list of floats from a string and return a list of floats +# e.g., "1.0,2.5,3.14" -> [1.0, 2.5, 3.14] +def parse_float_list(float_list_str): + if float_list_str is None: + return [None] + + return [float(x) for x in float_list_str.split(',')] \ No newline at end of file diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index da950c2..7fd3770 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -61,15 +61,16 @@ def parse_options(): # Approximate Minimum Spanning Tree (AMST) # -a: approximate bound, default 1.1 - parser.add_argument('-a', '--amst_approx_bound', type=float, default=1.1, - help='Approximation bound for AMST. Default is 1.1.') + parser.add_argument('-a', '--amst_approx_bound', default='1.1', + dest='amst_approx_bound_range', action='store', type=str, + help='Approximation bound for AMST. Single value or comma separated list of values (e.g., 1.1,1.2,1.4). Default is 1.1.') # For HDBSCAN # Min cluster size, conmma separated list of min cluster sizes, or range of min cluster sizes parser.add_argument('-m', '--min_cluster_size', default='5', dest='min_cluster_size_range', action='store', type=str, - help='Minimum cluster size. Comma separated list of values or range of values (e.g., 5,10,15 or 5-15)') + help='Minimum cluster size. Single value or, comma separated list of values or range of values (e.g., 5,10,15 or 5-15)') # Min samples parser.add_argument('-s', '--min_samples', type=int, default=-1, dest='min_samples', @@ -153,7 +154,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, dnnd_exe, nng_k, distance_func, points_file_format, point_path, mfc_exe, - amst_exe, amst_approx_bound, + amst_exe, amst_approx_bound_list, clustering_exe, evaluator, ygm_cluster_eval, verbose, @@ -193,7 +194,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, job_script.write("echo\n") job_script.write("date\n") job_script.write(f"echo ================================\n") - job_script.write(f"echo \"Connect components\"\n") + job_script.write(f"echo \"Running MFS\"\n") job_script.write(f"echo ================================\n") mfc_command = f"{mfc_exe} -d {dnnd_ds_path} -f {distance_func}" add_srun_cmd(num_tasks_per_node, mfc_command, job_script) @@ -211,53 +212,58 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, conv2coredist_cmd = f"./src/conv_knng_to_core_dist -i {dnnd_ds_path} -o {knng_coredist_dir}/knng.txt -m {min_samples}" add_cmd(conv2coredist_cmd, job_script) - # Run the AMST step - job_script.write("echo\n") - job_script.write("date\n") - job_script.write(f"echo ================================\n") - job_script.write(f"echo \"Running AMST\"\n") - job_script.write(f"echo ================================\n") - amst_ds_path = f"{work_dir}/amst_pm_datastore" - amst_command = f"{amst_exe} -d {dnnd_ds_path} -p {amst_ds_path} -e {amst_approx_bound}" - add_srun_cmd(num_tasks_per_node, amst_command, job_script) - - # Run the HPC Clustering step - job_script.write("echo\n") - job_script.write("date\n") - job_script.write(f"echo ================================\n") - job_script.write(f"echo \"Running CLAMS-HDBSCAN\"\n") - job_script.write(f"echo ================================\n") - for try_no, set_cmd in enumerate(min_cluster_size_set_cmnds): - add_cmd(set_cmd, job_script, False, False) - - job_script.write( - f"echo \"Min clustering size ${{MIN_CLUSTER_SIZE}}\"\n") - cluster_label_file = f"{work_dir}/cluster_labels-try{try_no}.txt" - cluster_tree_file = f"{work_dir}/cluster_tree-try{try_no}.txt" - hpc_clustering_command = (f"{clustering_exe} -i {amst_ds_path} -M " - f" -m ${{MIN_CLUSTER_SIZE}} " - f" -o {cluster_label_file} " - f" -c {cluster_tree_file} " - f" -P ") - add_cmd(hpc_clustering_command, job_script) - - # Run the evaluation step - if ground_truth_path: - if ygm_cluster_eval: - job_script.write( - f"echo \"Evaluating Clustering using YGM\"\n") - if verbose: - evaluation_command = f"{evaluator} -v -g {ground_truth_path} {cluster_label_file}" + try_no = 0 + for amst_approx_bound in amst_approx_bound_list: + # Run the AMST step + job_script.write("echo\n") + job_script.write("date\n") + job_script.write(f"echo ================================\n") + job_script.write(f"echo \"Running AMST, approx bound = {amst_approx_bound}\"\n") + job_script.write(f"echo ================================\n") + amst_ds_path = f"{work_dir}/amst_pm_datastore_a{amst_approx_bound}" + amst_command = f"{amst_exe} -d {dnnd_ds_path} -p {amst_ds_path} -e {amst_approx_bound}" + add_srun_cmd(num_tasks_per_node, amst_command, job_script) + + # Run the HPC Clustering step + for set_cmd in min_cluster_size_set_cmnds: + # Set the min cluster size environment variable for this run + add_cmd(set_cmd, job_script, False, False) + + job_script.write("echo\n") + job_script.write("date\n") + job_script.write(f"echo ================================\n") + job_script.write(f"echo \"Running CLAMS-HDBSCAN, min cluster size = ${{MIN_CLUSTER_SIZE}}\"\n") + job_script.write(f"echo ================================\n") + + job_script.write( + f"echo \"Min cluster size ${{MIN_CLUSTER_SIZE}}\"\n") + cluster_label_file = f"{work_dir}/cluster_labels_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}.txt" + cluster_tree_file = f"{work_dir}/cluster_tree_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}.txt" + hpc_clustering_command = (f"{clustering_exe} -i {amst_ds_path} -M " + f" -m ${{MIN_CLUSTER_SIZE}} " + f" -o {cluster_label_file} " + f" -c {cluster_tree_file} " + f" -P ") + add_cmd(hpc_clustering_command, job_script) + + # Run the evaluation step + if ground_truth_path: + if ygm_cluster_eval: + job_script.write( + f"echo \"Evaluating Clustering using YGM\"\n") + if verbose: + evaluation_command = f"{evaluator} -v -g {ground_truth_path} {cluster_label_file}" + else: + evaluation_command = f"{evaluator} -g {ground_truth_path} {cluster_label_file}" + add_srun_cmd(num_tasks_per_node, evaluation_command, + job_script) else: - evaluation_command = f"{evaluator} -g {ground_truth_path} {cluster_label_file}" - add_srun_cmd(num_tasks_per_node, evaluation_command, - job_script) - else: - job_script.write( - f"echo \"Evaluating Clustering using python script\"\n") - evaluation_command = f"python3 {evaluator} -c {cluster_label_file} -g {ground_truth_path}" - add_cmd(evaluation_command, job_script) - job_script.write(f"echo \"\" \n") + job_script.write( + f"echo \"Evaluating Clustering using python script\"\n") + evaluation_command = f"python3 {evaluator} -c {cluster_label_file} -g {ground_truth_path}" + add_cmd(evaluation_command, job_script) + job_script.write(f"echo \"\" \n") + try_no += 1 # If the file was not created, return an error if not os.path.exists(job_script_path): @@ -306,6 +312,8 @@ def main(): min_cluster_size_set_cmnds = [f'MIN_CLUSTER_SIZE={x}' for x in min_cluster_size_list] + amst_approx_bound_list = parse_float_list(opts.amst_approx_bound_range) + # Generate a benchmark batch script job_script = gen_clams_bench_script(job_name, job_dir, work_dir, opts.num_nodes, @@ -316,7 +324,7 @@ def main(): opts.point_path, opts.mfc_exe, opts.amst_exe, - opts.amst_approx_bound, + amst_approx_bound_list, opts.clustering_exe, evaluator, opts.ygm_cluster_eval, From 3ad72d64c94788307641a8cc0de01bc9ab92502a Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Sat, 8 Aug 2026 02:16:12 -0700 Subject: [PATCH 05/37] Update bench scripts --- script/benchmark/run_clams_bench.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index 7fd3770..42a9348 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -54,16 +54,21 @@ def parse_options(): 'a KNNG value.') parser.add_argument('--neodnnd', action='store_true', help='Use NEO-DNND (build_knng_neo) instead of build_knng.') + parser.add_argument('--neodnnd_threads', type=int, default=2, + help='Number of threads to use for NEO-DNND. Ignored if --neodnnd is not specified.') parser.add_argument('-G', '--input_dnnd_ds_path', default='', help='Path to an already constructed DNND PM datastore. If specified, skip the DNND step.') + # KNNG backup + parser.add_argument('-R', '--backup_knng', action='store_true', + help='Backup the KNNG datastore after DNND step.') + # Approximate Minimum Spanning Tree (AMST) - # -a: approximate bound, default 1.1 parser.add_argument('-a', '--amst_approx_bound', default='1.1', dest='amst_approx_bound_range', action='store', type=str, - help='Approximation bound for AMST. Single value or comma separated list of values (e.g., 1.1,1.2,1.4). Default is 1.1.') + help='Approximation bound for AMST. Single value or comma separated list of values (e.g., 1.1,1.2,1.4).') # For HDBSCAN # Min cluster size, conmma separated list of min cluster sizes, or range of min cluster sizes @@ -153,6 +158,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, num_nodes, num_tasks_per_node, dnnd_exe, nng_k, distance_func, points_file_format, point_path, + backup_knng, mfc_exe, amst_exe, amst_approx_bound_list, clustering_exe, @@ -185,6 +191,13 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, verbose_flag = '-v' if verbose else '' dnnd_command = f"{dnnd_exe} {verbose_flag} -k {nng_k} -f {distance_func} -o {dnnd_ds_path} -b {dnnd_batch_size} -p {points_file_format} {point_path}" add_srun_cmd(num_tasks_per_node, dnnd_command, job_script) + if backup_knng: + dnnd_ds_path_backup = f"{dnnd_ds_path}_backup" + job_script.write(f"echo\n") + job_script.write("date\n") + job_script.write(f"echo \"Backing up KNNG datastore\"\n") + backup_knng_command = f"cp -r {dnnd_ds_path} {dnnd_ds_path_backup}" + add_cmd(backup_knng_command, job_script) else: job_script.write( f"Using existing DNND datastore at {input_dnnd_ds_path}\n") @@ -302,6 +315,8 @@ def main(): default_dnnd_exe = f'{os.getcwd()}/src/knng/build_knng' if dnnd_exe == default_dnnd_exe: dnnd_exe = f'{os.getcwd()}/src/knng/build_knng_neo' + # This is not the best way to set the number of threads for NEO-DNND, but it is a simple way to do it for now. + dnnd_exe = f'{dnnd_exe} -T {opts.neodnnd_threads}' if opts.ygm_cluster_eval: evaluator = opts.ygm_evaluator_exe @@ -322,6 +337,7 @@ def main(): opts.distance_func, opts.points_file_format, opts.point_path, + opts.backup_knng, opts.mfc_exe, opts.amst_exe, amst_approx_bound_list, From e10e513ff4afd657054ac134053c2f46d3db9838 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Mon, 10 Aug 2026 12:30:38 -0700 Subject: [PATCH 06/37] Minor code brush ups --- src/common.hpp | 24 ++++++++++++------------ src/mst/build_amst.cpp | 2 +- src/mst/run_kruskal_mst.cpp | 20 +++++++++++++------- 3 files changed, 26 insertions(+), 20 deletions(-) diff --git a/src/common.hpp b/src/common.hpp index 7e01cd4..b75efc6 100644 --- a/src/common.hpp +++ b/src/common.hpp @@ -1,7 +1,6 @@ // Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS // Project Developers. See the top-level COPYRIGHT file for details. - #pragma once #include @@ -55,7 +54,7 @@ inline std::vector find_files( /// \brief Read edges. /// \param path A path to an edge file or a directory that contains edge files. inline void read_edges(const std::filesystem::path &path, - weighted_edge_list_t &edges) { + weighted_edge_list_t &edges) { const auto files = find_files(path); for (const auto &file : files) { std::ifstream ifs(file); @@ -66,8 +65,8 @@ inline void read_edges(const std::filesystem::path &path, std::string line; while (true) { - id_t src; - id_t dst; + id_t src; + id_t dst; distance_t dist; ifs >> src >> dst >> dist; if (ifs.eof()) { @@ -78,12 +77,13 @@ inline void read_edges(const std::filesystem::path &path, } } -inline void read_knn_edges(const std::vector &knng_files, - weighted_edge_list_t &edges) { +inline void read_knng_edges( + const std::vector &knng_files, + weighted_edge_list_t &edges) { std::size_t num_edges = 0; OMP_DIRECTIVE(parallel for reduction(+ : num_edges)) for (std::size_t fno = 0; fno < knng_files.size(); ++fno) { - const auto &file = knng_files[fno]; + const auto &file = knng_files[fno]; std::ifstream ifs(file); if (!ifs.is_open()) { std::cerr << "Cannot open file: " << file << std::endl; @@ -97,7 +97,7 @@ inline void read_knn_edges(const std::vector &knng_files, break; // End of file } std::istringstream iss(line); - id_t buf; + id_t buf; while (iss >> buf) { ++num_edges; } @@ -112,7 +112,7 @@ inline void read_knn_edges(const std::vector &knng_files, std::atomic cnt_edges{0}; OMP_DIRECTIVE(parallel for) for (std::size_t fno = 0; fno < knng_files.size(); ++fno) { - const auto &file = knng_files[fno]; + const auto &file = knng_files[fno]; std::ifstream ifs(file); if (!ifs.is_open()) { std::cerr << "Cannot open file: " << file << std::endl; @@ -125,7 +125,7 @@ inline void read_knn_edges(const std::vector &knng_files, { std::getline(ifs, line); std::istringstream iss(line); - id_t buf; + id_t buf; while (iss >> buf) { ids.push_back(buf); } @@ -135,7 +135,7 @@ inline void read_knn_edges(const std::vector &knng_files, { std::getline(ifs, line); std::istringstream iss(line); - distance_t buf; + distance_t buf; while (iss >> buf) { dists.push_back(buf); } @@ -154,7 +154,7 @@ inline void read_knn_edges(const std::vector &knng_files, const id_t src = ids[0]; for (std::size_t i = 1; i < ids.size(); ++i) { const auto index = cnt_edges.fetch_add(1); - edges[index] = weighted_edge_t{src, ids[i], dists[i]}; + edges[index] = weighted_edge_t{src, ids[i], dists[i]}; } } } diff --git a/src/mst/build_amst.cpp b/src/mst/build_amst.cpp index 73f6d02..be40d97 100644 --- a/src/mst/build_amst.cpp +++ b/src/mst/build_amst.cpp @@ -143,7 +143,7 @@ std::vector read_dnnd_output(const std::vector &filenames, distance_t dist; neighbors_ss >> src; - distances_ss >> dist; + distances_ss >> dist; // Dummy value, not used while (neighbors_ss >> ngbr) { distances_ss >> dist; diff --git a/src/mst/run_kruskal_mst.cpp b/src/mst/run_kruskal_mst.cpp index cebf442..0a7a5ca 100644 --- a/src/mst/run_kruskal_mst.cpp +++ b/src/mst/run_kruskal_mst.cpp @@ -1,7 +1,6 @@ // Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS // Project Developers. See the top-level COPYRIGHT file for details. - // Find the MST using Kruskal's algorithm // If GCC is used, @@ -11,11 +10,11 @@ #warning "No parallelism library is used." #endif +#include #include #include #include #include -#include #include #include @@ -27,14 +26,13 @@ namespace omp = metall::utility::omp; bool parse_option(int argc, char *argv[], std::filesystem::path &knng_dir, std::filesystem::path &output_dir) { - if (argc != 3) { std::cerr << "Usage: " << argv[0] << " " << std::endl; return false; } - knng_dir = std::filesystem::path(argv[1]); + knng_dir = std::filesystem::path(argv[1]); output_dir = std::filesystem::path(argv[2]); return true; @@ -60,7 +58,7 @@ int main(int argc, char *argv[]) { clams::weighted_edge_list_t edges; spdlog::info("Read KNNG edges"); - clams::read_knn_edges(knng_files, edges); + clams::read_knng_edges(knng_files, edges); spdlog::info("#of edges {}", edges.size()); spdlog::info("Sort edges"); @@ -76,7 +74,7 @@ int main(int argc, char *argv[]) { spdlog::info("Finished sorting edges"); std::unordered_map parent_tbl; - clams::weighted_edge_list_t mst_edges; + clams::weighted_edge_list_t mst_edges; spdlog::info("Start Kruskal's algorithm"); for (const auto &edge : edges) { const auto src = edge.ids[0]; @@ -103,9 +101,17 @@ int main(int argc, char *argv[]) { } spdlog::info("Finished Kruskal's algorithm"); - std::filesystem::create_directories(output_dir); + // Show the sum of distances of the MST edges + { + double total_distance = 0.0; + for (const auto &edge : mst_edges) { + total_distance += edge.distance; + } + spdlog::info("Total distance of MST edges: {}", total_distance); + } spdlog::info("Write MST edges int {}", output_dir.string()); + std::filesystem::create_directories(output_dir); OMP_DIRECTIVE(parallel) { std::string path = output_dir / From 5c17de2efefccb0dfc55caa90af223732f5efa02 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Mon, 10 Aug 2026 17:07:54 -0700 Subject: [PATCH 07/37] Update bench script --- script/benchmark/run_clams_bench.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index 42a9348..f6953a9 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -56,6 +56,10 @@ def parse_options(): help='Use NEO-DNND (build_knng_neo) instead of build_knng.') parser.add_argument('--neodnnd_threads', type=int, default=2, help='Number of threads to use for NEO-DNND. Ignored if --neodnnd is not specified.') + parser.add_argument('--nng_r', type=float, default=0.5, + help='The r (sampling) parameter for KNNG construction.') + parser.add_argument('--nng_delta', type=float, default=0.0001, + help='The delta parameter (termination condition) for KNNG construction.') parser.add_argument('-G', '--input_dnnd_ds_path', default='', @@ -158,6 +162,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, num_nodes, num_tasks_per_node, dnnd_exe, nng_k, distance_func, points_file_format, point_path, + nng_r, nng_delta, backup_knng, mfc_exe, amst_exe, amst_approx_bound_list, @@ -189,7 +194,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, dnnd_ds_path = f"{work_dir}/dnnd_pm_datastore" dnnd_batch_size = 2 ** 25 verbose_flag = '-v' if verbose else '' - dnnd_command = f"{dnnd_exe} {verbose_flag} -k {nng_k} -f {distance_func} -o {dnnd_ds_path} -b {dnnd_batch_size} -p {points_file_format} {point_path}" + dnnd_command = f"{dnnd_exe} {verbose_flag} -k {nng_k} -r {nng_r} -d {nng_delta} -f {distance_func} -o {dnnd_ds_path} -b {dnnd_batch_size} -p {points_file_format} {point_path}" add_srun_cmd(num_tasks_per_node, dnnd_command, job_script) if backup_knng: dnnd_ds_path_backup = f"{dnnd_ds_path}_backup" @@ -337,6 +342,8 @@ def main(): opts.distance_func, opts.points_file_format, opts.point_path, + opts.nng_r, + opts.nng_delta, opts.backup_knng, opts.mfc_exe, opts.amst_exe, From 8189ff9d9d66f597e0db9cf9d23ad2f52f5aa286 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Thu, 13 Aug 2026 07:02:40 -0700 Subject: [PATCH 08/37] Assign clusters to noise points --- script/benchmark/clustering_utilities.py | 30 ++- .../benchmark/evaluate_clustering_quality.py | 8 +- script/benchmark/run_clams_bench.py | 86 +++++-- src/clustering/CMakeLists.txt | 5 + src/clustering/cluster_noise_points.cpp | 218 ++++++++++++++++++ 5 files changed, 320 insertions(+), 27 deletions(-) create mode 100644 src/clustering/cluster_noise_points.cpp diff --git a/script/benchmark/clustering_utilities.py b/script/benchmark/clustering_utilities.py index d2cf2de..e72f8be 100644 --- a/script/benchmark/clustering_utilities.py +++ b/script/benchmark/clustering_utilities.py @@ -10,24 +10,35 @@ # Modified from: # https://gist.github.com/lmcinnes/24ed5c22c80125be5133811d677eae7b -def eval_clusters(cluster_labels, true_labels): +def eval_clusters(cluster_labels, true_labels, singleton_cluster_to_noise_points=False): + max_cluster_id = max(np.max(cluster_labels), np.max(true_labels)) + if np.any(true_labels < 0): print("Ground truth labels contain noise points") pct_clustered_gt = (np.sum(true_labels >= 0) / cluster_labels.shape[0]) print(f"GT clustered Points: {pct_clustered_gt * 100:.2f}%") print( "Assigning a singleton cluster to each noise point in the ground truth labels") - true_labels = assign_singleton_cluster_to_noise_point(true_labels) + true_labels = assign_singleton_cluster_to_noise_points(true_labels, max_cluster_id) + max_cluster_id = max(max_cluster_id, np.max(true_labels)) if np.any(cluster_labels < 0): # Has noise points clustered_points = (cluster_labels >= 0) - ari = adjusted_rand_score(true_labels[clustered_points], - cluster_labels[clustered_points]) - ami = adjusted_mutual_info_score(true_labels[clustered_points], - cluster_labels[clustered_points]) - # sil = silhouette_score(raw_data[clustered_points], cluster_labels[clustered_points]) pct_clustered = (np.sum(clustered_points) / cluster_labels.shape[0]) print(f"Cluster coverage (%): {pct_clustered * 100:.2f}") + + if singleton_cluster_to_noise_points: + print( + "Assigning a singleton cluster to each noise point in the clustering result") + cluster_labels = assign_singleton_cluster_to_noise_points(cluster_labels, max_cluster_id) + ari = adjusted_rand_score(true_labels, cluster_labels) + ami = adjusted_mutual_info_score(true_labels, cluster_labels) + else: + ari = adjusted_rand_score(true_labels[clustered_points], + cluster_labels[clustered_points]) + ami = adjusted_mutual_info_score(true_labels[clustered_points], + cluster_labels[clustered_points]) + # sil = silhouette_score(raw_data[clustered_points], cluster_labels[clustered_points]) else: ari = adjusted_rand_score(true_labels, cluster_labels) ami = adjusted_mutual_info_score(true_labels, cluster_labels) @@ -39,15 +50,14 @@ def eval_clusters(cluster_labels, true_labels): # Assign a cluster ID to each noise point -def assign_singleton_cluster_to_noise_point(cluster_labels): +def assign_singleton_cluster_to_noise_points(cluster_labels, noise_id_offset): new_labels = cluster_labels.copy() - max_cluster_id = np.max(cluster_labels) cnt_noise = 0 for i, label in enumerate(cluster_labels): if label == -1: cnt_noise += 1 - new_labels[i] = max_cluster_id + cnt_noise + new_labels[i] = cnt_noise + noise_id_offset return new_labels diff --git a/script/benchmark/evaluate_clustering_quality.py b/script/benchmark/evaluate_clustering_quality.py index 81bef75..0cc5813 100644 --- a/script/benchmark/evaluate_clustering_quality.py +++ b/script/benchmark/evaluate_clustering_quality.py @@ -51,7 +51,8 @@ def parse_options(): dest='gt_path', required=True, action='store', type=str, help='Path to file or directory containing ground truth labels') - + parser.add_argument('-s', '--singleton_noise_points', action='store_true', + help='Assign a singleton cluster to each noise point in the clustering result.') args = parser.parse_args() return args @@ -63,7 +64,10 @@ def main(): true_labels = read_label_data(opt.gt_path) try: - eval_clusters(cluster_labels, true_labels) + eval_clusters(cluster_labels, true_labels, False) + contains_noise_points = np.any(cluster_labels == -1) + if opt.singleton_noise_points and contains_noise_points: + eval_clusters(cluster_labels, true_labels, True) except Exception as e: print(f"Error: {e}") exit(1) diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index f6953a9..eff7811 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -92,6 +92,8 @@ def parse_options(): # Use the -y flag to use a YGM calculation instead parser.add_argument('-y', '--ygm_cluster_eval', action='store_true', help='Use YGM cluster evaluation calculation.') + parser.add_argument('-z', '--dummy_cluster_id', action='store_true', + help='Assign a singleton cluster to each noise point in the clustering result.') # For output parser.add_argument('-o', '--output_root_dir', @@ -138,6 +140,9 @@ def parse_options(): parser.add_argument('-Y', '--ygm_evaluator_exe', default=f'{cwd}/tpls/partition-comparison/build/src/clustering_metrics', help='Path to the YGM clustering evaluation executable.') + parser.add_argument('--noise_point_assigner_exe', + default=f'{cwd}/src/clustering/cluster_noise_points', + help='Path to the noise point assigner executable.') # Etc parser.add_argument('-v', '--verbose', action='store_true', @@ -152,6 +157,42 @@ def generate_job_name(): return f"job_{datetime.now().strftime('%Y%m%d_%H%M%S')}" +def add_clustering_evaluation(job_script, cluster_label_file, amst_ds_path, + ground_truth_path, evaluator, + noise_point_assigner_exe, ygm_cluster_eval, + num_tasks_per_node, verbose, + singleton_cluster_to_noise_points): + job_script.write("echo\n") + job_script.write("date\n") + job_script.write(f"echo ================================\n") + job_script.write(f"echo \"Evaluating Clustering\"\n") + job_script.write(f"echo ================================\n") + + if ygm_cluster_eval: + evaluated_label_file = f"{cluster_label_file}.noise_assigned" + assignment_command = (f"{noise_point_assigner_exe} -M " + f"-m {amst_ds_path} " + f"-c {cluster_label_file} " + f"-o {evaluated_label_file}") + add_cmd(assignment_command, job_script) + + job_script.write("echo \"Evaluating Clustering using YGM\"\n") + verbose_flag = " -v" if verbose else "" + evaluation_command = (f"{evaluator}{verbose_flag} " + f"-g {ground_truth_path} " + f"{evaluated_label_file}") + add_srun_cmd(num_tasks_per_node, evaluation_command, job_script) + else: + job_script.write( + "echo \"Evaluating Clustering using python script\"\n") + evaluation_command = (f"python3 {evaluator} " + f"-c {cluster_label_file} " + f"-g {ground_truth_path}") + if singleton_cluster_to_noise_points: + evaluation_command += " -s" + add_cmd(evaluation_command, job_script) + + # Function to generate the batch script for running a benchmark # # Each line in min_cluster_size_set_cmnds is a shell commands to @@ -170,9 +211,11 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, evaluator, ygm_cluster_eval, verbose, ground_truth_path, + singleton_cluster_to_noise_points, min_cluster_size_set_cmnds, min_samples, - input_dnnd_ds_path=''): + noise_point_assigner_exe, + input_dnnd_ds_path='',): create_dir(job_dir) job_script_path = f'{job_dir}/job.sh' @@ -266,20 +309,31 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, # Run the evaluation step if ground_truth_path: - if ygm_cluster_eval: - job_script.write( - f"echo \"Evaluating Clustering using YGM\"\n") - if verbose: - evaluation_command = f"{evaluator} -v -g {ground_truth_path} {cluster_label_file}" - else: - evaluation_command = f"{evaluator} -g {ground_truth_path} {cluster_label_file}" - add_srun_cmd(num_tasks_per_node, evaluation_command, - job_script) - else: - job_script.write( - f"echo \"Evaluating Clustering using python script\"\n") - evaluation_command = f"python3 {evaluator} -c {cluster_label_file} -g {ground_truth_path}" - add_cmd(evaluation_command, job_script) + add_clustering_evaluation( + job_script, cluster_label_file, amst_ds_path, + ground_truth_path, evaluator, noise_point_assigner_exe, + ygm_cluster_eval, num_tasks_per_node, verbose, + singleton_cluster_to_noise_points) + + job_script.write("echo\n") + job_script.write("date\n") + job_script.write(f"echo ================================\n") + job_script.write(f"echo \"Assign clusters to noise points\"\n") + job_script.write(f"echo ================================\n") + # Remove the .txt extension and add .noise_assigned.txt + cluster_label_file_no_noise = f"{cluster_label_file[:-4]}.noise_assigned.txt" + cluster_assign_command = (f"{noise_point_assigner_exe} -M " + f"-m {amst_ds_path} " + f"-c {cluster_label_file} " + f"-o {cluster_label_file_no_noise}") + add_cmd(cluster_assign_command, job_script) + add_clustering_evaluation( + job_script, cluster_label_file_no_noise, amst_ds_path, + ground_truth_path, evaluator, noise_point_assigner_exe, + ygm_cluster_eval, num_tasks_per_node, verbose, + singleton_cluster_to_noise_points) + + job_script.write(f"echo \"\" \n") try_no += 1 @@ -353,8 +407,10 @@ def main(): opts.ygm_cluster_eval, opts.verbose, opts.ground_truth_path, + opts.dummy_cluster_id, min_cluster_size_set_cmnds, opts.min_samples, + opts.noise_point_assigner_exe, opts.input_dnnd_ds_path) print(f"Generated batch script: {job_script}") diff --git a/src/clustering/CMakeLists.txt b/src/clustering/CMakeLists.txt index 34039cb..b604f17 100644 --- a/src/clustering/CMakeLists.txt +++ b/src/clustering/CMakeLists.txt @@ -6,3 +6,8 @@ add_basic_executable(run_hdbscan_clustering run_hdbscan_clustering.cpp) setup_omp_target(run_hdbscan_clustering) setup_metall_target(run_hdbscan_clustering) setup_spdlog_target(run_hdbscan_clustering) + +add_basic_executable(cluster_noise_points cluster_noise_points.cpp) +setup_omp_target(cluster_noise_points) +setup_metall_target(cluster_noise_points) +setup_spdlog_target(cluster_noise_points) \ No newline at end of file diff --git a/src/clustering/cluster_noise_points.cpp b/src/clustering/cluster_noise_points.cpp new file mode 100644 index 0000000..5842981 --- /dev/null +++ b/src/clustering/cluster_noise_points.cpp @@ -0,0 +1,218 @@ +// Assign cluster IDs to noise points by traversing the MST edges. +// Traverse the MST edges from each noise point in BFS manner until a point that +// belongs to a cluster is found. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "../common.hpp" + +using namespace clams; + +template +using map_t = boost::unordered::unordered_flat_map; + +static constexpr id_t k_noise_cluster_id = static_cast(-1); + +struct option { + std::filesystem::path mst_edges_path; + bool metall_mst{false}; + std::filesystem::path cluster_ids_input_path; + std::filesystem::path cluster_ids_out_path; +}; + +void show_help() { + std::cout + << "<>\n" + "Required arguments:\n" + " -m path to a file or a directory that contains input MST.\n" + " -c path to input cluster IDs.\n" + " -o path to output cluster IDs.\n" + "Optional arguments:\n" + " -M If specified, input is a Metall datastore.\n" + " -h Show help." + << std::endl; +} + +void parse_option(int argc, char* argv[], option& opt) { + int opt_char; + while ((opt_char = getopt(argc, argv, "m:c:o:Mh")) != -1) { + switch (opt_char) { + case 'm': + opt.mst_edges_path = optarg; + break; + case 'c': + opt.cluster_ids_input_path = optarg; + break; + case 'o': + opt.cluster_ids_out_path = optarg; + break; + case 'M': + opt.metall_mst = true; + break; + case 'h': + show_help(); + std::exit(EXIT_SUCCESS); + default: + show_help(); + std::exit(EXIT_FAILURE); + } + } +} + +void read_cluster_ids(const std::filesystem::path& input_path, + map_t& point_cluster_map) { + spdlog::info("Reading cluster IDs from {}", input_path.string()); + std::ifstream ifs(input_path); + if (!ifs) { + std::cerr << "Failed to open " << input_path << std::endl; + std::abort(); + } + + std::string line; + while (std::getline(ifs, line)) { + if (line.empty() || line[0] == '#') { + continue; // Skip empty lines and comments + } + std::istringstream iss(line); + id_t point_id, cluster_id; + if (!(iss >> point_id >> cluster_id)) { + std::cerr << "Error parsing line: " << line << std::endl; + std::abort(); + } + point_cluster_map[point_id] = cluster_id; + } +} + +void dump_point_cluster_ids(const map_t& cluster_id, + const std::filesystem::path& output_path) { + std::ofstream ofs(output_path); + if (!ofs) { + std::cerr << "Failed to open " << output_path << std::endl; + std::abort(); + } + + for (const auto& [i, final_cluster_id] : cluster_id) { + ofs << i << "\t" << final_cluster_id; + ofs << "\n"; + } + ofs.close(); + if (!ofs) { + std::cerr << "Failed to write to " << output_path << std::endl; + std::abort(); + } +} + +int main(int argc, char* argv[]) { + option opt; + parse_option(argc, argv, opt); + + map_t> mst; + if (opt.metall_mst) { + spdlog::info("Attaching MST in Metall datastore"); + metall::manager metall_manager(metall::open_read_only, opt.mst_edges_path); + auto* input_mst_edges = + metall_manager.find(metall::unique_instance) + .first; + if (!input_mst_edges) { + spdlog::critical("Failed to find MST edges in Metall datastore at {}", + opt.mst_edges_path.string()); + std::abort(); + } + spdlog::info("Copying MST edges from Metall datastore"); + for (const auto& edge : *input_mst_edges) { + mst[edge.ids[0]].push_back(edge.ids[1]); + mst[edge.ids[1]].push_back(edge.ids[0]); + } + spdlog::info("#of MST edges: {}", mst.size()); + } else { + spdlog::info("Reading MST edges"); + weighted_edge_list_t input_mst_edges; + read_edges(opt.mst_edges_path, input_mst_edges); + assert(opt.min_cluster_size > 1); // Current implementation restriction. + spdlog::info("#of MST edges: {}", input_mst_edges.size()); + for (const auto& edge : input_mst_edges) { + mst[edge.ids[0]].push_back(edge.ids[1]); + mst[edge.ids[1]].push_back(edge.ids[0]); + } + } + + if (mst.empty()) { + spdlog::warn("No MST edges found in the input file or directory: {}", + opt.mst_edges_path.string()); + return EXIT_SUCCESS; + } + + map_t point_cluster_map; + read_cluster_ids(opt.cluster_ids_input_path, point_cluster_map); + spdlog::info("Read {} point cluster IDs from {}", point_cluster_map.size(), + opt.cluster_ids_input_path.string()); + + std::vector point_ids; + point_ids.reserve(point_cluster_map.size()); + for (const auto& [point_id, cluster_id] : point_cluster_map) { + point_ids.push_back(point_id); + } + + spdlog::info( + "Assigning cluster IDs to noise points by traversing the MST edges"); + std::size_t n_noise_points = 0; + std::size_t n_assigned_points = 0; + OMP_DIRECTIVE(parallel for reduction(+ : n_noise_points, n_assigned_points)) + for (size_t i = 0; i < point_ids.size(); ++i) { + const auto point_id = point_ids.at(i); + if (point_cluster_map.at(point_id) != k_noise_cluster_id) { + continue; // Already assigned to a cluster + } + ++n_noise_points; + + std::deque bfs_queue; + bfs_queue.push_back(point_id); + bool found_cluster = false; + + while (!bfs_queue.empty() && !found_cluster) { + const auto current_point_id = bfs_queue.front(); + bfs_queue.pop_front(); + + for (const auto neighbor_id : mst.at(current_point_id)) { + if (point_cluster_map.at(neighbor_id) != k_noise_cluster_id) { + point_cluster_map[point_id] = point_cluster_map.at(neighbor_id); + found_cluster = true; + ++n_assigned_points; + break; + } else { + bfs_queue.push_back(neighbor_id); + } + } + } + + if (!found_cluster) { + spdlog::warn("Point {} could not be assigned to any cluster.", point_id); + } + } + spdlog::info("Finished assigning cluster IDs to noise points"); + spdlog::info("Number of noise points: {}", n_noise_points); + spdlog::info("Number of noise points assigned to clusters: {}", + n_assigned_points); + + dump_point_cluster_ids(point_cluster_map, opt.cluster_ids_out_path); + spdlog::info("Dumped point cluster IDs to {}", + opt.cluster_ids_out_path.string()); + + return EXIT_SUCCESS; +} From 495a90385896df021ebb64be7e172b808d6bc943 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Thu, 13 Aug 2026 07:03:08 -0700 Subject: [PATCH 09/37] Add copy_pm_datastore --- src/tools/CMakeLists.txt | 6 ++- src/tools/copy_pm_datastore.cpp | 78 +++++++++++++++++++++++++++++++++ 2 files changed, 82 insertions(+), 2 deletions(-) create mode 100644 src/tools/copy_pm_datastore.cpp diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt index 0a5e994..68e2516 100644 --- a/src/tools/CMakeLists.txt +++ b/src/tools/CMakeLists.txt @@ -1,6 +1,8 @@ # Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS # Project Developers. See the top-level COPYRIGHT file for details. - add_basic_executable(conv_knng_to_mreach_dist conv_knng_to_mreach_dist.cpp) -setup_saltatlas_target(conv_knng_to_mreach_dist) \ No newline at end of file +setup_saltatlas_target(conv_knng_to_mreach_dist) + +add_basic_executable(copy_pm_datastore copy_pm_datastore.cpp) +setup_saltatlas_target(copy_pm_datastore) \ No newline at end of file diff --git a/src/tools/copy_pm_datastore.cpp b/src/tools/copy_pm_datastore.cpp new file mode 100644 index 0000000..2ff01cd --- /dev/null +++ b/src/tools/copy_pm_datastore.cpp @@ -0,0 +1,78 @@ +// Copyright 2026 Lawrence Livermore National Security, LLC and other Metall +// Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (Apache-2.0 OR MIT) + +// This program copies a Metall datastore created by the metall_mpi_adaptor +// utility to a new location. It is useful for copying a large datastore because +// it can copy the datastore in parallel using multiple MPI processes. + +#include + +#include +#include +#include + +#include + +#include + +namespace { + +bool parse_options(int argc, char **argv, std::filesystem::path &datastore_path, + std::filesystem::path ©_path, bool &overwrite) { + overwrite = false; + int option; + while ((option = ::getopt(argc, argv, "s:t:o")) != -1) { + switch (option) { + case 's': + datastore_path = optarg; + break; + case 't': + copy_path = optarg; + break; + case 'o': + overwrite = true; + break; + default: + return false; + } + } + + return !datastore_path.empty() && !copy_path.empty() && optind == argc; +} + +} // namespace + +int main(int argc, char **argv) { + ygm::comm comm(&argc, &argv); + { + std::filesystem::path datastore_path; + std::filesystem::path copy_path; + bool overwrite = false; + + if (!parse_options(argc, argv, datastore_path, copy_path, overwrite)) { + comm.cerr0() << "Usage: " << argv[0] + << " -s -t " + " [-o (overwrite)]" + << std::endl; + return EXIT_FAILURE; + } + + comm.cout0() << "Copying a Metall datastore from " << datastore_path + << " to " << copy_path << ", overwrite=" << overwrite + << std::endl; + + const bool ret = metall::utility::metall_mpi_adaptor::copy( + datastore_path, copy_path, comm.get_mpi_comm(), overwrite); + comm.cf_barrier(); + + if (ret) { + comm.cout0() << "Done." << std::endl; + } else { + comm.cerr0() << "Failed to copy the Metall datastore." << std::endl; + } + } + + return 0; +} From a7f510db98a7b42e91ed829fb9df80f4c2f1748c Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Thu, 13 Aug 2026 07:16:24 -0700 Subject: [PATCH 10/37] Fix build error in cluster_noise_points.cpp --- src/clustering/cluster_noise_points.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/clustering/cluster_noise_points.cpp b/src/clustering/cluster_noise_points.cpp index 5842981..b36ab14 100644 --- a/src/clustering/cluster_noise_points.cpp +++ b/src/clustering/cluster_noise_points.cpp @@ -144,7 +144,6 @@ int main(int argc, char* argv[]) { spdlog::info("Reading MST edges"); weighted_edge_list_t input_mst_edges; read_edges(opt.mst_edges_path, input_mst_edges); - assert(opt.min_cluster_size > 1); // Current implementation restriction. spdlog::info("#of MST edges: {}", input_mst_edges.size()); for (const auto& edge : input_mst_edges) { mst[edge.ids[0]].push_back(edge.ids[1]); From 726511fbd3cdb8f58d8d3436111f179840e0f02b Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Thu, 13 Aug 2026 09:33:47 -0700 Subject: [PATCH 11/37] Improve run_clams_bench.py --- script/benchmark/run_clams_bench.py | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index eff7811..13b416a 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -159,7 +159,7 @@ def generate_job_name(): def add_clustering_evaluation(job_script, cluster_label_file, amst_ds_path, ground_truth_path, evaluator, - noise_point_assigner_exe, ygm_cluster_eval, + ygm_cluster_eval, num_tasks_per_node, verbose, singleton_cluster_to_noise_points): job_script.write("echo\n") @@ -169,18 +169,11 @@ def add_clustering_evaluation(job_script, cluster_label_file, amst_ds_path, job_script.write(f"echo ================================\n") if ygm_cluster_eval: - evaluated_label_file = f"{cluster_label_file}.noise_assigned" - assignment_command = (f"{noise_point_assigner_exe} -M " - f"-m {amst_ds_path} " - f"-c {cluster_label_file} " - f"-o {evaluated_label_file}") - add_cmd(assignment_command, job_script) - job_script.write("echo \"Evaluating Clustering using YGM\"\n") verbose_flag = " -v" if verbose else "" - evaluation_command = (f"{evaluator}{verbose_flag} " + evaluation_command = (f"{evaluator} {verbose_flag} " f"-g {ground_truth_path} " - f"{evaluated_label_file}") + f"{cluster_label_file}") add_srun_cmd(num_tasks_per_node, evaluation_command, job_script) else: job_script.write( @@ -311,7 +304,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, if ground_truth_path: add_clustering_evaluation( job_script, cluster_label_file, amst_ds_path, - ground_truth_path, evaluator, noise_point_assigner_exe, + ground_truth_path, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, singleton_cluster_to_noise_points) @@ -327,9 +320,11 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, f"-c {cluster_label_file} " f"-o {cluster_label_file_no_noise}") add_cmd(cluster_assign_command, job_script) + + add_clustering_evaluation( job_script, cluster_label_file_no_noise, amst_ds_path, - ground_truth_path, evaluator, noise_point_assigner_exe, + ground_truth_path, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, singleton_cluster_to_noise_points) From 756c207604b0b79cadbabfc605e9e8f852920f1d Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Thu, 3 Sep 2026 20:44:50 -0700 Subject: [PATCH 12/37] Enhance clustering evaluation and correlation analysis - Updated eval_clusters function to include an option to ignore true noise points. - Added new tool for evaluating correlation between kNNG and clusters. --- script/benchmark/clustering_utilities.py | 16 +- .../benchmark/evaluate_clustering_quality.py | 12 +- script/benchmark/run_hdbscan_bench.py | 2 +- src/clustering/cluster_noise_points.cpp | 60 ++----- src/common.hpp | 60 ++++++- src/tools/CMakeLists.txt | 8 +- ...evaluate_correlation_knng_and_clusters.cpp | 156 ++++++++++++++++++ 7 files changed, 251 insertions(+), 63 deletions(-) create mode 100644 src/tools/evaluate_correlation_knng_and_clusters.cpp diff --git a/script/benchmark/clustering_utilities.py b/script/benchmark/clustering_utilities.py index e72f8be..2dbc633 100644 --- a/script/benchmark/clustering_utilities.py +++ b/script/benchmark/clustering_utilities.py @@ -10,17 +10,23 @@ # Modified from: # https://gist.github.com/lmcinnes/24ed5c22c80125be5133811d677eae7b -def eval_clusters(cluster_labels, true_labels, singleton_cluster_to_noise_points=False): +def eval_clusters(cluster_labels, true_labels, singleton_cluster_to_noise_points=False, ignore_true_noise_points=False): max_cluster_id = max(np.max(cluster_labels), np.max(true_labels)) if np.any(true_labels < 0): print("Ground truth labels contain noise points") pct_clustered_gt = (np.sum(true_labels >= 0) / cluster_labels.shape[0]) print(f"GT clustered Points: {pct_clustered_gt * 100:.2f}%") - print( - "Assigning a singleton cluster to each noise point in the ground truth labels") - true_labels = assign_singleton_cluster_to_noise_points(true_labels, max_cluster_id) - max_cluster_id = max(max_cluster_id, np.max(true_labels)) + + if ignore_true_noise_points: + print("Remove true noise points from the evaluation") + print(f"Before filtering: {len(true_labels)} points") + mask = true_labels >= 0 + true_labels = true_labels[mask] + cluster_labels = cluster_labels[mask] + max_cluster_id = max(np.max(cluster_labels), np.max(true_labels)) + print(f"After filtering: {len(true_labels)} points") + if np.any(cluster_labels < 0): # Has noise points clustered_points = (cluster_labels >= 0) diff --git a/script/benchmark/evaluate_clustering_quality.py b/script/benchmark/evaluate_clustering_quality.py index 0cc5813..6dc700a 100644 --- a/script/benchmark/evaluate_clustering_quality.py +++ b/script/benchmark/evaluate_clustering_quality.py @@ -41,7 +41,7 @@ def parse_options(): parser = argparse.ArgumentParser( - description='Evaluate kNN index') + description='Evaluate clustering quality') parser.add_argument('-c', '--cluster', dest='cluster_path', @@ -52,7 +52,10 @@ def parse_options(): required=True, action='store', type=str, help='Path to file or directory containing ground truth labels') parser.add_argument('-s', '--singleton_noise_points', action='store_true', - help='Assign a singleton cluster to each noise point in the clustering result.') + help='Assign a singleton cluster to each noise point in the clustering labels to evaluate.') + parser.add_argument('-T', '--ignore_true_noise_points', action='store_true', + help='Remove true noise points from the evaluation.') + args = parser.parse_args() return args @@ -64,10 +67,7 @@ def main(): true_labels = read_label_data(opt.gt_path) try: - eval_clusters(cluster_labels, true_labels, False) - contains_noise_points = np.any(cluster_labels == -1) - if opt.singleton_noise_points and contains_noise_points: - eval_clusters(cluster_labels, true_labels, True) + eval_clusters(cluster_labels, true_labels, opt.singleton_noise_points, opt.ignore_true_noise_points) except Exception as e: print(f"Error: {e}") exit(1) diff --git a/script/benchmark/run_hdbscan_bench.py b/script/benchmark/run_hdbscan_bench.py index d2454de..88b40d0 100644 --- a/script/benchmark/run_hdbscan_bench.py +++ b/script/benchmark/run_hdbscan_bench.py @@ -13,7 +13,7 @@ import argparse from clustering_utilities import * from script.benchmark.bench_utilities import * -from hdbscan.run_hdbscan import run_hdbscan +from script.benchmark.hdbscan.run_hdbscan import run_hdbscan def parse_options(): diff --git a/src/clustering/cluster_noise_points.cpp b/src/clustering/cluster_noise_points.cpp index b36ab14..faa18db 100644 --- a/src/clustering/cluster_noise_points.cpp +++ b/src/clustering/cluster_noise_points.cpp @@ -27,8 +27,6 @@ using namespace clams; template using map_t = boost::unordered::unordered_flat_map; -static constexpr id_t k_noise_cluster_id = static_cast(-1); - struct option { std::filesystem::path mst_edges_path; bool metall_mst{false}; @@ -75,49 +73,6 @@ void parse_option(int argc, char* argv[], option& opt) { } } -void read_cluster_ids(const std::filesystem::path& input_path, - map_t& point_cluster_map) { - spdlog::info("Reading cluster IDs from {}", input_path.string()); - std::ifstream ifs(input_path); - if (!ifs) { - std::cerr << "Failed to open " << input_path << std::endl; - std::abort(); - } - - std::string line; - while (std::getline(ifs, line)) { - if (line.empty() || line[0] == '#') { - continue; // Skip empty lines and comments - } - std::istringstream iss(line); - id_t point_id, cluster_id; - if (!(iss >> point_id >> cluster_id)) { - std::cerr << "Error parsing line: " << line << std::endl; - std::abort(); - } - point_cluster_map[point_id] = cluster_id; - } -} - -void dump_point_cluster_ids(const map_t& cluster_id, - const std::filesystem::path& output_path) { - std::ofstream ofs(output_path); - if (!ofs) { - std::cerr << "Failed to open " << output_path << std::endl; - std::abort(); - } - - for (const auto& [i, final_cluster_id] : cluster_id) { - ofs << i << "\t" << final_cluster_id; - ofs << "\n"; - } - ofs.close(); - if (!ofs) { - std::cerr << "Failed to write to " << output_path << std::endl; - std::abort(); - } -} - int main(int argc, char* argv[]) { option opt; parse_option(argc, argv, opt); @@ -180,8 +135,10 @@ int main(int argc, char* argv[]) { } ++n_noise_points; - std::deque bfs_queue; + std::deque bfs_queue; + map_t visited; bfs_queue.push_back(point_id); + visited[point_id] = true; bool found_cluster = false; while (!bfs_queue.empty() && !found_cluster) { @@ -189,12 +146,16 @@ int main(int argc, char* argv[]) { bfs_queue.pop_front(); for (const auto neighbor_id : mst.at(current_point_id)) { + if (visited.find(neighbor_id) != visited.end()) { + continue; // Already visited, e.g., the node we came from + } if (point_cluster_map.at(neighbor_id) != k_noise_cluster_id) { point_cluster_map[point_id] = point_cluster_map.at(neighbor_id); found_cluster = true; ++n_assigned_points; break; } else { + visited[neighbor_id] = true; bfs_queue.push_back(neighbor_id); } } @@ -205,9 +166,10 @@ int main(int argc, char* argv[]) { } } spdlog::info("Finished assigning cluster IDs to noise points"); - spdlog::info("Number of noise points: {}", n_noise_points); - spdlog::info("Number of noise points assigned to clusters: {}", - n_assigned_points); + spdlog::info("Number of noise points in the original data: {}", + n_noise_points); + spdlog::info("Number of remaining noise points: {}", + n_noise_points - n_assigned_points); dump_point_cluster_ids(point_cluster_map, opt.cluster_ids_out_path); spdlog::info("Dumped point cluster IDs to {}", diff --git a/src/common.hpp b/src/common.hpp index b75efc6..5c21616 100644 --- a/src/common.hpp +++ b/src/common.hpp @@ -22,6 +22,8 @@ namespace clams { +static constexpr id_t k_noise_cluster_id = static_cast(-1); + inline std::vector find_files( const std::filesystem::path &path) { std::vector files; @@ -159,4 +161,60 @@ inline void read_knng_edges( } } } -} // namespace clams \ No newline at end of file + +/// Reads point-to-cluster assignments from `input_path` into +/// `point_cluster_map`. +/// +/// Each nonempty, noncomment input line must contain a point ID followed by its +/// cluster ID, separated by whitespace. Lines whose first character is `#` are +/// ignored. The output maps each point ID to its cluster ID; if an ID occurs +/// more than once, the last assignment wins. +template +void read_cluster_ids(const std::filesystem::path &input_path, + cluster_id_table_t &point_cluster_map) { + using id_t = typename cluster_id_table_t::key_type; + using cluster_id_t = typename cluster_id_table_t::mapped_type; + + std::ifstream ifs(input_path); + if (!ifs) { + std::cerr << "Failed to open " << input_path << std::endl; + std::abort(); + } + + std::string line; + while (std::getline(ifs, line)) { + if (line.empty() || line[0] == '#') { + continue; // Skip empty lines and comments + } + std::istringstream iss(line); + id_t point_id; + cluster_id_t cluster_id; + if (!(iss >> point_id >> cluster_id)) { + std::cerr << "Error parsing line: " << line << std::endl; + std::abort(); + } + point_cluster_map[point_id] = cluster_id; + } +} + +template +void dump_point_cluster_ids(const cluster_id_table_t &cluster_id, + const std::filesystem::path &output_path) { + std::ofstream ofs(output_path); + if (!ofs) { + std::cerr << "Failed to open " << output_path << std::endl; + std::abort(); + } + + for (const auto &[i, final_cluster_id] : cluster_id) { + ofs << i << "\t" << final_cluster_id; + ofs << "\n"; + } + ofs.close(); + if (!ofs) { + std::cerr << "Failed to write to " << output_path << std::endl; + std::abort(); + } +} + +} // namespace clams diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt index 68e2516..8c9744f 100644 --- a/src/tools/CMakeLists.txt +++ b/src/tools/CMakeLists.txt @@ -5,4 +5,10 @@ add_basic_executable(conv_knng_to_mreach_dist conv_knng_to_mreach_dist.cpp) setup_saltatlas_target(conv_knng_to_mreach_dist) add_basic_executable(copy_pm_datastore copy_pm_datastore.cpp) -setup_saltatlas_target(copy_pm_datastore) \ No newline at end of file +setup_saltatlas_target(copy_pm_datastore) + +add_basic_executable(evaluate_correlation_knng_and_clusters + evaluate_correlation_knng_and_clusters.cpp) +setup_omp_target(evaluate_correlation_knng_and_clusters) +setup_metall_target(evaluate_correlation_knng_and_clusters) +setup_spdlog_target(evaluate_correlation_knng_and_clusters) \ No newline at end of file diff --git a/src/tools/evaluate_correlation_knng_and_clusters.cpp b/src/tools/evaluate_correlation_knng_and_clusters.cpp new file mode 100644 index 0000000..0600bdf --- /dev/null +++ b/src/tools/evaluate_correlation_knng_and_clusters.cpp @@ -0,0 +1,156 @@ +// Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS +// Project Developers. See the top-level COPYRIGHT file for details. + +// Read kNNG and compute how well the kNNG and clusters are correlated. + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "../common.hpp" +#include "../details/shm_graph.hpp" + +namespace omp = metall::utility::omp; + +using cluster_id_table_t = boost::unordered::unordered_flat_map; + +// todo: +// add 'max_k' option (-k): the maximum number of neighbors to consider for +// correlation analysis. +bool parse_option(int argc, char *argv[], + std::filesystem::path &input_knng_path, + std::filesystem::path &input_clusters_path, size_t &max_k) { + input_knng_path.clear(); + + int opt; + while ((opt = ::getopt(argc, argv, "g:c:k:")) != -1) { + switch (opt) { + case 'g': { + input_knng_path = std::filesystem::path(optarg); + break; + } + case 'c': { + input_clusters_path = std::filesystem::path(optarg); + break; + } + case 'k': { + max_k = std::stoul(optarg); + break; + } + default: { + std::cerr << "Unknown option: " << opt << std::endl; + return false; + } + } + } + + if (input_knng_path.empty()) { + std::cerr << "No input kNNG path is specified" << std::endl; + return false; + } + + if (input_clusters_path.empty()) { + std::cerr << "No input clusters path is specified" << std::endl; + return false; + } + + return true; +} + +int main(int argc, char *argv[]) { + std::filesystem::path input_knng_path; + std::filesystem::path input_clusters_path; + size_t max_k = 100; // default value + + if (!parse_option(argc, argv, input_knng_path, input_clusters_path, max_k)) { + return EXIT_FAILURE; + } + + const auto knng_files = clams::find_files(input_knng_path); + + clams::shm_graph_t graph; + + std::cout << "Read knng" << std::endl; + clams::read_knng(knng_files, graph); + std::cout << "#of points: " << graph.num_keys() << std::endl; + std::cout << "#of neighbors: " << graph.num_values() << std::endl; + + cluster_id_table_t cluster_id_table; + clams::read_cluster_ids(input_clusters_path, cluster_id_table); + std::cout << "#of points with clusters: " << cluster_id_table.size() + << std::endl; + const auto &cluster_ids = cluster_id_table; + + size_t max_n_neighbors = 0; + std::vector point_ids; + point_ids.reserve(graph.num_keys()); + for (auto itr = graph.keys_begin(); itr != graph.keys_end(); ++itr) { + point_ids.push_back(itr->first); + max_n_neighbors = std::max(max_n_neighbors, graph.num_values(itr->first)); + } + std::cout << "max_n_neighbors: " << max_n_neighbors << std::endl; + + size_t n_neighbors_with_same_cluster = 0; + size_t n_neighbors_with_different_cluster = 0; + size_t n_missing_neighbors = 0; + size_t n_missing_cluster_ids = 0; + size_t n_noise_cluster_points = 0; + std::cout + << "k\tsame_cluster\tdifferent_cluster\tcorrelation_rate(%)\tmissing_" + "neighbors\tmissing_cluster_ids" + << std::endl; + for (size_t k = 0; k < std::min(max_n_neighbors, max_k); ++k) { + OMP_DIRECTIVE(parallel for reduction(+ : n_neighbors_with_same_cluster, n_neighbors_with_different_cluster, n_missing_neighbors, n_missing_cluster_ids)) + for (size_t pi = 0; pi < point_ids.size(); ++pi) { + const auto pid = point_ids[pi]; + if (graph.num_values(pid) <= k) { + ++n_missing_neighbors; + continue; + } + + const auto neighbor_id = (graph.values_begin(pid) + k)->first; + const auto point_cluster_itr = cluster_ids.find(pid); + const auto neighbor_cluster_itr = cluster_ids.find(neighbor_id); + if (point_cluster_itr == cluster_ids.end() || + neighbor_cluster_itr == cluster_ids.end()) { + ++n_missing_cluster_ids; + continue; + } + const auto point_cluster_id = point_cluster_itr->second; + const auto neighbor_cluster_id = neighbor_cluster_itr->second; + if (point_cluster_id == clams::k_noise_cluster_id || + neighbor_cluster_id == clams::k_noise_cluster_id) { + ++n_noise_cluster_points; + continue; + } + + if (point_cluster_itr->second == neighbor_cluster_itr->second) { + ++n_neighbors_with_same_cluster; + } else { + ++n_neighbors_with_different_cluster; + } + } + // Show statistics for each k + const auto n_neighbors_with_clusters = + n_neighbors_with_same_cluster + n_neighbors_with_different_cluster; + const auto correlation_rate = + n_neighbors_with_clusters == 0 + ? 0.0 + : static_cast(n_neighbors_with_same_cluster) / + n_neighbors_with_clusters; + std::cout << k + 1 << "\t" << n_neighbors_with_same_cluster << "\t" + << n_neighbors_with_different_cluster << "\t" + << correlation_rate * 100.0 << "\t" << n_missing_neighbors << "\t" + << n_missing_cluster_ids << std::endl; + } + + return EXIT_SUCCESS; +} From 0c57721e92669327758ecc9e1ebade84d768a490 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Thu, 3 Sep 2026 21:58:59 -0700 Subject: [PATCH 13/37] Enhance documentation for kNNG functions and correlation evaluation --- src/common.hpp | 4 ++ src/details/shm_graph.hpp | 37 ++++++++++++------- ...evaluate_correlation_knng_and_clusters.cpp | 8 +++- 3 files changed, 34 insertions(+), 15 deletions(-) diff --git a/src/common.hpp b/src/common.hpp index 5c21616..705e429 100644 --- a/src/common.hpp +++ b/src/common.hpp @@ -79,6 +79,10 @@ inline void read_edges(const std::filesystem::path &path, } } +/// \brief Read a k-nearest-neighbor graph (kNNG) from files and store it as +/// an edge list. +/// \param knng_files A list of knng files. +/// \param graph A graph to store the knng. inline void read_knng_edges( const std::vector &knng_files, weighted_edge_list_t &edges) { diff --git a/src/details/shm_graph.hpp b/src/details/shm_graph.hpp index 7aff55e..cc682ac 100644 --- a/src/details/shm_graph.hpp +++ b/src/details/shm_graph.hpp @@ -1,7 +1,6 @@ // Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS // Project Developers. See the top-level COPYRIGHT file for details. - #pragma once #include @@ -12,21 +11,31 @@ #include #include -#include "multithread_adjacency_list.hpp" #include "data_types.hpp" +#include "multithread_adjacency_list.hpp" namespace clams { -using shm_graph_t = multithread_adjacency_list>; +using shm_graph_t = + multithread_adjacency_list>; /// \brief Read knng files. /// \param knng_files A list of knng files. /// \param graph A graph to store the knng. +/// \details Each knng file is expected to have the following format: +/// - Each even numbered line (0-based) contains a source vertex ID followed +/// by its neighbor vertex IDs, separated by whitespace. +/// - Each odd numbered line contains the corresponding distances to the +/// neighbors, separated by whitespace. +/// - The first neighbor ID in each even numbered line is the source vertex ID +/// itself, and the corresponding distance in the odd numbered line is 0.0 +/// (dummy distance). +/// - The number of neighbor IDs in each even numbered line must match the +/// number of distances in the corresponding odd numbered void read_knng(const std::vector &knng_files, - shm_graph_t &graph) { + shm_graph_t &graph) { OMP_DIRECTIVE(parallel for) for (std::size_t fno = 0; fno < knng_files.size(); ++fno) { - const auto &file = knng_files[fno]; + const auto &file = knng_files[fno]; std::ifstream ifs(file); if (!ifs.is_open()) { std::cerr << "Cannot open file: " << file << std::endl; @@ -39,7 +48,7 @@ void read_knng(const std::vector &knng_files, { std::getline(ifs, line); std::istringstream iss(line); - id_t buf; + id_t buf; while (iss >> buf) { ids.push_back(buf); } @@ -49,7 +58,7 @@ void read_knng(const std::vector &knng_files, { std::getline(ifs, line); std::istringstream iss(line); - distance_t buf; + distance_t buf; while (iss >> buf) { dists.push_back(buf); } @@ -92,7 +101,7 @@ void dump_graph(const shm_graph_t &graph, const std::filesystem::path &file) { for (auto eit = graph.values_begin(vid); eit != graph.values_end(vid); ++eit) { ofs << " " << eit->first; - } + } ofs << std::endl; // Distances @@ -100,7 +109,7 @@ void dump_graph(const shm_graph_t &graph, const std::filesystem::path &file) { for (auto eit = graph.values_begin(vid); eit != graph.values_end(vid); ++eit) { ofs << eit->second << " "; - } + } ofs << std::endl; } } @@ -118,10 +127,10 @@ void make_undirected_graph(shm_graph_t &graph) { const auto vid = vertices[i]; for (auto eit = graph.values_begin(vid); eit != graph.values_end(vid); ++eit) { - const auto nid = eit->first; + const auto nid = eit->first; const auto dist = eit->second; r_graph.add(nid, std::make_pair(vid, dist)); - } + } } // Merge the reversed graph to the original graph @@ -134,7 +143,7 @@ void make_undirected_graph(shm_graph_t &graph) { for (auto eit = r_graph.values_begin(vid); eit != r_graph.values_end(vid); ++eit) { graph.add(vid, *eit); - } + } } } -} \ No newline at end of file +} // namespace clams \ No newline at end of file diff --git a/src/tools/evaluate_correlation_knng_and_clusters.cpp b/src/tools/evaluate_correlation_knng_and_clusters.cpp index 0600bdf..35178a4 100644 --- a/src/tools/evaluate_correlation_knng_and_clusters.cpp +++ b/src/tools/evaluate_correlation_knng_and_clusters.cpp @@ -1,7 +1,13 @@ // Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS // Project Developers. See the top-level COPYRIGHT file for details. -// Read kNNG and compute how well the kNNG and clusters are correlated. +// This program evaluates how strongly a k-nearest-neighbor graph (kNNG) agrees +// with clustering result. For each point, it checks whether its nearest k +// neighbors belong to the same cluster as the point itself. It reports the +// correlation rate for each k, which is the fraction of neighbors that are in +// the same cluster as the point. A high correlation rate indicates that the +// clustering result is consistent with the local neighborhood structure of the +// data. #include From b4f47fd1dfa8f190d7c5bec2f7c866baa4e1695a Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Thu, 3 Sep 2026 23:37:18 -0700 Subject: [PATCH 14/37] Refactor clustering evaluation output and streamline echo statements in job script --- script/benchmark/clustering_utilities.py | 16 -------- script/benchmark/run_clams_bench.py | 52 +++++++----------------- 2 files changed, 14 insertions(+), 54 deletions(-) diff --git a/script/benchmark/clustering_utilities.py b/script/benchmark/clustering_utilities.py index 2ccf18e..91e85f4 100644 --- a/script/benchmark/clustering_utilities.py +++ b/script/benchmark/clustering_utilities.py @@ -33,21 +33,6 @@ def eval_clusters(cluster_labels, true_labels, singleton_cluster_to_noise_points pct_clustered = (np.sum(clustered_points) / cluster_labels.shape[0]) print(f"Cluster coverage (%): {pct_clustered * 100:.2f}") - if singleton_cluster_to_noise_points: - print( - "Assigning a singleton cluster to each noise point in the clustering result") - cluster_labels = assign_singleton_cluster_to_noise_points(cluster_labels, max_cluster_id) - ari = adjusted_rand_score(true_labels, cluster_labels) - ami = adjusted_mutual_info_score(true_labels, cluster_labels) - else: - ari = adjusted_rand_score(true_labels[clustered_points], - cluster_labels[clustered_points]) - ami = adjusted_mutual_info_score(true_labels[clustered_points], - cluster_labels[clustered_points]) - # sil = silhouette_score(raw_data[clustered_points], cluster_labels[clustered_points]) - pct_clustered = (np.sum(clustered_points) / cluster_labels.shape[0]) - print(f"Cluster coverage (%): {pct_clustered * 100:.2f}") - if singleton_cluster_to_noise_points: print( "Assigning a singleton cluster to each noise point in the clustering result") @@ -71,7 +56,6 @@ def eval_clusters(cluster_labels, true_labels, singleton_cluster_to_noise_points # Assign a cluster ID to each noise point -def assign_singleton_cluster_to_noise_points(cluster_labels, noise_id_offset): def assign_singleton_cluster_to_noise_points(cluster_labels, noise_id_offset): new_labels = cluster_labels.copy() diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index c67a5c1..9e3d45e 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -159,18 +159,18 @@ def generate_job_name(): time.sleep(2) return f"job_{datetime.now().strftime('%Y%m%d_%H%M%S')}" +def echo_stage_name(job_script, stage_name): + job_script.write("echo\n") + job_script.write("date\n") + job_script.write(f"echo ================================\n") + job_script.write(f"echo \"{stage_name}\"\n") + job_script.write(f"echo ================================\n") def add_clustering_evaluation(job_script, cluster_label_file, amst_ds_path, ground_truth_path, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, singleton_cluster_to_noise_points): - job_script.write("echo\n") - job_script.write("date\n") - job_script.write(f"echo ================================\n") - job_script.write(f"echo \"Evaluating Clustering\"\n") - job_script.write(f"echo ================================\n") - if ygm_cluster_eval: job_script.write("echo \"Evaluating Clustering using YGM\"\n") verbose_flag = " -v" if verbose else "" @@ -227,9 +227,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, job_script.write("echo\n") job_script.write("date\n") if len(input_dnnd_ds_path) == 0: - job_script.write(f"echo ================================\n") - job_script.write(f"echo \"Building KNNG\"\n") - job_script.write(f"echo ================================\n") + echo_stage_name(job_script, "Building KNNG") dnnd_ds_path = f"{work_dir}/dnnd_pm_datastore" dnnd_batch_size = 2 ** 25 verbose_flag = '-v' if verbose else '' @@ -237,8 +235,6 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, add_srun_cmd(num_tasks_per_node, dnnd_command, job_script) if backup_knng: dnnd_ds_path_backup = f"{dnnd_ds_path}_backup" - job_script.write(f"echo\n") - job_script.write("date\n") job_script.write(f"echo \"Backing up KNNG datastore\"\n") backup_knng_command = f"cp -r {dnnd_ds_path} {dnnd_ds_path_backup}" add_cmd(backup_knng_command, job_script) @@ -248,22 +244,14 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, dnnd_ds_path = input_dnnd_ds_path # Connect the CCs - job_script.write("echo\n") - job_script.write("date\n") - job_script.write(f"echo ================================\n") - job_script.write(f"echo \"Running MFC\"\n") - job_script.write(f"echo ================================\n") + echo_stage_name(job_script, "Connecting the CCs using MFC") mfc_command = f"{mfc_exe} -d {dnnd_ds_path} -f {distance_func}" add_srun_cmd(num_tasks_per_node, mfc_command, job_script) # Convert to core distance # TODO: Implement if False and min_samples > 0: - job_script.write("echo\n") - job_script.write("date\n") - job_script.write(f"echo ================================\n") - job_script.write(f"echo Convert to core distance kNNG\n") - job_script.write(f"echo ================================\n") + echo_stage_name(job_script, "Convert to core distance kNNG") knng_coredist_dir = f"{work_dir}/knng_coredist/" add_cmd(f'mkdir -p {knng_coredist_dir}', job_script) conv2coredist_cmd = f"./src/conv_knng_to_core_dist -i {dnnd_ds_path} -o {knng_coredist_dir}/knng.txt -m {min_samples}" @@ -272,11 +260,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, try_no = 0 for amst_approx_bound in amst_approx_bound_list: # Run the AMST step - job_script.write("echo\n") - job_script.write("date\n") - job_script.write(f"echo ================================\n") - job_script.write(f"echo \"Running AMST, approx bound = {amst_approx_bound}\"\n") - job_script.write(f"echo ================================\n") + echo_stage_name(job_script, f"Running AMST, approx bound = {amst_approx_bound}") amst_ds_path = f"{work_dir}/amst_pm_datastore_a{amst_approx_bound}" amst_command = f"{amst_exe} -d {dnnd_ds_path} -p {amst_ds_path} -e {amst_approx_bound}" add_srun_cmd(num_tasks_per_node, amst_command, job_script) @@ -286,12 +270,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, # Set the min cluster size environment variable for this run add_cmd(set_cmd, job_script, False, False) - job_script.write("echo\n") - job_script.write("date\n") - job_script.write(f"echo ================================\n") - job_script.write(f"echo \"Running CLAMS-HDBSCAN, min cluster size = ${{MIN_CLUSTER_SIZE}}\"\n") - job_script.write(f"echo ================================\n") - + echo_stage_name(job_script, f"Running CLAMS-HDBSCAN") job_script.write( f"echo \"Min cluster size ${{MIN_CLUSTER_SIZE}}\"\n") if distributed_hdbscan: @@ -316,17 +295,14 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, # Run the evaluation step if ground_truth_path: + echo_stage_name(job_script, "Evaluating Clustering Results") add_clustering_evaluation( job_script, cluster_label_file, amst_ds_path, ground_truth_path, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, singleton_cluster_to_noise_points) - job_script.write("echo\n") - job_script.write("date\n") - job_script.write(f"echo ================================\n") - job_script.write(f"echo \"Assign clusters to noise points\"\n") - job_script.write(f"echo ================================\n") + echo_stage_name(job_script, "Assign clusters to noise points") # Remove the .txt extension and add .noise_assigned.txt cluster_label_file_no_noise = f"{cluster_label_file[:-4]}.noise_assigned.txt" cluster_assign_command = (f"{noise_point_assigner_exe} -M " @@ -335,7 +311,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, f"-o {cluster_label_file_no_noise}") add_cmd(cluster_assign_command, job_script) - + echo_stage_name(job_script, "Evaluate clustering results after assigning clusters to noise points") add_clustering_evaluation( job_script, cluster_label_file_no_noise, amst_ds_path, ground_truth_path, evaluator, From c2000683987c85d95b94fefa2d367bc43af6d51a Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 4 Sep 2026 01:42:40 -0700 Subject: [PATCH 15/37] Refactor clustering execution logic into a dedicated function for improved readability and maintainability --- script/benchmark/run_clams_bench.py | 115 +++++++++++++++------------- 1 file changed, 62 insertions(+), 53 deletions(-) diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index 9e3d45e..fd9c0bd 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -188,6 +188,62 @@ def add_clustering_evaluation(job_script, cluster_label_file, amst_ds_path, evaluation_command += " -s" add_cmd(evaluation_command, job_script) +def run_clustering(job_script, set_cmd, work_dir, amst_approx_bound, + amst_ds_path, clustering_exe, distributed_hdbscan, + evaluator, ygm_cluster_eval, num_tasks_per_node, + verbose, ground_truth_path, + singleton_cluster_to_noise_points, + noise_point_assigner_exe): + # Set the min cluster size environment variable for this run + add_cmd(set_cmd, job_script, False, False) + + echo_stage_name(job_script, "Running CLAMS-HDBSCAN") + job_script.write(f"echo \"Min cluster size ${{MIN_CLUSTER_SIZE}}\"\n") + if distributed_hdbscan: + cluster_label_file = f"{work_dir}/cluster_labels_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}/" + cluster_tree_file = f"{work_dir}/cluster_tree_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}/" + verbose_flag = '-v' if verbose else '' + hpc_clustering_command = (f"{clustering_exe} {verbose_flag} -i {amst_ds_path} -M " + f" -m ${{MIN_CLUSTER_SIZE}} " + f" -o {cluster_label_file} " + f" -c {cluster_tree_file} " + f" -n {num_tasks_per_node}") + add_srun_cmd(num_tasks_per_node, hpc_clustering_command, job_script) + else: + cluster_label_file = f"{work_dir}/cluster_labels_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}.txt" + cluster_tree_file = f"{work_dir}/cluster_tree_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}.txt" + hpc_clustering_command = (f"{clustering_exe} -i {amst_ds_path} -M " + f" -m ${{MIN_CLUSTER_SIZE}} " + f" -o {cluster_label_file} " + f" -c {cluster_tree_file} " + f" -P ") + add_cmd(hpc_clustering_command, job_script) + + if ground_truth_path: + echo_stage_name(job_script, "Evaluating Clustering Results") + add_clustering_evaluation( + job_script, cluster_label_file, amst_ds_path, ground_truth_path, + evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, + singleton_cluster_to_noise_points) + + echo_stage_name(job_script, "Assign clusters to noise points") + cluster_label_file_no_noise = f"{cluster_label_file[:-4]}.noise_assigned.txt" + cluster_assign_command = (f"{noise_point_assigner_exe} -M " + f"-m {amst_ds_path} " + f"-c {cluster_label_file} " + f"-o {cluster_label_file_no_noise}") + add_cmd(cluster_assign_command, job_script) + + echo_stage_name( + job_script, + "Evaluate clustering results after assigning clusters to noise points") + add_clustering_evaluation( + job_script, cluster_label_file_no_noise, amst_ds_path, + ground_truth_path, evaluator, ygm_cluster_eval, + num_tasks_per_node, verbose, singleton_cluster_to_noise_points) + + job_script.write("echo \"\" \n") + # Function to generate the batch script for running a benchmark # @@ -267,59 +323,12 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, # Run the HPC Clustering step for set_cmd in min_cluster_size_set_cmnds: - # Set the min cluster size environment variable for this run - add_cmd(set_cmd, job_script, False, False) - - echo_stage_name(job_script, f"Running CLAMS-HDBSCAN") - job_script.write( - f"echo \"Min cluster size ${{MIN_CLUSTER_SIZE}}\"\n") - if distributed_hdbscan: - cluster_label_file = f"{work_dir}/cluster_labels_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}/" - cluster_tree_file = f"{work_dir}/cluster_tree_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}/" - verbose_flag = '-v' if verbose else '' - hpc_clustering_command = (f"{clustering_exe} {verbose_flag} -i {amst_ds_path} -M " - f" -m ${{MIN_CLUSTER_SIZE}} " - f" -o {cluster_label_file} " - f" -c {cluster_tree_file} " - f" -n {num_tasks_per_node}") - add_srun_cmd(num_tasks_per_node, hpc_clustering_command, job_script) - else: - cluster_label_file = f"{work_dir}/cluster_labels_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}.txt" - cluster_tree_file = f"{work_dir}/cluster_tree_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}.txt" - hpc_clustering_command = (f"{clustering_exe} -i {amst_ds_path} -M " - f" -m ${{MIN_CLUSTER_SIZE}} " - f" -o {cluster_label_file} " - f" -c {cluster_tree_file} " - f" -P ") - add_cmd(hpc_clustering_command, job_script) - - # Run the evaluation step - if ground_truth_path: - echo_stage_name(job_script, "Evaluating Clustering Results") - add_clustering_evaluation( - job_script, cluster_label_file, amst_ds_path, - ground_truth_path, evaluator, - ygm_cluster_eval, num_tasks_per_node, verbose, - singleton_cluster_to_noise_points) - - echo_stage_name(job_script, "Assign clusters to noise points") - # Remove the .txt extension and add .noise_assigned.txt - cluster_label_file_no_noise = f"{cluster_label_file[:-4]}.noise_assigned.txt" - cluster_assign_command = (f"{noise_point_assigner_exe} -M " - f"-m {amst_ds_path} " - f"-c {cluster_label_file} " - f"-o {cluster_label_file_no_noise}") - add_cmd(cluster_assign_command, job_script) - - echo_stage_name(job_script, "Evaluate clustering results after assigning clusters to noise points") - add_clustering_evaluation( - job_script, cluster_label_file_no_noise, amst_ds_path, - ground_truth_path, evaluator, - ygm_cluster_eval, num_tasks_per_node, verbose, - singleton_cluster_to_noise_points) - - - job_script.write(f"echo \"\" \n") + run_clustering( + job_script, set_cmd, work_dir, amst_approx_bound, + amst_ds_path, clustering_exe, distributed_hdbscan, + evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, + ground_truth_path, singleton_cluster_to_noise_points, + noise_point_assigner_exe) try_no += 1 # If the file was not created, return an error From 2f02c26a2775c8bd3b207747141cfc2ecfc92704 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 4 Sep 2026 02:15:05 -0700 Subject: [PATCH 16/37] Minor bugfix in cluster_noise_points.cpp --- src/clustering/cluster_noise_points.cpp | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/clustering/cluster_noise_points.cpp b/src/clustering/cluster_noise_points.cpp index faa18db..1cf891d 100644 --- a/src/clustering/cluster_noise_points.cpp +++ b/src/clustering/cluster_noise_points.cpp @@ -1,3 +1,6 @@ +// Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS +// Project Developers. See the top-level COPYRIGHT file for details. + // Assign cluster IDs to noise points by traversing the MST edges. // Traverse the MST edges from each noise point in BFS manner until a point that // belongs to a cluster is found. @@ -77,7 +80,7 @@ int main(int argc, char* argv[]) { option opt; parse_option(argc, argv, opt); - map_t> mst; + map_t> mst_graph; if (opt.metall_mst) { spdlog::info("Attaching MST in Metall datastore"); metall::manager metall_manager(metall::open_read_only, opt.mst_edges_path); @@ -89,24 +92,24 @@ int main(int argc, char* argv[]) { opt.mst_edges_path.string()); std::abort(); } + spdlog::info("#of MST edges: {}", input_mst_edges.size()); spdlog::info("Copying MST edges from Metall datastore"); for (const auto& edge : *input_mst_edges) { - mst[edge.ids[0]].push_back(edge.ids[1]); - mst[edge.ids[1]].push_back(edge.ids[0]); + mst_graph[edge.ids[0]].push_back(edge.ids[1]); + mst_graph[edge.ids[1]].push_back(edge.ids[0]); } - spdlog::info("#of MST edges: {}", mst.size()); } else { spdlog::info("Reading MST edges"); weighted_edge_list_t input_mst_edges; read_edges(opt.mst_edges_path, input_mst_edges); spdlog::info("#of MST edges: {}", input_mst_edges.size()); for (const auto& edge : input_mst_edges) { - mst[edge.ids[0]].push_back(edge.ids[1]); - mst[edge.ids[1]].push_back(edge.ids[0]); + mst_graph[edge.ids[0]].push_back(edge.ids[1]); + mst_graph[edge.ids[1]].push_back(edge.ids[0]); } } - if (mst.empty()) { + if (mst_graph.empty()) { spdlog::warn("No MST edges found in the input file or directory: {}", opt.mst_edges_path.string()); return EXIT_SUCCESS; @@ -145,7 +148,7 @@ int main(int argc, char* argv[]) { const auto current_point_id = bfs_queue.front(); bfs_queue.pop_front(); - for (const auto neighbor_id : mst.at(current_point_id)) { + for (const auto neighbor_id : mst_graph.at(current_point_id)) { if (visited.find(neighbor_id) != visited.end()) { continue; // Already visited, e.g., the node we came from } From c030037bf1a537aeb19e347ea59296f2e8261337 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 4 Sep 2026 02:36:18 -0700 Subject: [PATCH 17/37] Minor bugfix in cluster_noise_points.cpp --- src/clustering/cluster_noise_points.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/clustering/cluster_noise_points.cpp b/src/clustering/cluster_noise_points.cpp index 1cf891d..8a9f11e 100644 --- a/src/clustering/cluster_noise_points.cpp +++ b/src/clustering/cluster_noise_points.cpp @@ -92,7 +92,7 @@ int main(int argc, char* argv[]) { opt.mst_edges_path.string()); std::abort(); } - spdlog::info("#of MST edges: {}", input_mst_edges.size()); + spdlog::info("#of MST edges: {}", input_mst_edges->size()); spdlog::info("Copying MST edges from Metall datastore"); for (const auto& edge : *input_mst_edges) { mst_graph[edge.ids[0]].push_back(edge.ids[1]); @@ -117,7 +117,7 @@ int main(int argc, char* argv[]) { map_t point_cluster_map; read_cluster_ids(opt.cluster_ids_input_path, point_cluster_map); - spdlog::info("Read {} point cluster IDs from {}", point_cluster_map.size(), + spdlog::info("Read {} points' cluster IDs from {}", point_cluster_map.size(), opt.cluster_ids_input_path.string()); std::vector point_ids; From 4dda3cae9405fd18ff44b6325db71b4e6575705d Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 4 Sep 2026 03:09:28 -0700 Subject: [PATCH 18/37] Brush up benchmark scripts --- script/benchmark/clustering_utilities.py | 29 ++++++++++++------------ src/clustering/cluster_noise_points.cpp | 4 ++++ 2 files changed, 19 insertions(+), 14 deletions(-) diff --git a/script/benchmark/clustering_utilities.py b/script/benchmark/clustering_utilities.py index 91e85f4..8ef7b79 100644 --- a/script/benchmark/clustering_utilities.py +++ b/script/benchmark/clustering_utilities.py @@ -16,23 +16,18 @@ def eval_clusters(cluster_labels, true_labels, singleton_cluster_to_noise_points if np.any(true_labels < 0): print("Ground truth labels contain noise points") pct_clustered_gt = (np.sum(true_labels >= 0) / cluster_labels.shape[0]) - print(f"GT clustered Points: {pct_clustered_gt * 100:.2f}%") + print(f"Ground truth cluster coverage (%): {pct_clustered_gt * 100:.2f}%") if ignore_true_noise_points: - print("Remove true noise points from the evaluation") - print(f"Before filtering: {len(true_labels)} points") + print("Remove noise points in the ground truth from the evaluation") + print(f"Before filtering: {len(true_labels)} ground truth points") mask = true_labels >= 0 true_labels = true_labels[mask] cluster_labels = cluster_labels[mask] max_cluster_id = max(np.max(cluster_labels), np.max(true_labels)) - print(f"After filtering: {len(true_labels)} points") - + print(f"After filtering: {len(true_labels)} ground truth points") if np.any(cluster_labels < 0): # Has noise points - clustered_points = (cluster_labels >= 0) - pct_clustered = (np.sum(clustered_points) / cluster_labels.shape[0]) - print(f"Cluster coverage (%): {pct_clustered * 100:.2f}") - if singleton_cluster_to_noise_points: print( "Assigning a singleton cluster to each noise point in the clustering result") @@ -40,11 +35,15 @@ def eval_clusters(cluster_labels, true_labels, singleton_cluster_to_noise_points ari = adjusted_rand_score(true_labels, cluster_labels) ami = adjusted_mutual_info_score(true_labels, cluster_labels) else: - ari = adjusted_rand_score(true_labels[clustered_points], - cluster_labels[clustered_points]) - ami = adjusted_mutual_info_score(true_labels[clustered_points], - cluster_labels[clustered_points]) - # sil = silhouette_score(raw_data[clustered_points], cluster_labels[clustered_points]) + print("Noise points are ignored in the evaluation") + non_noise_points = (cluster_labels >= 0) + pct_clustered = (np.sum(non_noise_points) / cluster_labels.shape[0]) + print(f"Cluster coverage (%): {pct_clustered * 100:.2f}") + ari = adjusted_rand_score(true_labels[non_noise_points], + cluster_labels[non_noise_points]) + ami = adjusted_mutual_info_score(true_labels[non_noise_points], + cluster_labels[non_noise_points]) + # sil = silhouette_score(raw_data[non_noise_points], cluster_labels[non_noise_points]) else: ari = adjusted_rand_score(true_labels, cluster_labels) ami = adjusted_mutual_info_score(true_labels, cluster_labels) @@ -66,6 +65,8 @@ def assign_singleton_cluster_to_noise_points(cluster_labels, noise_id_offset): new_labels[i] = cnt_noise + noise_id_offset new_labels[i] = cnt_noise + noise_id_offset + print(f"Assigned singleton clusters to {cnt_noise} noise points") + return new_labels diff --git a/src/clustering/cluster_noise_points.cpp b/src/clustering/cluster_noise_points.cpp index 8a9f11e..bf7c0a8 100644 --- a/src/clustering/cluster_noise_points.cpp +++ b/src/clustering/cluster_noise_points.cpp @@ -130,6 +130,9 @@ int main(int argc, char* argv[]) { "Assigning cluster IDs to noise points by traversing the MST edges"); std::size_t n_noise_points = 0; std::size_t n_assigned_points = 0; + // NOTE: this algorithm is not determinstic because threads update the shared + // point_cluster_map concurrently. + // We employ this algorithm because it is simple and fast. OMP_DIRECTIVE(parallel for reduction(+ : n_noise_points, n_assigned_points)) for (size_t i = 0; i < point_ids.size(); ++i) { const auto point_id = point_ids.at(i); @@ -153,6 +156,7 @@ int main(int argc, char* argv[]) { continue; // Already visited, e.g., the node we came from } if (point_cluster_map.at(neighbor_id) != k_noise_cluster_id) { + // Found a neighbor that belongs to a cluster point_cluster_map[point_id] = point_cluster_map.at(neighbor_id); found_cluster = true; ++n_assigned_points; From 15acaa44984df3672f1e53ca0477f04c46ab5412 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 4 Sep 2026 07:00:24 -0700 Subject: [PATCH 19/37] Brush up benchmark scripts --- script/benchmark/clustering_utilities.py | 23 +++--- .../extract_clams_clustering_qualities.sh | 49 ++++++++++++ .../extract_clams_clustering_time.sh | 79 +++++++++++++++++++ .../extract_hpc_clustering_bench_result.sh | 19 ----- script/benchmark/run_clams_bench.py | 34 ++++++-- src/clustering/cluster_noise_points.cpp | 8 +- 6 files changed, 175 insertions(+), 37 deletions(-) create mode 100755 script/benchmark/extract_clams_clustering_qualities.sh create mode 100755 script/benchmark/extract_clams_clustering_time.sh delete mode 100644 script/benchmark/extract_hpc_clustering_bench_result.sh diff --git a/script/benchmark/clustering_utilities.py b/script/benchmark/clustering_utilities.py index 8ef7b79..6cb80ae 100644 --- a/script/benchmark/clustering_utilities.py +++ b/script/benchmark/clustering_utilities.py @@ -11,6 +11,8 @@ # Modified from: # https://gist.github.com/lmcinnes/24ed5c22c80125be5133811d677eae7b def eval_clusters(cluster_labels, true_labels, singleton_cluster_to_noise_points=False, ignore_true_noise_points=False): + print(f"Assigning singleton clusters to noise points: {singleton_cluster_to_noise_points}") + print(f"Ignore true noise points in the ground truth: {ignore_true_noise_points}") max_cluster_id = max(np.max(cluster_labels), np.max(true_labels)) if np.any(true_labels < 0): @@ -27,7 +29,11 @@ def eval_clusters(cluster_labels, true_labels, singleton_cluster_to_noise_points max_cluster_id = max(np.max(cluster_labels), np.max(true_labels)) print(f"After filtering: {len(true_labels)} ground truth points") - if np.any(cluster_labels < 0): # Has noise points + non_noise_points_mask = (cluster_labels >= 0) + pct_clustered = (np.sum(non_noise_points_mask) / cluster_labels.shape[0]) + print(f"Cluster coverage (%): {pct_clustered * 100:.2f}") + + if len(non_noise_points_mask) < cluster_labels.shape[0]: # Has noise points if singleton_cluster_to_noise_points: print( "Assigning a singleton cluster to each noise point in the clustering result") @@ -36,19 +42,16 @@ def eval_clusters(cluster_labels, true_labels, singleton_cluster_to_noise_points ami = adjusted_mutual_info_score(true_labels, cluster_labels) else: print("Noise points are ignored in the evaluation") - non_noise_points = (cluster_labels >= 0) - pct_clustered = (np.sum(non_noise_points) / cluster_labels.shape[0]) - print(f"Cluster coverage (%): {pct_clustered * 100:.2f}") - ari = adjusted_rand_score(true_labels[non_noise_points], - cluster_labels[non_noise_points]) - ami = adjusted_mutual_info_score(true_labels[non_noise_points], - cluster_labels[non_noise_points]) - # sil = silhouette_score(raw_data[non_noise_points], cluster_labels[non_noise_points]) + ari = adjusted_rand_score(true_labels[non_noise_points_mask], + cluster_labels[non_noise_points_mask]) + ami = adjusted_mutual_info_score(true_labels[non_noise_points_mask], + cluster_labels[non_noise_points_mask]) + # sil = silhouette_score(raw_data[non_noise_points_mask], cluster_labels[non_noise_points_mask]) else: + print(f"No noise points in the clustering result") ari = adjusted_rand_score(true_labels, cluster_labels) ami = adjusted_mutual_info_score(true_labels, cluster_labels) # sil = silhouette_score(raw_data, cluster_labels) - print(f"No noise points in the clustering result") print(f"ARI: {ari:.4f}") print(f"AMI: {ami:.4f}") diff --git a/script/benchmark/extract_clams_clustering_qualities.sh b/script/benchmark/extract_clams_clustering_qualities.sh new file mode 100755 index 0000000..e044d42 --- /dev/null +++ b/script/benchmark/extract_clams_clustering_qualities.sh @@ -0,0 +1,49 @@ +#!/bin/bash +# Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS +# Project Developers. See the top-level COPYRIGHT file for details. +# +# Extracts clustering parameters and results (ARI/AMI) from ClaMS benchmark +# log files (e.g., out.log). Values are tracked as a running state while +# scanning the log top to bottom, so parameter lines that appear less often +# than ARI/AMI lines are repeated/shared across all following ARI/AMI rows +# until the next occurrence of that parameter updates the state. +# +# Usage: extract_clustering_results.sh [ ...] + +set -euo pipefail + +if [[ $# -eq 0 ]]; then + echo "Usage: $0 [ ...]" >&2 + exit 1 +fi + +printf "kNNG_k,min_cluster_size,final_clusters,singleton_to_noise,cluster_coverage,ARI,AMI,file\n" + +for file in "$@"; do + awk -v fname="$file" ' + # Return the last whitespace/colon separated token on the line, i.e. its value + function lastfield(s, n, arr) { + n = split(s, arr, /[[:space:]:]+/) + return arr[n] + } + + BEGIN { + k = "NA"; mcs = "NA"; fc = "NA"; singleton = "NA"; coverage = "NA"; have_ari = 0; ari = "NA" + } + + /kNNG k:/ || /^[[:space:]]*k:[[:space:]]*[0-9]+[[:space:]]*$/ { k = lastfield($0); next } + /Min cluster size/ { mcs = lastfield($0); next } + /Evaluating Clustering Results, no noise points =/ { singleton = lastfield($0); next } + /#of final clusters/ { fc = lastfield($0); next } + /Cluster coverage \(%\)/ { coverage = lastfield($0); next } + /ARI/ { ari = lastfield($0); have_ari = 1; next } + /AMI/ { + if (have_ari) { + ami = lastfield($0) + printf "%s,%s,%s,%s,%s,%s,%s,%s\n", k, mcs, fc, singleton, coverage, ari, ami, fname + have_ari = 0 + } + next + } + ' "$file" +done diff --git a/script/benchmark/extract_clams_clustering_time.sh b/script/benchmark/extract_clams_clustering_time.sh new file mode 100755 index 0000000..6cc017f --- /dev/null +++ b/script/benchmark/extract_clams_clustering_time.sh @@ -0,0 +1,79 @@ +#!/bin/bash +# Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS +# Project Developers. See the top-level COPYRIGHT file for details. +# +# Extracts clustering execution times from ClaMS benchmark log files +# (e.g., out.log). Each stage is delimited by a ": " line +# and a matching "Finished : " line; the elapsed time +# between the two is reported in seconds. Only the first occurrence of each +# stage pair is used (e.g., "Running CLAMS-HDBSCAN" and +# "Assign clusters to noise points" can repeat per min-cluster-size value). +# +# Usage: extract_clams_clustering_time.sh [ ...] + +set -euo pipefail + +if [[ $# -eq 0 ]]; then + echo "Usage: $0 [ ...]" >&2 + exit 1 +fi + +printf "kNNG k,nodes,tasks/node,kNNG (s),MFC (s),AMST (s),CLAMS-HDBSCAN (s),Assigning noise points (s),file\n" + +for file in "$@"; do + awk -v fname="$file" ' + # Convert a "YYYY/MM/DD HH:MM:SS" timestamp found at the end of the line to epoch seconds + function to_epoch(line, ts, n, a) { + if (!match(line, /[0-9]{4}\/[0-9]{2}\/[0-9]{2} [0-9]{2}:[0-9]{2}:[0-9]{2}$/)) return -1 + ts = substr(line, RSTART, RLENGTH) + n = split(ts, a, /[\/ :]+/) + return mktime(a[1] " " a[2] " " a[3] " " a[4] " " a[5] " " a[6]) + } + + # Return the last whitespace/colon separated token on the line, i.e. its value + function lastfield(s, n, arr) { + n = split(s, arr, /[[:space:]:]+/) + return arr[n] + } + + function elapsed(start, finish) { + if (start == -1 || finish == -1) return "NA" + return finish - start + } + + BEGIN { + k = "NA"; nodes = "NA"; tasks = "NA" + knng_start = -1; knng_end = -1 + mfc_start = -1; mfc_end = -1 + amst_start = -1; amst_end = -1 + hdbscan_start = -1; hdbscan_end = -1 + noise_start = -1; noise_end = -1 + } + + /^kNNG k:/ && k == "NA" { k = lastfield($0) } + /^Compute nodes:/ && nodes == "NA" { nodes = lastfield($0) } + /^Tasks per node:/ && tasks == "NA" { tasks = lastfield($0) } + + /^Building KNNG:/ && knng_start == -1 { knng_start = to_epoch($0) } + /^Finished Building KNNG:/ && knng_start != -1 && knng_end == -1 { knng_end = to_epoch($0) } + + /^Connecting the CCs using MFC:/ && mfc_start == -1 { mfc_start = to_epoch($0) } + /^Finished Connecting the CCs using MFC:/ && mfc_start != -1 && mfc_end == -1 { mfc_end = to_epoch($0) } + + /^Running AMST, approx bound/ && amst_start == -1 { amst_start = to_epoch($0) } + /^Finished Running AMST, approx bound/ && amst_start != -1 && amst_end == -1 { amst_end = to_epoch($0) } + + /^Running CLAMS-HDBSCAN:/ && hdbscan_start == -1 { hdbscan_start = to_epoch($0) } + /^Finished Running CLAMS-HDBSCAN:/ && hdbscan_start != -1 && hdbscan_end == -1 { hdbscan_end = to_epoch($0) } + + /^Assign clusters to noise points:/ && noise_start == -1 { noise_start = to_epoch($0) } + /^Finished Assign clusters to noise points:/ && noise_start != -1 && noise_end == -1 { noise_end = to_epoch($0) } + + END { + printf "%s,%s,%s,%s,%s,%s,%s,%s,%s\n", k, nodes, tasks, \ + elapsed(knng_start, knng_end), elapsed(mfc_start, mfc_end), \ + elapsed(amst_start, amst_end), elapsed(hdbscan_start, hdbscan_end), \ + elapsed(noise_start, noise_end), fname + } + ' "$file" +done diff --git a/script/benchmark/extract_hpc_clustering_bench_result.sh b/script/benchmark/extract_hpc_clustering_bench_result.sh deleted file mode 100644 index b9a0f1c..0000000 --- a/script/benchmark/extract_hpc_clustering_bench_result.sh +++ /dev/null @@ -1,19 +0,0 @@ -JOBS=(job_20240925_085603 job_20240925_085709) - -# Execute a command and print the execution result without the newline -function exe() { - ret=$("$@") - echo -n $ret -} - -printf "Min Cluster Size, Clustered %%, ARI, AMI\n" -for job in "${JOBS[@]}"; do - echo "$job" - for i in {0..4}; do - echo -n $(cat bench_outputs/${job}/job_${i}/out.log | grep "Minimum cluster size: " | awk '{print $7}'); printf ",\t" - echo -n $(cat bench_outputs/${job}/job_${i}/out.log | grep % | awk '{print $4}'); printf ",\t" - echo -n $(cat bench_outputs/${job}/job_${i}/out.log | grep ARI | awk '{print $2}'); printf ",\t" - echo -n $(cat bench_outputs/${job}/job_${i}/out.log | grep AMI | awk '{print $2}') - echo "" - done -done \ No newline at end of file diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index fd9c0bd..dd0b964 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -56,6 +56,8 @@ def parse_options(): help='Use NEO-DNND (build_knng_neo) instead of build_knng.') parser.add_argument('--neodnnd_threads', type=int, default=2, help='Number of threads to use for NEO-DNND. Ignored if --neodnnd is not specified.') + parser.add_argument('--neodnnd_replicate_rate', type=float, default=0.0, + help='The replicate rate for NEO-DNND. Ignored if --neodnnd is not specified. See NEO-DNND\'s documentation for details.') parser.add_argument('--nng_r', type=float, default=0.5, help='The r (sampling) parameter for KNNG construction.') parser.add_argument('--nng_delta', type=float, default=0.0001, @@ -159,12 +161,19 @@ def generate_job_name(): time.sleep(2) return f"job_{datetime.now().strftime('%Y%m%d_%H%M%S')}" +stage_name_stack = [] # Stack to store stage names for logging def echo_stage_name(job_script, stage_name): job_script.write("echo\n") - job_script.write("date\n") job_script.write(f"echo ================================\n") - job_script.write(f"echo \"{stage_name}\"\n") + job_script.write(f"echo \"{stage_name}\": $(date \"+%Y/%m/%d %H:%M:%S\")\n") job_script.write(f"echo ================================\n") + stage_name_stack.append(stage_name) # Add the stage name to the queue + +def finish_stage(job_script): + if stage_name_stack: + stage_name = stage_name_stack.pop() # Get the last stage name from the stack + job_script.write(f"echo \"Finished {stage_name}: $(date \"+%Y/%m/%d %H:%M:%S\")\"\n") + job_script.write("echo\n") def add_clustering_evaluation(job_script, cluster_label_file, amst_ds_path, ground_truth_path, evaluator, @@ -218,13 +227,15 @@ def run_clustering(job_script, set_cmd, work_dir, amst_approx_bound, f" -c {cluster_tree_file} " f" -P ") add_cmd(hpc_clustering_command, job_script) + finish_stage(job_script) if ground_truth_path: - echo_stage_name(job_script, "Evaluating Clustering Results") + echo_stage_name(job_script, "Evaluating Clustering Results, no noise points = False") add_clustering_evaluation( job_script, cluster_label_file, amst_ds_path, ground_truth_path, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, singleton_cluster_to_noise_points) + finish_stage(job_script) echo_stage_name(job_script, "Assign clusters to noise points") cluster_label_file_no_noise = f"{cluster_label_file[:-4]}.noise_assigned.txt" @@ -233,14 +244,16 @@ def run_clustering(job_script, set_cmd, work_dir, amst_approx_bound, f"-c {cluster_label_file} " f"-o {cluster_label_file_no_noise}") add_cmd(cluster_assign_command, job_script) + finish_stage(job_script) echo_stage_name( job_script, - "Evaluate clustering results after assigning clusters to noise points") + "Evaluate clustering results, no noise points = True") add_clustering_evaluation( job_script, cluster_label_file_no_noise, amst_ds_path, ground_truth_path, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, singleton_cluster_to_noise_points) + finish_stage(job_script) job_script.write("echo \"\" \n") @@ -276,12 +289,13 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, # Set up the batch script header set_up_batch_header(job_script, job_name, job_dir, num_nodes) + job_script.write(f"echo \"Compute nodes: {num_nodes}\"\n") + job_script.write(f"echo \"Tasks per node: {num_tasks_per_node}\"\n") add_cmd(f'mkdir -p {work_dir}', job_script, True, False) # Run the DNND step - job_script.write("echo\n") - job_script.write("date\n") + job_script.write(f"echo \"kNNG k: {nng_k}\"\n") if len(input_dnnd_ds_path) == 0: echo_stage_name(job_script, "Building KNNG") dnnd_ds_path = f"{work_dir}/dnnd_pm_datastore" @@ -289,6 +303,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, verbose_flag = '-v' if verbose else '' dnnd_command = f"{dnnd_exe} {verbose_flag} -k {nng_k} -r {nng_r} -d {nng_delta} -f {distance_func} -o {dnnd_ds_path} -b {dnnd_batch_size} -p {points_file_format} {point_path}" add_srun_cmd(num_tasks_per_node, dnnd_command, job_script) + finish_stage(job_script) if backup_knng: dnnd_ds_path_backup = f"{dnnd_ds_path}_backup" job_script.write(f"echo \"Backing up KNNG datastore\"\n") @@ -303,6 +318,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, echo_stage_name(job_script, "Connecting the CCs using MFC") mfc_command = f"{mfc_exe} -d {dnnd_ds_path} -f {distance_func}" add_srun_cmd(num_tasks_per_node, mfc_command, job_script) + finish_stage(job_script) # Convert to core distance # TODO: Implement @@ -312,6 +328,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, add_cmd(f'mkdir -p {knng_coredist_dir}', job_script) conv2coredist_cmd = f"./src/conv_knng_to_core_dist -i {dnnd_ds_path} -o {knng_coredist_dir}/knng.txt -m {min_samples}" add_cmd(conv2coredist_cmd, job_script) + finish_stage(job_script) try_no = 0 for amst_approx_bound in amst_approx_bound_list: @@ -320,6 +337,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, amst_ds_path = f"{work_dir}/amst_pm_datastore_a{amst_approx_bound}" amst_command = f"{amst_exe} -d {dnnd_ds_path} -p {amst_ds_path} -e {amst_approx_bound}" add_srun_cmd(num_tasks_per_node, amst_command, job_script) + finish_stage(job_script) # Run the HPC Clustering step for set_cmd in min_cluster_size_set_cmnds: @@ -331,6 +349,8 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, noise_point_assigner_exe) try_no += 1 + finish_stage(job_script) + # If the file was not created, return an error if not os.path.exists(job_script_path): print(f"Error: Could not create the shell script {job_script_path}") @@ -369,7 +389,7 @@ def main(): if dnnd_exe == default_dnnd_exe: dnnd_exe = f'{os.getcwd()}/src/knng/build_knng_neo' # This is not the best way to set the number of threads for NEO-DNND, but it is a simple way to do it for now. - dnnd_exe = f'{dnnd_exe} -T {opts.neodnnd_threads}' + dnnd_exe = f'{dnnd_exe} -T {opts.neodnnd_threads} -R {opts.neodnnd_replicate_rate}' # Select HDBSCAN executable # If --distributed_hdbscan is sepecified and the user did not override --clustering_exe, diff --git a/src/clustering/cluster_noise_points.cpp b/src/clustering/cluster_noise_points.cpp index bf7c0a8..316268a 100644 --- a/src/clustering/cluster_noise_points.cpp +++ b/src/clustering/cluster_noise_points.cpp @@ -128,6 +128,10 @@ int main(int argc, char* argv[]) { spdlog::info( "Assigning cluster IDs to noise points by traversing the MST edges"); + // Start a timer to measure the time taken for assigning cluster IDs to noise + // points + auto kernel_timer = spdlog::stopwatch(); + std::size_t n_noise_points = 0; std::size_t n_assigned_points = 0; // NOTE: this algorithm is not determinstic because threads update the shared @@ -172,7 +176,9 @@ int main(int argc, char* argv[]) { spdlog::warn("Point {} could not be assigned to any cluster.", point_id); } } - spdlog::info("Finished assigning cluster IDs to noise points"); + const auto kernel_elapsed_time = kernel_timer.elapsed(); + spdlog::info("Finished assigning cluster IDs to noise points {}s", + kernel_elapsed_time.count()); spdlog::info("Number of noise points in the original data: {}", n_noise_points); spdlog::info("Number of remaining noise points: {}", From b0b78d9f6251c3d5c4d865e4269d6201cda5a21b Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 4 Sep 2026 08:07:12 -0700 Subject: [PATCH 20/37] Brush up benchmark scripts --- script/benchmark/extract_clams_clustering_qualities.sh | 8 ++++---- script/benchmark/run_clams_bench.py | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/script/benchmark/extract_clams_clustering_qualities.sh b/script/benchmark/extract_clams_clustering_qualities.sh index e044d42..74ff9d9 100755 --- a/script/benchmark/extract_clams_clustering_qualities.sh +++ b/script/benchmark/extract_clams_clustering_qualities.sh @@ -17,7 +17,7 @@ if [[ $# -eq 0 ]]; then exit 1 fi -printf "kNNG_k,min_cluster_size,final_clusters,singleton_to_noise,cluster_coverage,ARI,AMI,file\n" +printf "kNNG k,Min Cluster Size,Found Clusters,MST clustering,Cluster coverage,ARI,AMI,file\n" for file in "$@"; do awk -v fname="$file" ' @@ -28,19 +28,19 @@ for file in "$@"; do } BEGIN { - k = "NA"; mcs = "NA"; fc = "NA"; singleton = "NA"; coverage = "NA"; have_ari = 0; ari = "NA" + k = "NA"; mcs = "NA"; fc = "NA"; mst_clustering = "NA"; coverage = "NA"; have_ari = 0; ari = "NA" } /kNNG k:/ || /^[[:space:]]*k:[[:space:]]*[0-9]+[[:space:]]*$/ { k = lastfield($0); next } /Min cluster size/ { mcs = lastfield($0); next } - /Evaluating Clustering Results, no noise points =/ { singleton = lastfield($0); next } + /MST-based cluster guess/ { mst_clustering = lastfield($0); next } /#of final clusters/ { fc = lastfield($0); next } /Cluster coverage \(%\)/ { coverage = lastfield($0); next } /ARI/ { ari = lastfield($0); have_ari = 1; next } /AMI/ { if (have_ari) { ami = lastfield($0) - printf "%s,%s,%s,%s,%s,%s,%s,%s\n", k, mcs, fc, singleton, coverage, ari, ami, fname + printf "%s,%s,%s,%s,%s,%s,%s,%s\n", k, mcs, fc, mst_clustering, coverage, ari, ami, fname have_ari = 0 } next diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index dd0b964..66696fd 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -230,7 +230,8 @@ def run_clustering(job_script, set_cmd, work_dir, amst_approx_bound, finish_stage(job_script) if ground_truth_path: - echo_stage_name(job_script, "Evaluating Clustering Results, no noise points = False") + echo_stage_name(job_script, "Evaluating Clustering Results") + job_script.write(f"echo \"MST-based cluster guess: False\"\n") add_clustering_evaluation( job_script, cluster_label_file, amst_ds_path, ground_truth_path, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, @@ -246,9 +247,8 @@ def run_clustering(job_script, set_cmd, work_dir, amst_approx_bound, add_cmd(cluster_assign_command, job_script) finish_stage(job_script) - echo_stage_name( - job_script, - "Evaluate clustering results, no noise points = True") + echo_stage_name(job_script, "Evaluate clustering results") + job_script.write(f"echo \"MST-based cluster guess: True\"\n") add_clustering_evaluation( job_script, cluster_label_file_no_noise, amst_ds_path, ground_truth_path, evaluator, ygm_cluster_eval, From 17ec57456960b6bb76763c91f68b070824e25c1e Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Sat, 5 Sep 2026 04:52:56 -0700 Subject: [PATCH 21/37] Update output header in extract_clams_clustering_time.sh --- script/benchmark/extract_clams_clustering_time.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/script/benchmark/extract_clams_clustering_time.sh b/script/benchmark/extract_clams_clustering_time.sh index 6cc017f..6e4d828 100755 --- a/script/benchmark/extract_clams_clustering_time.sh +++ b/script/benchmark/extract_clams_clustering_time.sh @@ -18,7 +18,7 @@ if [[ $# -eq 0 ]]; then exit 1 fi -printf "kNNG k,nodes,tasks/node,kNNG (s),MFC (s),AMST (s),CLAMS-HDBSCAN (s),Assigning noise points (s),file\n" +printf "kNNG k,nodes,tasks/node,kNNG (s),MFC (s),AMST (s),CLAMS-HDBSCAN (s),Noise clustering (s),file\n" for file in "$@"; do awk -v fname="$file" ' From 3ab6ad5aba2d12bdc6dc5f42e569cf7d8117e92d Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Sun, 6 Sep 2026 06:00:42 -0700 Subject: [PATCH 22/37] Bugfix: fix to detect noise points correctly. --- script/benchmark/clustering_utilities.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/script/benchmark/clustering_utilities.py b/script/benchmark/clustering_utilities.py index 6cb80ae..81646e9 100644 --- a/script/benchmark/clustering_utilities.py +++ b/script/benchmark/clustering_utilities.py @@ -30,10 +30,14 @@ def eval_clusters(cluster_labels, true_labels, singleton_cluster_to_noise_points print(f"After filtering: {len(true_labels)} ground truth points") non_noise_points_mask = (cluster_labels >= 0) - pct_clustered = (np.sum(non_noise_points_mask) / cluster_labels.shape[0]) - print(f"Cluster coverage (%): {pct_clustered * 100:.2f}") + num_non_noise_points = np.sum(non_noise_points_mask) + pct_clustered = (num_non_noise_points / cluster_labels.shape[0]) + print(f"Cluster coverage (%): {pct_clustered * 100:.2f}%") - if len(non_noise_points_mask) < cluster_labels.shape[0]: # Has noise points + # Make sure always num_non_noise_points <= cluster_labels.shape[0] + assert num_non_noise_points <= cluster_labels.shape[0] + + if num_non_noise_points != cluster_labels.shape[0]: # Has noise points if singleton_cluster_to_noise_points: print( "Assigning a singleton cluster to each noise point in the clustering result") @@ -99,7 +103,7 @@ def find_files_in_dir(dir_path, ext=''): # If the first column is not a point ID, set has_ids to False. # If there are multiple files in a directory, the point IDs must be present. def read_point_data(data_path, has_ids=True): - print(f"Loading data from {data_path}") + print(f"Loading data from {data_path}", flush=True) files = find_files_in_dir(data_path) if len(files) == 0: @@ -137,7 +141,7 @@ def read_point_data(data_path, has_ids=True): points_table[pid] = list(map(float, items)) - print(f"Loaded {len(points_table)} items from {len(files)} files") + print(f"Loaded {len(points_table)} items from {len(files)} files", flush=True) # numpy array of feature vectors # if the IDs are not continuous, fill the missing IDs with -1 @@ -167,7 +171,7 @@ def read_point_data(data_path, has_ids=True): # # Both File types can also contain comment lines, which must start from #. def read_label_data(data_path): - print(f"Loading data from {data_path}") + print(f"Loading data from {data_path}", flush=True) labels_dict = {} files = [] if os.path.isdir(data_path): @@ -186,9 +190,9 @@ def read_label_data(data_path): break if contains_ids: - print("Loading point IDs and labels") + print("Loading point IDs and labels", flush=True) else: - print("Loading only labels") + print("Loading only labels", flush=True) if len(files) > 1 and not contains_ids: print("Multiple files are provided," @@ -220,8 +224,8 @@ def read_label_data(data_path): exit(1) labels_dict[pid] = label - print(f"Loaded {len(labels_dict)} items from {len(files)} files") - print(f"Max ID: {max_id}") + print(f"Loaded {len(labels_dict)} items from {len(files)} files", flush=True) + print(f"Max ID: {max_id}", flush=True) # numpy array of labels # if the IDs are not continuous, fill the missing IDs with -1 From 1ca4060a45916679a03acd789d064e9e14db61d0 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Sun, 6 Sep 2026 06:01:08 -0700 Subject: [PATCH 23/37] Brush up on the HDBSCAN script --- README.md | 4 +- script/benchmark/hdbscan/run_hdbscan.py | 188 ------------------ .../hdbscan/test_hdbscan_clustering.sh | 35 ---- script/benchmark/run_hdbscan_bench.py | 92 ++++++++- 4 files changed, 89 insertions(+), 230 deletions(-) delete mode 100644 script/benchmark/hdbscan/run_hdbscan.py delete mode 100755 script/benchmark/hdbscan/test_hdbscan_clustering.sh diff --git a/README.md b/README.md index e088b5f..ce87bd6 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,7 @@ general approach taken is to take the standard HDBSCAN algorithm and swap compon non-Euclidean data for scalable primitives. This often requires resorting to algorithms that are approximations of what is done in HDBSCAN, sometimes without approximation guarantees. -Note: This code works best if the input data is de-duplicated, with exact duplicate points removed. +Note: This code works best if the input data is de-duplicated, with exact duplicate points removed. If not, HDBSCAN will identify clusters with duplicate points as highly stable and always select them. @@ -84,7 +84,7 @@ source ./venv/bin/activate # In clams/build # -m: min cluster size # -s: min samples -python3 ../script/benchmark/hdbscan/run_hdbscan.py -m 10 -s 5 -p ../dataset/fashion-mnist/points.txt -g ../dataset/fashion-mnist/labels.txt +python3 ../script/benchmark/hdbscan_benchmark/run_hdbscan.py -m 10 -s 5 -p ../dataset/fashion-mnist/points.txt -g ../dataset/fashion-mnist/labels.txt ``` # License diff --git a/script/benchmark/hdbscan/run_hdbscan.py b/script/benchmark/hdbscan/run_hdbscan.py deleted file mode 100644 index e7a5ace..0000000 --- a/script/benchmark/hdbscan/run_hdbscan.py +++ /dev/null @@ -1,188 +0,0 @@ -# Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS -# Project Developers. See the top-level COPYRIGHT file for details. - - -# Description: This script generates a synthetic dataset and runs HDBSCAN on it. -# It also can read a point file and run HDBSCAN on it. -# Usage: -# python run_hdbscan.py -n 1000 -m 5 -# Or -# python run_hdbscan.py -p points.txt -m 5 - - -from sklearn.datasets import make_blobs -import hdbscan -import os -import argparse -import numpy as np -from clustering_utilities import * - - -def parse_options(): - parser = argparse.ArgumentParser( - description='Evaluate kNN index') - - # For input point (feature) data - parser.add_argument('-p', '--point_data_path', - dest='point_data_path', - required=False, action='store', type=str, - help='Input point file path') - - parser.add_argument('-I', '--has_ids', - dest='has_ids', - required=False, action='store_true', - help='If specified, the input point file has point IDs') - - # For generate synthetic data - parser.add_argument('-n', '--n_samples', - dest='n_samples', - required=False, action='store', type=int, - default=1000, - help='#of samples to generate by the data generator') - parser.add_argument('-f', '--n_features', - dest='n_features', - required=False, action='store', type=int, - default=32, - help='#of features to generate by the data generator') - - # HDBSCAN parameters - parser.add_argument('-m', '--min_cluster_size', - dest='min_cluster_size', - required=False, action='store', type=int, - default=5, - help='Minimum cluster size.') - parser.add_argument('-s', '--min_samples', - dest='min_samples', - required=False, action='store', type=int, - default=None, - help='#of samples in a neighborhood for a point to be considered as a core point.' - 'When None, defaults to min_cluster_size') - - # Output file paths - parser.add_argument('-o', '--cluster_labels_out_path', - dest='cluster_labels_out_path', - required=False, action='store', type=str, - default='out_clusters.txt', - help='File path to store computed cluster IDs') - parser.add_argument('-M', '--mst_out_path', - dest='mst_out_path', - required=False, action='store', type=str, - help='Output file path for intermediate MST data') - parser.add_argument('-c', '--condensed_tree_out_path', - dest='condensed_tree_out_path', - required=False, action='store', type=str, - default=None, - help='File path to store internal cluster data (mainly for debugging)') - parser.add_argument('-C', '--cluster_persistence_out_path', - dest='cluster_persistence_out_path', - required=False, action='store', type=str, - default=None, - help='File path to store cluster persistence data') - - # Ground truth file path - parser.add_argument('-g', '--ground_truth', - dest='gt_file', - required=False, action='store', type=str, - help='Ground truth file path') - - args = parser.parse_args() - return args - - -def run_hdbscan_kernel(points, min_cluster_size, min_samples, - gen_min_span_tree=False): - clusters = hdbscan.HDBSCAN(gen_min_span_tree=gen_min_span_tree, - min_cluster_size=min_cluster_size, - min_samples=min_samples, - core_dist_n_jobs=-1) - print(f'\nStart clustering') - show_time_now() - clusters.fit(points) - print(f'Finish clustering') - show_time_now() - - return clusters - - -def run_hdbscan(points, min_cluster_size, min_samples, - cluster_labels_out_path='out_clusters.txt', - gt_file=None, - mst_out_path=None, - condensed_tree_out_path=None, - cluster_persistence_out_path=None, - assign_cluster_to_noise=False): - clusters = run_hdbscan_kernel(points, min_cluster_size, min_samples, - gen_min_span_tree=(mst_out_path is not None)) - - # Evaluate the clustering quality - if gt_file: - print('\nLoading ground truth data') - gt_labels = read_label_data(gt_file) - - print('\nEvaluating clustering quality') - eval_clusters(clusters.labels_, gt_labels) - - if assign_cluster_to_noise: - print('\nAssigning a cluster ID to every noise point') - no_noise_labels = assign_singleton_cluster_to_noise_point( - clusters.labels_) - eval_clusters(no_noise_labels, gt_labels) - - # Save the condensed tree data - if condensed_tree_out_path: - print(f'\nSaving condensed tree data in {condensed_tree_out_path}') - clusters.condensed_tree_.to_pandas().to_csv(condensed_tree_out_path) - - # Save the MST data - if mst_out_path: - print(f'\nSaving MST data in {mst_out_path}') - with open(mst_out_path, 'w') as fout_mst: - mst = clusters.minimum_spanning_tree_.to_numpy() - # Format: Point0, Point1, Distance - for edge in mst: - fout_mst.write(f'{int(edge[0])}\t{int(edge[1])}\t{edge[2]}\n') - fout_mst.close() - show_time_now() - - # Save the cluster IDs. - if cluster_labels_out_path: - print(f'\nSaving cluster IDs in {cluster_labels_out_path}') - with open(cluster_labels_out_path, 'w') as fout: - fout.write(f'# Node ID\tCluster ID\n') - for i, label in enumerate(clusters.labels_): - fout.write(f'{i}\t{label}\n') - show_time_now() - print(f'Cluster IDs are saved in {cluster_labels_out_path}') - - # Save the cluster persistence data - if cluster_persistence_out_path: - print( - f'\nSaving cluster persistence data in {cluster_persistence_out_path}') - with open(cluster_persistence_out_path, 'w') as fout: - fout.write('Cluster ID\tPersistence\n') - for i, persistence in enumerate(clusters.cluster_persistence_): - fout.write(f'{i}\t{persistence}\n') - - -def main(): - opt = parse_options() - - show_time_now() - if opt.point_data_path: - points = read_point_data(opt.point_data_path, opt.has_ids) - else: - points, _ = make_blobs(n_samples=opt.n_samples, - n_features=opt.n_features, - centers=10) - print(f'points data shape: {points.shape}') - - run_hdbscan(points, opt.min_cluster_size, opt.min_samples, - gt_file=opt.gt_file, - cluster_labels_out_path=opt.cluster_labels_out_path, - mst_out_path=opt.mst_out_path, - condensed_tree_out_path=opt.condensed_tree_out_path, - cluster_persistence_out_path=opt.cluster_persistence_out_path) - - -if __name__ == '__main__': - main() diff --git a/script/benchmark/hdbscan/test_hdbscan_clustering.sh b/script/benchmark/hdbscan/test_hdbscan_clustering.sh deleted file mode 100755 index 8efab4c..0000000 --- a/script/benchmark/hdbscan/test_hdbscan_clustering.sh +++ /dev/null @@ -1,35 +0,0 @@ -#!/bin/bash - -# Usage: -# cd /build -# ../script/test_hdbscan_clustering.sh - -mst_file="./blobs_mst.txt" -hdbscan_cluster_ids_file="./hdbscan_cluster_ids_file.txt" -num_samples=$((2**15)) -min_cluster_size=5 - -echo "Running HDBSCAN" -python3 ../script/run_hdbscan.py --n_samples $num_samples --min_cluster_size $min_cluster_size --mst_out_path ${mst_file} --cluster_labels_out_path ${hdbscan_cluster_ids_file} - - -cluster_ids_file="out_cluster_ids.txt" -echo "" -./src/run_hdbscan_clustering -i ${mst_file} -o ${cluster_ids_file} -m $min_cluster_size - -echo "" -echo "Reassigning cluster ids" -python3 ../script/reassign_sequential_cluster_ids.py -i ${hdbscan_cluster_ids_file} -o ${hdbscan_cluster_ids_file} - -echo "" -echo "Reassigning cluster ids" -python3 ../script/reassign_sequential_cluster_ids.py -i ${cluster_ids_file} -o ${cluster_ids_file} - -echo "" -echo "Running diff" -diff ${hdbscan_cluster_ids_file} ${cluster_ids_file} > diff.txt - -echo "#of lines in diff.txt:" -wc -l diff.txt - -echo "Done" \ No newline at end of file diff --git a/script/benchmark/run_hdbscan_bench.py b/script/benchmark/run_hdbscan_bench.py index 88b40d0..94b8a6c 100644 --- a/script/benchmark/run_hdbscan_bench.py +++ b/script/benchmark/run_hdbscan_bench.py @@ -11,9 +11,13 @@ """ import argparse +import hdbscan +import os +import numpy as np +import time + from clustering_utilities import * -from script.benchmark.bench_utilities import * -from script.benchmark.hdbscan.run_hdbscan import run_hdbscan +from bench_utilities import * def parse_options(): @@ -67,9 +71,87 @@ def parse_options(): return args +def run_hdbscan_kernel(points, min_cluster_size, min_samples, + gen_min_span_tree=False): + clusters = hdbscan.HDBSCAN(gen_min_span_tree=gen_min_span_tree, + min_cluster_size=min_cluster_size, + min_samples=min_samples, + core_dist_n_jobs=-1) + print(f'\nStart clustering', flush=True) + show_time_now() + # start timer + start_time = time.time() + clusters.fit(points) + print(f'Finish clustering in {time.time() - start_time:.2f} seconds', flush=True) + show_time_now() + + return clusters + + +def run_hdbscan(points, min_cluster_size, min_samples, + cluster_labels_out_path='out_clusters.txt', + gt_file=None, + mst_out_path=None, + condensed_tree_out_path=None, + cluster_persistence_out_path=None, + assign_cluster_to_noise=False): + clusters = run_hdbscan_kernel(points, min_cluster_size, min_samples, + gen_min_span_tree=(mst_out_path is not None)) + + # Evaluate the clustering quality + if gt_file: + print('\nLoading ground truth data') + gt_labels = read_label_data(gt_file) + + print('\nEvaluating clustering quality') + eval_clusters(clusters.labels_, gt_labels) + + if assign_cluster_to_noise: + print('\nAssigning a cluster ID to every noise point') + no_noise_labels = assign_singleton_cluster_to_noise_point( + clusters.labels_) + eval_clusters(no_noise_labels, gt_labels) + + # Save the condensed tree data + if condensed_tree_out_path: + print(f'\nSaving condensed tree data in {condensed_tree_out_path}') + clusters.condensed_tree_.to_pandas().to_csv(condensed_tree_out_path) + + # Save the MST data + if mst_out_path: + print(f'\nSaving MST data in {mst_out_path}') + with open(mst_out_path, 'w') as fout_mst: + mst = clusters.minimum_spanning_tree_.to_numpy() + # Format: Point0, Point1, Distance + for edge in mst: + fout_mst.write(f'{int(edge[0])}\t{int(edge[1])}\t{edge[2]}\n') + fout_mst.close() + show_time_now() + + # Save the cluster IDs. + if cluster_labels_out_path: + print(f'\nSaving cluster IDs in {cluster_labels_out_path}') + with open(cluster_labels_out_path, 'w') as fout: + fout.write(f'# Node ID\tCluster ID\n') + for i, label in enumerate(clusters.labels_): + fout.write(f'{i}\t{label}\n') + show_time_now() + print(f'Cluster IDs are saved in {cluster_labels_out_path}') + + # Save the cluster persistence data + if cluster_persistence_out_path: + print( + f'\nSaving cluster persistence data in {cluster_persistence_out_path}') + with open(cluster_persistence_out_path, 'w') as fout: + fout.write('Cluster ID\tPersistence\n') + for i, persistence in enumerate(clusters.cluster_persistence_): + fout.write(f'{i}\t{persistence}\n') + + def main(): opts = parse_options() + points = read_point_data(opts.point_data_path, opts.has_ids) min_cluster_size_list = parse_range(opts.min_cluster_size_range) @@ -77,9 +159,9 @@ def main(): for min_cluster_size in min_cluster_size_list: for min_samples in min_samples_list: - print(f'\n--------------------------------') - print((f"min_cluster_size: {min_cluster_size}")) - print((f"min_samples: {min_samples}")) + print(f'\n--------------------------------', flush=True) + print((f"min_cluster_size: {min_cluster_size}"), flush=True) + print((f"min_samples: {min_samples}"), flush=True) run_hdbscan(points, min_cluster_size, min_samples, gt_file=opts.gt_file, cluster_labels_out_path=opts.cluster_labels_out_path, From 13db1e767a6ef6c04db51086796ca7b836d09a65 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 11 Sep 2026 11:53:10 -0700 Subject: [PATCH 24/37] Add knng and clusters correlation evaluation script --- .../evaluate_correlation_knng_and_clusters.py | 299 ++++++++++++++++++ 1 file changed, 299 insertions(+) create mode 100644 script/benchmark/evaluate_correlation_knng_and_clusters.py diff --git a/script/benchmark/evaluate_correlation_knng_and_clusters.py b/script/benchmark/evaluate_correlation_knng_and_clusters.py new file mode 100644 index 0000000..b41784b --- /dev/null +++ b/script/benchmark/evaluate_correlation_knng_and_clusters.py @@ -0,0 +1,299 @@ +#!/usr/bin/env python3 +""" +Evaluate the correlation between KNNG neighbors and clusters. + +Take a KNNG dump and a cluster label file, then count how many of the first k +neighbors of each point belong to the same cluster as the point itself. + +Usage: + python evaluate_correlation_knng_and_clusters.py -g knng_dir -c clusters.txt -k 10 + +Arguments: + -g/--graph: KNNG file or directory containing KNNG dump files. + -c/--cluster: cluster file or directory containing cluster-label files. + -k/--neighbors: number of neighbors to consider for each point. + +KNNG file format: + Each record is encoded as two lines: + ... + ... + The leading distance is a dummy value corresponding to the source ID and will be ignored. + +Cluster file format: + Each line contains a point and its cluster ID separated by whitespace. + The first line can be a comment starting with '#'. + A file containing only cluster IDs is also accepted; in that case point IDs are + assumed to be 0, 1, 2, ... in order. +""" + +from __future__ import annotations + +import argparse +import glob +from pathlib import Path +from typing import Iterable, Iterator, Sequence + + +def resolve_files(paths: Sequence[str]) -> list[Path]: + """Resolve file, directory, and glob inputs to a unique sorted file list.""" + resolved: set[Path] = set() + for value in paths: + path = Path(value) + if path.is_dir(): + matches = set(path.glob("*.txt")) | set(path.glob("knng-*.txt")) + elif path.is_file(): + matches = {path} + else: + matches = {Path(match) for match in glob.glob(value)} + + resolved.update(match.resolve() for match in matches if match.is_file()) + + if not resolved: + raise FileNotFoundError(f"No files matched: {paths}") + + return sorted(resolved) + + +def iter_knng_records(path: Path, contains_distance: bool) -> Iterator[tuple[int, list[int]]]: + """Yield (source_id, neighbor_ids) pairs from one KNNG dump file.""" + with path.open("r", encoding="utf-8") as knng_file: + while True: + ids_line = knng_file.readline() + if not ids_line: + return + + stripped = ids_line.strip() + if not stripped: + # This is a critical error. + raise ValueError(f"{path}: failed to parse row: {ids_line!r}") + + if contains_distance: + # Skip the distance line if it exists, as it will be processed later. + knng_file.readline() + + src_id = int(stripped.split()[0]) + neighbor_ids = [int(value) for value in stripped.split()[1:]] + + yield src_id, neighbor_ids + +def read_cluster_labels(cluster_path: str | Path) -> dict[int, int]: + """Read cluster labels from a file or directory into a point->cluster map.""" + labels: dict[int, int] = {} + files = resolve_files([str(cluster_path)]) + + if not files: + raise FileNotFoundError(f"No cluster files found in {cluster_path}") + + # Heuristic: determine whether the file format contains point IDs by scanning + # the first non-comment, non-empty line. + first_file = files[0] + contains_ids = False + with first_file.open("r") as handle: + for line in handle: + stripped = line.strip() + if not stripped or stripped.startswith("#"): + continue + if len(stripped.split()) >= 2: + contains_ids = True + break + + for file in files: + with file.open("r") as handle: + for line in handle: + stripped = line.strip() + if not stripped or stripped.startswith("#"): + continue + + items = stripped.split() + if contains_ids: + if len(items) < 2: + raise ValueError(f"{file}: invalid cluster row: {line.rstrip()}") + point_id, cluster_id = int(items[0]), int(items[1]) + else: + point_id = len(labels) + cluster_id = int(items[0]) + + if point_id in labels: + raise ValueError(f"Duplicate point ID {point_id} in cluster file {file}") + labels[point_id] = cluster_id + + if not labels: + raise ValueError(f"No cluster labels were loaded from {cluster_path}") + + return labels + + +def evaluate_knng_clusters(graph_path: str | Path, cluster_path: str | Path, k: int = 10, no_distance: bool = False): + """Return per-point and aggregate same-cluster neighbor statistics.""" + if k <= 0: + raise ValueError("-k/--neighbors must be positive") + + labels = read_cluster_labels(cluster_path) + print(f"Loaded {len(labels)} cluster labels from {cluster_path}") + + files = resolve_files([str(graph_path)]) + + per_point: list[tuple[int, int, int, int, float]] = [] + total_same = 0 + total_neighbors = 0 + covered_points = 0 + noise_points = 0 + + for knng_file in files: + print(f"Evaluating KNNG file {knng_file} with k={k}") + for source_id, neighbor_ids in iter_knng_records(knng_file, not no_distance): + if source_id not in labels: + continue + + source_label = labels[source_id] + if source_label == -1: + # Skip points that are not assigned to any cluster. + noise_points += 1 + continue + + limited_neighbors = neighbor_ids[:k] + if not limited_neighbors: + continue + + same_count = 0 + for neighbor_id in limited_neighbors: + if labels.get(neighbor_id) == source_label: + same_count += 1 + + ratio = same_count / len(limited_neighbors) + per_point.append((source_id, source_label, same_count, len(limited_neighbors), ratio)) + total_same += same_count + total_neighbors += len(limited_neighbors) + covered_points += 1 + + avg_ratio = (total_same / total_neighbors) if total_neighbors > 0 else 0.0 + return per_point, { + "num_points": covered_points, + "num_noise_points": noise_points, + "num_same_cluster_neighbors": total_same, + "num_considered_neighbors": total_neighbors, + "avg_same_cluster_fraction": avg_ratio, + "num_clustered_points": len(labels), + } + + +def _parse_k_values(raw_values: Sequence[Sequence[str]] | None) -> list[int]: + """Normalize one or more k values from CLI input (supports repeated or comma-separated lists).""" + values: list[int] = [] + raw_entries = raw_values or [] + + for group in raw_entries: + if isinstance(group, str): + group_tokens = [group] + else: + group_tokens = list(group) + + for token in group_tokens: + for item in str(token).split(","): + item = item.strip() + if not item: + continue + try: + value = int(item) + except ValueError as exc: + raise argparse.ArgumentTypeError(f"invalid neighbor count: {item!r}") from exc + if value <= 0: + raise argparse.ArgumentTypeError(f"neighbor count must be positive: {value}") + values.append(value) + + if not values: + return [10] + return values + + +def _parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description=( + "Evaluate how often a point's nearest neighbors belong to the same " + "cluster as the point itself." + ) + ) + parser.add_argument( + "-g", + "--graph", + required=True, + help="KNNG file or directory containing KNNG dump files", + ) + parser.add_argument( + "-N", + "--no-distance", + action="store_true", + help="Indicate if the KNNG files *do not* contain distance information", + ) + parser.add_argument( + "-c", + "--cluster", + required=True, + help="Cluster label file or directory containing cluster-label files", + ) + parser.add_argument( + "-k", + "--neighbors", + action="append", + nargs="+", + default=None, + help=( + "Neighbor count(s) to evaluate; accepts repeated uses or comma-separated " + "values, e.g. -k 5 10 or -k 3,5,10" + ), + ) + parser.add_argument( + "-o", + "--output", + type=str, + default=None, + help="Optional output file path for the summary table", + ) + args = parser.parse_args() + args.neighbors = _parse_k_values(args.neighbors or [[10]]) + return args + + +def main() -> None: + args = _parse_args() + + summary_lines: list[str] = [] + per_k_output: list[tuple[int, list[str]]] = [] + + for k in args.neighbors: + try: + per_point, summary = evaluate_knng_clusters(args.graph, args.cluster, k, no_distance=args.no_distance) + except Exception as exc: # pragma: no cover - CLI error path + raise SystemExit(f"Error: {exc}") from exc + + summary_lines.append(f"=== k={k} ===") + summary_lines.append( + f"Evaluated {summary['num_points']} source points; " + f"Skipped {summary['num_noise_points']} noise points" + ) + summary_lines.append( + f"Compared the first {k} neighbors per point; " + f"{summary['num_same_cluster_neighbors']}/{summary['num_considered_neighbors']} " + f"were in the same cluster (fraction={summary['avg_same_cluster_fraction']:.4f})" + ) + + per_k_lines = ["point_id cluster_id same_cluster total_checked fraction"] + for source_id, cluster_id, same_count, checked_count, ratio in per_point: + line = f"{source_id:>8} {cluster_id:>10} {same_count:>12} {checked_count:>12} {ratio:>8.4f}" + per_k_lines.append(line) + per_k_output.append((k, per_k_lines)) + + for line in summary_lines: + print(line) + + if args.output: + output_path = Path(args.output) + output_path.parent.mkdir(parents=True, exist_ok=True) + with output_path.open("w", encoding="utf-8") as handle: + handle.write("\n".join(line for _, lines in per_k_output for line in lines) + "\n") + print(f"\nPer-k statistics written to {output_path}") + + +if __name__ == "__main__": + main() + From 0164955e7006bd7b5f85f80463d391beacc72594 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 11 Sep 2026 11:54:27 -0700 Subject: [PATCH 25/37] Add MST-based clusetring complement v2 --- src/clustering/CMakeLists.txt | 5 + src/clustering/cluster_noise_points_v2.cpp | 203 +++++++++++++++++++++ src/details/multithread_adjacency_list.hpp | 61 +++---- 3 files changed, 230 insertions(+), 39 deletions(-) create mode 100644 src/clustering/cluster_noise_points_v2.cpp diff --git a/src/clustering/CMakeLists.txt b/src/clustering/CMakeLists.txt index 59a8ed2..de738e1 100644 --- a/src/clustering/CMakeLists.txt +++ b/src/clustering/CMakeLists.txt @@ -12,6 +12,11 @@ setup_omp_target(cluster_noise_points) setup_metall_target(cluster_noise_points) setup_spdlog_target(cluster_noise_points) +add_basic_executable(cluster_noise_points_v2 cluster_noise_points_v2.cpp) +setup_omp_target(cluster_noise_points_v2) +setup_metall_target(cluster_noise_points_v2) +setup_spdlog_target(cluster_noise_points_v2) + add_basic_executable(run_distributed_hdbscan_clustering run_distributed_hdbscan_clustering.cpp) setup_ygm_target(run_distributed_hdbscan_clustering) setup_omp_target(run_distributed_hdbscan_clustering) diff --git a/src/clustering/cluster_noise_points_v2.cpp b/src/clustering/cluster_noise_points_v2.cpp new file mode 100644 index 0000000..e2cbb8b --- /dev/null +++ b/src/clustering/cluster_noise_points_v2.cpp @@ -0,0 +1,203 @@ +// Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS +// Project Developers. See the top-level COPYRIGHT file for details. + +// Assign cluster IDs to noise points by traversing the MST edges. +// Traverse the MST edges from each noise point in BFS manner until a point that +// belongs to a cluster is found. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "../common.hpp" +#include "../details/multithread_adjacency_list.hpp" + +using namespace clams; + +template +using map_t = boost::unordered::unordered_flat_map; +using graph_t = multithread_adjacency_list>; + +struct option { + std::filesystem::path mst_edges_path; + bool metall_mst{false}; + std::filesystem::path cluster_ids_input_path; + std::filesystem::path cluster_ids_out_path; +}; + +void show_help() { + std::cout + << "<>\n" + "Required arguments:\n" + " -m path to a file or a directory that contains input MST.\n" + " -c path to input cluster IDs.\n" + " -o path to output cluster IDs.\n" + "Optional arguments:\n" + " -M If specified, input is a Metall datastore.\n" + " -h Show help." + << std::endl; +} + +void parse_option(int argc, char* argv[], option& opt) { + int opt_char; + while ((opt_char = getopt(argc, argv, "m:c:o:Mh")) != -1) { + switch (opt_char) { + case 'm': + opt.mst_edges_path = optarg; + break; + case 'c': + opt.cluster_ids_input_path = optarg; + break; + case 'o': + opt.cluster_ids_out_path = optarg; + break; + case 'M': + opt.metall_mst = true; + break; + case 'h': + show_help(); + std::exit(EXIT_SUCCESS); + default: + show_help(); + std::exit(EXIT_FAILURE); + } + } +} + +int main(int argc, char* argv[]) { + option opt; + parse_option(argc, argv, opt); + + graph_t mst_graph; + if (opt.metall_mst) { + spdlog::info("Attaching MST in Metall datastore"); + metall::manager metall_manager(metall::open_read_only, opt.mst_edges_path); + auto* input_mst_edges = + metall_manager.find(metall::unique_instance) + .first; + if (!input_mst_edges) { + spdlog::critical("Failed to find MST edges in Metall datastore at {}", + opt.mst_edges_path.string()); + std::abort(); + } + spdlog::info("#of MST edges: {}", input_mst_edges->size()); + spdlog::info("Constructing MST graph from Metall datastore"); + OMP_DIRECTIVE(parallel for) + for (size_t i = 0; i < input_mst_edges->size(); ++i) { + const auto& edge = input_mst_edges->at(i); + mst_graph.add(edge.ids[0], {edge.ids[1], edge.distance}); + mst_graph.add(edge.ids[1], {edge.ids[0], edge.distance}); + } + } else { + spdlog::info("Reading MST edges"); + weighted_edge_list_t input_mst_edges; + read_edges(opt.mst_edges_path, input_mst_edges); + spdlog::info("#of MST edges: {}", input_mst_edges.size()); + spdlog::info("Constructing MST graph from input file"); + OMP_DIRECTIVE(parallel for) + for (size_t i = 0; i < input_mst_edges.size(); ++i) { + const auto& edge = input_mst_edges.at(i); + mst_graph.add(edge.ids[0], {edge.ids[1], edge.distance}); + mst_graph.add(edge.ids[1], {edge.ids[0], edge.distance}); + } + } + spdlog::info("Finished constructing MST graph"); + + if (mst_graph.empty()) { + spdlog::warn("No MST edges found in the input file or directory: {}", + opt.mst_edges_path.string()); + return EXIT_SUCCESS; + } + + map_t point_cluster_map; + read_cluster_ids(opt.cluster_ids_input_path, point_cluster_map); + spdlog::info("Read {} points' cluster IDs from {}", point_cluster_map.size(), + opt.cluster_ids_input_path.string()); + + spdlog::info("Assigning cluster IDs to noise points by edge weight"); + auto kernel_timer = spdlog::stopwatch(); + + using heap_edge_t = std::tuple; + std::priority_queue, std::greater<>> + edge_heap; + + // Initialize the heap with edges connected to non-noise points + std::size_t n_noise_points = 0; + std::size_t n_assigned_points = 0; + for (const auto& [point_id, cluster_id] : point_cluster_map) { + if (cluster_id == k_noise_cluster_id) { + ++n_noise_points; + continue; + } + + for (auto edge_itr = mst_graph.values_begin(point_id); + edge_itr != mst_graph.values_end(point_id); ++edge_itr) { + const auto& neighbor_id = edge_itr->first; + const auto& distance = edge_itr->second; + if (point_cluster_map.at(neighbor_id) == k_noise_cluster_id) { + edge_heap.emplace(distance, point_id, neighbor_id); + } + } + } + + while (!edge_heap.empty()) { + const auto edge = edge_heap.top(); + edge_heap.pop(); + const auto predecessor_id = std::get<1>(edge); + const auto current_point_id = std::get<2>(edge); + + if (point_cluster_map.at(current_point_id) != k_noise_cluster_id) { + // Skip this edge if the neighbor point has already been visited. + continue; + } + + const auto source_cluster_id = point_cluster_map.at(predecessor_id); + point_cluster_map.at(current_point_id) = source_cluster_id; + ++n_assigned_points; + + for (auto edge_itr = mst_graph.values_begin(current_point_id); + edge_itr != mst_graph.values_end(current_point_id); ++edge_itr) { + const auto& neighbor_id = edge_itr->first; + const auto& distance = edge_itr->second; + if (point_cluster_map.at(neighbor_id) != k_noise_cluster_id) { + continue; + } + edge_heap.emplace(distance, current_point_id, neighbor_id); + } + } + + // for (const auto& [point_id, cluster_id] : point_cluster_map) { + // if (cluster_id == k_noise_cluster_id) { + // spdlog::warn("Point {} could not be assigned to any cluster.", + // point_id); + // } + // } + const auto kernel_elapsed_time = kernel_timer.elapsed(); + spdlog::info("Finished assigning cluster IDs to noise points {}s", + kernel_elapsed_time.count()); + spdlog::info("Number of noise points in the original data: {}", + n_noise_points); + spdlog::info("Number of remaining noise points: {}", + n_noise_points - n_assigned_points); + + dump_point_cluster_ids(point_cluster_map, opt.cluster_ids_out_path); + spdlog::info("Dumped point cluster IDs to {}", + opt.cluster_ids_out_path.string()); + + return EXIT_SUCCESS; +} diff --git a/src/details/multithread_adjacency_list.hpp b/src/details/multithread_adjacency_list.hpp index d6d36f4..3a6eb33 100644 --- a/src/details/multithread_adjacency_list.hpp +++ b/src/details/multithread_adjacency_list.hpp @@ -12,22 +12,15 @@ #pragma once -#include -#include #include +#include #include #include +#include -#define CLAMS_USE_STL_CONTAINERS_IN_ADJLIST 0 -#if CLAMS_USE_STL_CONTAINERS_IN_ADJLIST -#include -#include -#include -#else #include -#include +#include #include -#endif #include #include @@ -35,19 +28,15 @@ namespace clams { namespace { -#if CLAMS_USE_STL_CONTAINERS_IN_ADJLIST -namespace container = std; -#else namespace container = metall::container; -#endif -} +} // namespace template > class multithread_adjacency_list { public: - using key_type = _key_type; - using value_type = _value_type; + using key_type = _key_type; + using value_type = _value_type; static constexpr std::size_t k_num_banks = 1024; private: @@ -60,10 +49,9 @@ class multithread_adjacency_list { using key_table_allocator_type = container::scoped_allocator_adaptor< other_allocator_type>>; - using key_table_type = - container::unordered_map, - std::equal_to, - key_table_allocator_type>; + using key_table_type = container::unordered_node_map< + key_type, list_type, metall::utility::hash<>, std::equal_to, + key_table_allocator_type>; using bank_table_allocator_type = container::scoped_allocator_adaptor>; @@ -74,10 +62,10 @@ class multithread_adjacency_list { class impl_const_key_iterator; public: - using const_key_iterator = impl_const_key_iterator; - using value_iterator = typename list_type::iterator; - using const_value_iterator = typename list_type::const_iterator; - using local_key_iterator = typename key_table_type::iterator; + using const_key_iterator = impl_const_key_iterator; + using value_iterator = typename list_type::iterator; + using const_value_iterator = typename list_type::const_iterator; + using local_key_iterator = typename key_table_type::iterator; using const_local_key_iterator = typename key_table_type::const_iterator; explicit multithread_adjacency_list( @@ -90,14 +78,7 @@ class multithread_adjacency_list { auto guard = metall::utility::mutex::mutex_lock(bank_index(key)); #ifdef __clang__ -#if CLAMS_USE_STL_CONTAINERS_IN_ADJLIST - m_bank_table[bank_index(key)][key].emplace_back(std::move(value)); -#else m_bank_table[bank_index(key)][key].emplace_back(std::move(value)); - // m_bank_table[bank_index(key)].try_emplace(key, - // list_allocator_type(m_bank_table.get_allocator())); - // m_bank_table[bank_index(key)].at(key).emplace_back(std::move(value)); -#endif #else // MEMO: GCC does not work with STL Containers (tested with GCC 10.2.0 on // MacOS) @@ -106,6 +87,8 @@ class multithread_adjacency_list { return true; } + bool empty() const { return num_keys() == 0; } + std::size_t num_keys() const { std::size_t count = 0; for (auto &table : m_bank_table) { @@ -189,9 +172,9 @@ class multithread_adjacency_list< public: using difference_type = typename local_iterator_type::difference_type; - using value_type = typename local_iterator_type::value_type; - using pointer = typename local_iterator_type::pointer; - using reference = typename local_iterator_type::reference; + using value_type = typename local_iterator_type::value_type; + using pointer = typename local_iterator_type::pointer; + using reference = typename local_iterator_type::reference; explicit impl_const_key_iterator( const adjacency_list_type *const adjacency_list) @@ -210,10 +193,10 @@ class multithread_adjacency_list< } impl_const_key_iterator(const impl_const_key_iterator &) = default; - impl_const_key_iterator(impl_const_key_iterator &&) = default; + impl_const_key_iterator(impl_const_key_iterator &&) = default; impl_const_key_iterator &operator=(const impl_const_key_iterator &) = default; - impl_const_key_iterator &operator=(impl_const_key_iterator &&) = default; + impl_const_key_iterator &operator=(impl_const_key_iterator &&) = default; void move_to_end() { m_current_bank_index = m_ptr_adjacency_list->m_bank_table.size(); @@ -281,8 +264,8 @@ class multithread_adjacency_list< } const adjacency_list_type *m_ptr_adjacency_list; - std::size_t m_current_bank_index; - local_iterator_type m_local_iterator; + std::size_t m_current_bank_index; + local_iterator_type m_local_iterator; }; } // namespace clams From b97dde947252d6c7f6c3c03bb165c42e45b93489 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Tue, 15 Sep 2026 18:48:05 -0700 Subject: [PATCH 26/37] Change to use parallel PM datastore copy program when back up KNNG --- script/benchmark/run_clams_bench.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index 4aec0eb..b572393 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -148,6 +148,9 @@ def parse_options(): parser.add_argument('--noise_point_assigner_exe', default=f'{cwd}/src/clustering/cluster_noise_points', help='Path to the noise point assigner executable.') + parser.add_argument('--pm_datastore_copy_exe', + default=f'{cwd}/src/pm_datastore/copy_pm_datastore', + help='Path to the PM datastore copy executable.') # Etc parser.add_argument('-v', '--verbose', action='store_true', @@ -268,6 +271,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, dnnd_exe, nng_k, distance_func, points_file_format, point_path, nng_r, nng_delta, + pm_datastore_copy_exe, backup_knng, mfc_exe, amst_exe, amst_approx_bound_list, @@ -306,8 +310,8 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, if backup_knng: dnnd_ds_path_backup = f"{dnnd_ds_path}_backup" job_script.write(f"echo \"Backing up KNNG datastore\"\n") - backup_knng_command = f"cp -r {dnnd_ds_path} {dnnd_ds_path_backup}" - add_cmd(backup_knng_command, job_script) + backup_knng_command = f"{pm_datastore_copy_exe} -s {dnnd_ds_path} -t {dnnd_ds_path_backup}" + add_srun_cmd(num_tasks_per_node, backup_knng_command, job_script) else: job_script.write( f"Using existing DNND datastore at {input_dnnd_ds_path}\n") @@ -420,6 +424,7 @@ def main(): opts.point_path, opts.nng_r, opts.nng_delta, + opts.pm_datastore_copy_exe, opts.backup_knng, opts.mfc_exe, opts.amst_exe, From f489b5524d7b48926ff81e61a8c32f5ef2dd2527 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Tue, 15 Sep 2026 18:48:39 -0700 Subject: [PATCH 27/37] Bugfix: Boost component linking in CMakeLists.txt --- CMakeLists.txt | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1945028..2a353aa 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -81,11 +81,6 @@ else () COMPONENTS ${BOOST_COMPS}) endif () -# Boost Components to link -foreach (lib ${BOOST_INCLUDE_LIBRARIES}) - list(APPEND BOOST_COMPONENTS "Boost::${lib}") -endforeach () - # spdlog # FetchContent_Declare( @@ -200,7 +195,7 @@ function(add_basic_executable target_name src_name) # target_include_directories(${target_name} # PRIVATE ${PROJECT_SOURCE_DIR}/include) add_common_compile_options(${target_name}) - target_link_libraries(${target_name} PRIVATE Threads::Threads ${BOOST_COMPONENTS}) + target_link_libraries(${target_name} PRIVATE Threads::Threads ${BOOST_LIBS}) endfunction() # From 8c4c6439ea57d461cdf4c40fa80b07eeeb1efad8 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Tue, 15 Sep 2026 21:59:22 -0700 Subject: [PATCH 28/37] Use the distributed HDBSCAN by default --- script/benchmark/run_clams_bench.py | 82 ++++++++++--------------- src/clustering/CMakeLists.txt | 5 -- src/clustering/cluster_noise_points.cpp | 3 +- src/common.hpp | 45 ++++++++------ 4 files changed, 59 insertions(+), 76 deletions(-) diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index b572393..edf2228 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -86,9 +86,6 @@ def parse_options(): parser.add_argument('-s', '--min_samples', type=int, default=-1, dest='min_samples', help='min_sample value for calculating core distance. If -1 is given, core distance is not calculated.') - # Use distributed HDBSCAN - parser.add_argument('--distributed_hdbscan', action='store_true', - help='Use distributed HDBSCAN (run_distributed_hdbscan_clustering) instead of serial version (run_hdbscan_clustering).') # For evaluation parser.add_argument('-g', '--ground_truth_path', @@ -137,7 +134,7 @@ def parse_options(): default=f'{cwd}/src/mst/build_amst', help='Path to the AMST executable.') parser.add_argument('-C', '--clustering_exe', - default=f'{cwd}/src/clustering/run_hdbscan_clustering', + default=f'{cwd}/src/clustering/run_distributed_hdbscan_clustering', help='Path to the HPC clustering executable.') parser.add_argument('-E', '--python_evaluator', default=f'{cwd}/script/benchmark/evaluate_clustering_quality.py', @@ -165,20 +162,23 @@ def generate_job_name(): return f"job_{datetime.now().strftime('%Y%m%d_%H%M%S')}" stage_name_stack = [] # Stack to store stage names for logging -def echo_stage_name(job_script, stage_name): +# Show stage start message and stack the name in stage_name_stack +# This function must be paired with a call to finish_stage() to properly log the end of the stage. +def start_stage(job_script, stage_name): job_script.write("echo\n") job_script.write(f"echo ================================\n") job_script.write(f"echo \"{stage_name}\": $(date \"+%Y/%m/%d %H:%M:%S\")\n") job_script.write(f"echo ================================\n") stage_name_stack.append(stage_name) # Add the stage name to the queue +# Show stage finish message and pop the name from stage_name_stack def finish_stage(job_script): if stage_name_stack: stage_name = stage_name_stack.pop() # Get the last stage name from the stack job_script.write(f"echo \"Finished {stage_name}: $(date \"+%Y/%m/%d %H:%M:%S\")\"\n") job_script.write("echo\n") -def add_clustering_evaluation(job_script, cluster_label_file, amst_ds_path, +def add_clustering_evaluation(job_script, cluster_label_dir, amst_ds_path, ground_truth_path, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, @@ -188,20 +188,20 @@ def add_clustering_evaluation(job_script, cluster_label_file, amst_ds_path, verbose_flag = " -v" if verbose else "" evaluation_command = (f"{evaluator} {verbose_flag} " f"-g {ground_truth_path} " - f"{cluster_label_file}") + f"{cluster_label_dir}") add_srun_cmd(num_tasks_per_node, evaluation_command, job_script) else: job_script.write( "echo \"Evaluating Clustering using python script\"\n") evaluation_command = (f"python3 {evaluator} " - f"-c {cluster_label_file} " + f"-c {cluster_label_dir} " f"-g {ground_truth_path}") if singleton_cluster_to_noise_points: evaluation_command += " -s" add_cmd(evaluation_command, job_script) def run_clustering(job_script, set_cmd, work_dir, amst_approx_bound, - amst_ds_path, clustering_exe, distributed_hdbscan, + amst_ds_path, clustering_exe, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, ground_truth_path, singleton_cluster_to_noise_points, @@ -209,47 +209,37 @@ def run_clustering(job_script, set_cmd, work_dir, amst_approx_bound, # Set the min cluster size environment variable for this run add_cmd(set_cmd, job_script, False, False) - echo_stage_name(job_script, "Running CLAMS-HDBSCAN") + start_stage(job_script, "Running CLAMS-HDBSCAN") job_script.write(f"echo \"Min cluster size ${{MIN_CLUSTER_SIZE}}\"\n") - if distributed_hdbscan: - cluster_label_file = f"{work_dir}/cluster_labels_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}/" - cluster_tree_file = f"{work_dir}/cluster_tree_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}/" - verbose_flag = '-v' if verbose else '' - hpc_clustering_command = (f"{clustering_exe} {verbose_flag} -i {amst_ds_path} -M " - f" -m ${{MIN_CLUSTER_SIZE}} " - f" -o {cluster_label_file} " - f" -c {cluster_tree_file}") - add_srun_cmd(num_tasks_per_node, hpc_clustering_command, job_script) - else: - cluster_label_file = f"{work_dir}/cluster_labels_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}.txt" - cluster_tree_file = f"{work_dir}/cluster_tree_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}.txt" - hpc_clustering_command = (f"{clustering_exe} -i {amst_ds_path} -M " - f" -m ${{MIN_CLUSTER_SIZE}} " - f" -o {cluster_label_file} " - f" -c {cluster_tree_file} " - f" -P ") - add_cmd(hpc_clustering_command, job_script) + cluster_label_dir = f"{work_dir}/cluster_labels_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}/" + cluster_tree_dir = f"{work_dir}/cluster_tree_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}/" + verbose_flag = '-v' if verbose else '' + hpc_clustering_command = (f"{clustering_exe} {verbose_flag} -i {amst_ds_path} -M " + f" -m ${{MIN_CLUSTER_SIZE}} " + f" -o {cluster_label_dir} " + f" -c {cluster_tree_dir}") + add_srun_cmd(num_tasks_per_node, hpc_clustering_command, job_script) finish_stage(job_script) if ground_truth_path: - echo_stage_name(job_script, "Evaluating Clustering Results") + start_stage(job_script, "Evaluating Clustering Results") job_script.write(f"echo \"MST-based cluster guess: False\"\n") add_clustering_evaluation( - job_script, cluster_label_file, amst_ds_path, ground_truth_path, + job_script, cluster_label_dir, amst_ds_path, ground_truth_path, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, singleton_cluster_to_noise_points) finish_stage(job_script) - echo_stage_name(job_script, "Assign clusters to noise points") - cluster_label_file_no_noise = f"{cluster_label_file[:-4]}.noise_assigned.txt" + start_stage(job_script, "Assign clusters to noise points") + cluster_label_file_no_noise = f"{cluster_label_dir[:-4]}-noise_assigned.txt" cluster_assign_command = (f"{noise_point_assigner_exe} -M " f"-m {amst_ds_path} " - f"-c {cluster_label_file} " + f"-c {cluster_label_dir} " f"-o {cluster_label_file_no_noise}") add_cmd(cluster_assign_command, job_script) finish_stage(job_script) - echo_stage_name(job_script, "Evaluate clustering results") + start_stage(job_script, "Evaluate clustering results") job_script.write(f"echo \"MST-based cluster guess: True\"\n") add_clustering_evaluation( job_script, cluster_label_file_no_noise, amst_ds_path, @@ -275,7 +265,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, backup_knng, mfc_exe, amst_exe, amst_approx_bound_list, - clustering_exe, distributed_hdbscan, + clustering_exe, evaluator, ygm_cluster_eval, verbose, ground_truth_path, @@ -300,7 +290,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, # Run the DNND step job_script.write(f"echo \"kNNG k: {nng_k}\"\n") if len(input_dnnd_ds_path) == 0: - echo_stage_name(job_script, "Building KNNG") + start_stage(job_script, "Building KNNG") dnnd_ds_path = f"{work_dir}/dnnd_pm_datastore" dnnd_batch_size = 2 ** 25 verbose_flag = '-v' if verbose else '' @@ -318,7 +308,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, dnnd_ds_path = input_dnnd_ds_path # Connect the CCs - echo_stage_name(job_script, "Connecting the CCs using MFC") + start_stage(job_script, "Connecting the CCs using MFC") mfc_command = f"{mfc_exe} -d {dnnd_ds_path} -f {distance_func}" add_srun_cmd(num_tasks_per_node, mfc_command, job_script) finish_stage(job_script) @@ -326,7 +316,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, # Convert to core distance # TODO: Implement if False and min_samples > 0: - echo_stage_name(job_script, "Convert to core distance kNNG") + start_stage(job_script, "Convert to core distance kNNG") knng_coredist_dir = f"{work_dir}/knng_coredist/" add_cmd(f'mkdir -p {knng_coredist_dir}', job_script) conv2coredist_cmd = f"./src/conv_knng_to_core_dist -i {dnnd_ds_path} -o {knng_coredist_dir}/knng.txt -m {min_samples}" @@ -336,7 +326,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, try_no = 0 for amst_approx_bound in amst_approx_bound_list: # Run the AMST step - echo_stage_name(job_script, f"Running AMST, approx bound = {amst_approx_bound}") + start_stage(job_script, f"Running AMST, approx bound = {amst_approx_bound}") amst_ds_path = f"{work_dir}/amst_pm_datastore_a{amst_approx_bound}" amst_command = f"{amst_exe} -d {dnnd_ds_path} -p {amst_ds_path} -e {amst_approx_bound}" add_srun_cmd(num_tasks_per_node, amst_command, job_script) @@ -346,7 +336,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, for set_cmd in min_cluster_size_set_cmnds: run_clustering( job_script, set_cmd, work_dir, amst_approx_bound, - amst_ds_path, clustering_exe, distributed_hdbscan, + amst_ds_path, clustering_exe, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, ground_truth_path, singleton_cluster_to_noise_points, noise_point_assigner_exe) @@ -394,15 +384,6 @@ def main(): # This is not the best way to set the number of threads for NEO-DNND, but it is a simple way to do it for now. dnnd_exe = f'{dnnd_exe} -T {opts.neodnnd_threads} -R {opts.neodnnd_replicate_rate}' - # Select HDBSCAN executable - # If --distributed_hdbscan is sepecified and the user did not override --clustering_exe, - # switch the default to run_distributed_hdbscan_clustering - clustering_exe = opts.clustering_exe - if opts.distributed_hdbscan: - default_clustering_exe = f'{os.getcwd()}/src/clustering/run_hdbscan_clustering' - if clustering_exe == default_clustering_exe: - clustering_exe = f'{os.getcwd()}/src/clustering/run_distributed_hdbscan_clustering' - if opts.ygm_cluster_eval: evaluator = opts.ygm_evaluator_exe else: @@ -429,8 +410,7 @@ def main(): opts.mfc_exe, opts.amst_exe, amst_approx_bound_list, - clustering_exe, - opts.distributed_hdbscan, + opts.clustering_exe, evaluator, opts.ygm_cluster_eval, opts.verbose, diff --git a/src/clustering/CMakeLists.txt b/src/clustering/CMakeLists.txt index de738e1..59a8ed2 100644 --- a/src/clustering/CMakeLists.txt +++ b/src/clustering/CMakeLists.txt @@ -12,11 +12,6 @@ setup_omp_target(cluster_noise_points) setup_metall_target(cluster_noise_points) setup_spdlog_target(cluster_noise_points) -add_basic_executable(cluster_noise_points_v2 cluster_noise_points_v2.cpp) -setup_omp_target(cluster_noise_points_v2) -setup_metall_target(cluster_noise_points_v2) -setup_spdlog_target(cluster_noise_points_v2) - add_basic_executable(run_distributed_hdbscan_clustering run_distributed_hdbscan_clustering.cpp) setup_ygm_target(run_distributed_hdbscan_clustering) setup_omp_target(run_distributed_hdbscan_clustering) diff --git a/src/clustering/cluster_noise_points.cpp b/src/clustering/cluster_noise_points.cpp index 316268a..43ab552 100644 --- a/src/clustering/cluster_noise_points.cpp +++ b/src/clustering/cluster_noise_points.cpp @@ -116,7 +116,8 @@ int main(int argc, char* argv[]) { } map_t point_cluster_map; - read_cluster_ids(opt.cluster_ids_input_path, point_cluster_map); + read_cluster_ids(clams::find_files(opt.cluster_ids_input_path), + point_cluster_map); spdlog::info("Read {} points' cluster IDs from {}", point_cluster_map.size(), opt.cluster_ids_input_path.string()); diff --git a/src/common.hpp b/src/common.hpp index 705e429..7a7b257 100644 --- a/src/common.hpp +++ b/src/common.hpp @@ -174,30 +174,37 @@ inline void read_knng_edges( /// ignored. The output maps each point ID to its cluster ID; if an ID occurs /// more than once, the last assignment wins. template -void read_cluster_ids(const std::filesystem::path &input_path, - cluster_id_table_t &point_cluster_map) { +void read_cluster_ids(const std::vector &input_paths, + cluster_id_table_t &point_cluster_map) { using id_t = typename cluster_id_table_t::key_type; using cluster_id_t = typename cluster_id_table_t::mapped_type; - std::ifstream ifs(input_path); - if (!ifs) { - std::cerr << "Failed to open " << input_path << std::endl; - std::abort(); - } - - std::string line; - while (std::getline(ifs, line)) { - if (line.empty() || line[0] == '#') { - continue; // Skip empty lines and comments - } - std::istringstream iss(line); - id_t point_id; - cluster_id_t cluster_id; - if (!(iss >> point_id >> cluster_id)) { - std::cerr << "Error parsing line: " << line << std::endl; + for (const auto &input_path : input_paths) { + std::ifstream ifs(input_path); + if (!ifs) { + std::cerr << "Failed to open " << input_path << std::endl; std::abort(); } - point_cluster_map[point_id] = cluster_id; + + std::string line; + while (std::getline(ifs, line)) { + if (line.empty() || line[0] == '#') { + continue; // Skip empty lines and comments + } + std::istringstream iss(line); + id_t point_id; + cluster_id_t cluster_id; + if (!(iss >> point_id >> cluster_id)) { + std::cerr << "Error parsing line: " << line << std::endl; + std::abort(); + } + if (point_cluster_map.find(point_id) != point_cluster_map.end()) { + std::cerr << "Error: point ID " << point_id + << " is assigned to multiple clusters." << std::endl; + std::abort(); + } + point_cluster_map[point_id] = cluster_id; + } } } From 6e6f2bdcffa0cca4883e182922e839c42abe7259 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Tue, 15 Sep 2026 22:31:22 -0700 Subject: [PATCH 29/37] Enhance KNNG evaluation scripts with no-distance option and improve output reporting --- .../evaluate_correlation_knng_and_clusters.py | 6 +++++ src/knng/serialize_pm_knng.cpp | 26 ++++++++++++------- ...evaluate_correlation_knng_and_clusters.cpp | 20 +++++++++----- 3 files changed, 36 insertions(+), 16 deletions(-) diff --git a/script/benchmark/evaluate_correlation_knng_and_clusters.py b/script/benchmark/evaluate_correlation_knng_and_clusters.py index b41784b..7ddc7e4 100644 --- a/script/benchmark/evaluate_correlation_knng_and_clusters.py +++ b/script/benchmark/evaluate_correlation_knng_and_clusters.py @@ -10,14 +10,17 @@ Arguments: -g/--graph: KNNG file or directory containing KNNG dump files. + -N/--no-distance: indicate that the KNNG files do not contain distance lines. -c/--cluster: cluster file or directory containing cluster-label files. -k/--neighbors: number of neighbors to consider for each point. + -h/--help: show this help message and exit. KNNG file format: Each record is encoded as two lines: ... ... The leading distance is a dummy value corresponding to the source ID and will be ignored. + If the KNNG files do not contain distance lines, use the -N/--no-distance flag. Cluster file format: Each line contains a point and its cluster ID separated by whitespace. @@ -263,6 +266,9 @@ def main() -> None: for k in args.neighbors: try: per_point, summary = evaluate_knng_clusters(args.graph, args.cluster, k, no_distance=args.no_distance) + print(f"\nResults for k={k}:") + print(summary["num_points"], "points evaluated,", summary["num_noise_points"], "noise points skipped.") + print(summary["num_same_cluster_neighbors"], "same-cluster neighbors out of", summary["num_considered_neighbors"], "considered neighbors,", "fraction =", summary["avg_same_cluster_fraction"]) except Exception as exc: # pragma: no cover - CLI error path raise SystemExit(f"Error: {exc}") from exc diff --git a/src/knng/serialize_pm_knng.cpp b/src/knng/serialize_pm_knng.cpp index 93e05b8..48facc9 100644 --- a/src/knng/serialize_pm_knng.cpp +++ b/src/knng/serialize_pm_knng.cpp @@ -1,7 +1,6 @@ // Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS // Project Developers. See the top-level COPYRIGHT file for details. - #define CLAMS_USE_SALTATLAS #define METALL_DISABLE_CONCURRENCY @@ -9,11 +8,11 @@ #include #include #include -#include -#include #include +#include #include #include +#include #include @@ -22,9 +21,9 @@ using namespace clams; void parse_option(int argc, char *argv[], std::filesystem::path &pm_knng_path, - std::filesystem::path &output_path) { + std::filesystem::path &output_path, bool &no_distance_dump) { int opt_char; - while ((opt_char = getopt(argc, argv, "i:o:h")) != -1) { + while ((opt_char = getopt(argc, argv, "i:o:Nh")) != -1) { switch (opt_char) { case 'i': pm_knng_path = std::filesystem::path(optarg); @@ -32,6 +31,10 @@ void parse_option(int argc, char *argv[], std::filesystem::path &pm_knng_path, case 'o': output_path = std::filesystem::path(optarg); break; + case 'N': + no_distance_dump = true; + break; + case 'h': std::cout << "Usage: " << argv[0] << " -i -o " << std::endl; @@ -59,7 +62,8 @@ int main(int argc, char *argv[]) { std::filesystem::path pm_knng_path; std::filesystem::path output_path; - parse_option(argc, argv, pm_knng_path, output_path); + bool no_distance_dump = false; + parse_option(argc, argv, pm_knng_path, output_path, no_distance_dump); dist_pm_knng_t pm_knng(comm.get_mpi_comm()); pm_knng.open_read_only(pm_knng_path); @@ -83,11 +87,13 @@ int main(int argc, char *argv[]) { } ofs << "\n"; - ofs << "0.0"; - for (const auto &neighbor : neighbors) { - ofs << " " << neighbor.distance; + if (!no_distance_dump) { + ofs << "0.0"; + for (const auto &neighbor : neighbors) { + ofs << " " << neighbor.distance; + } + ofs << "\n"; } - ofs << "\n"; } comm.cf_barrier(); diff --git a/src/tools/evaluate_correlation_knng_and_clusters.cpp b/src/tools/evaluate_correlation_knng_and_clusters.cpp index 35178a4..903a22e 100644 --- a/src/tools/evaluate_correlation_knng_and_clusters.cpp +++ b/src/tools/evaluate_correlation_knng_and_clusters.cpp @@ -33,11 +33,12 @@ using cluster_id_table_t = boost::unordered::unordered_flat_map; // correlation analysis. bool parse_option(int argc, char *argv[], std::filesystem::path &input_knng_path, - std::filesystem::path &input_clusters_path, size_t &max_k) { + std::filesystem::path &input_clusters_path, size_t &max_k, + bool &no_distance) { input_knng_path.clear(); int opt; - while ((opt = ::getopt(argc, argv, "g:c:k:")) != -1) { + while ((opt = ::getopt(argc, argv, "g:c:k:N")) != -1) { switch (opt) { case 'g': { input_knng_path = std::filesystem::path(optarg); @@ -51,6 +52,10 @@ bool parse_option(int argc, char *argv[], max_k = std::stoul(optarg); break; } + case 'N': { + no_distance = true; + break; + } default: { std::cerr << "Unknown option: " << opt << std::endl; return false; @@ -74,9 +79,11 @@ bool parse_option(int argc, char *argv[], int main(int argc, char *argv[]) { std::filesystem::path input_knng_path; std::filesystem::path input_clusters_path; - size_t max_k = 100; // default value + size_t max_k = 100; // default value + bool no_distance = false; - if (!parse_option(argc, argv, input_knng_path, input_clusters_path, max_k)) { + if (!parse_option(argc, argv, input_knng_path, input_clusters_path, max_k, + no_distance)) { return EXIT_FAILURE; } @@ -85,12 +92,13 @@ int main(int argc, char *argv[]) { clams::shm_graph_t graph; std::cout << "Read knng" << std::endl; - clams::read_knng(knng_files, graph); + clams::read_knng(knng_files, graph, !no_distance); std::cout << "#of points: " << graph.num_keys() << std::endl; std::cout << "#of neighbors: " << graph.num_values() << std::endl; cluster_id_table_t cluster_id_table; - clams::read_cluster_ids(input_clusters_path, cluster_id_table); + clams::read_cluster_ids(clams::find_files(input_clusters_path), + cluster_id_table); std::cout << "#of points with clusters: " << cluster_id_table.size() << std::endl; const auto &cluster_ids = cluster_id_table; From 918f0f8f223f32a05587b91d3178c7a23f43c993 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Tue, 15 Sep 2026 22:38:46 -0700 Subject: [PATCH 30/37] Add support for optional distance reading in knng file processing --- src/details/shm_graph.hpp | 79 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 76 insertions(+), 3 deletions(-) diff --git a/src/details/shm_graph.hpp b/src/details/shm_graph.hpp index cc682ac..280c7a3 100644 --- a/src/details/shm_graph.hpp +++ b/src/details/shm_graph.hpp @@ -11,6 +11,8 @@ #include #include +#include + #include "data_types.hpp" #include "multithread_adjacency_list.hpp" @@ -21,6 +23,8 @@ using shm_graph_t = /// \brief Read knng files. /// \param knng_files A list of knng files. /// \param graph A graph to store the knng. +/// \param has_distance If true, read the distances in the knng files; if false, +/// ignore them. /// \details Each knng file is expected to have the following format: /// - Each even numbered line (0-based) contains a source vertex ID followed /// by its neighbor vertex IDs, separated by whitespace. @@ -30,9 +34,9 @@ using shm_graph_t = /// itself, and the corresponding distance in the odd numbered line is 0.0 /// (dummy distance). /// - The number of neighbor IDs in each even numbered line must match the -/// number of distances in the corresponding odd numbered +/// number of distances in the corresponding odd numbered line. void read_knng(const std::vector &knng_files, - shm_graph_t &graph) { + shm_graph_t &graph, const bool has_distance = true) { OMP_DIRECTIVE(parallel for) for (std::size_t fno = 0; fno < knng_files.size(); ++fno) { const auto &file = knng_files[fno]; @@ -55,16 +59,82 @@ void read_knng(const std::vector &knng_files, } std::vector dists; + if (has_distance) { + std::getline(ifs, line); + std::istringstream iss(line); + distance_t buf; + while (iss >> buf) { + dists.push_back(buf); + } + } else { + dists.resize(ids.size()); + } + + if (has_distance && ids.size() != dists.size()) { + std::cerr << "Invalid file: " << file << std::endl; + std::cerr << "#of IDs and distances do not match" << std::endl; + std::exit(1); + } + + if (ids.empty()) { + break; // End of file + } + + const id_t src = ids[0]; + for (std::size_t i = 1; i < ids.size(); ++i) { + graph.add(src, std::make_pair(ids[i], dists[i])); + } + } + } +} + +/// \brief Read knng files, keeping only the rows whose source vertex ID is +/// in `keep_ids`. +/// \param knng_files A list of knng files. +/// \param keep_ids The set of source vertex IDs to keep; rows for other +/// source IDs are skipped. +/// \param graph A graph to store the knng. +/// \param has_distance If true, read the distances in the knng files; if +/// false, ignore them. +/// \details Uses the same file format as \ref read_knng. +void read_knng_filtered( + const std::vector &knng_files, + const boost::unordered::unordered_flat_set &keep_ids, + shm_graph_t &graph, const bool has_distance = true) { + OMP_DIRECTIVE(parallel for) + for (std::size_t fno = 0; fno < knng_files.size(); ++fno) { + const auto &file = knng_files[fno]; + std::ifstream ifs(file); + if (!ifs.is_open()) { + std::cerr << "Cannot open file: " << file << std::endl; + std::exit(1); + } + + std::string line; + while (true) { + std::vector ids; { + std::getline(ifs, line); + std::istringstream iss(line); + id_t buf; + while (iss >> buf) { + ids.push_back(buf); + } + } + + std::vector dists; + if (has_distance) { std::getline(ifs, line); std::istringstream iss(line); distance_t buf; while (iss >> buf) { dists.push_back(buf); } + } else { + dists.resize(ids.size()); } - if (ids.size() != dists.size()) { + if (has_distance && ids.size() != dists.size()) { std::cerr << "Invalid file: " << file << std::endl; std::cerr << "#of IDs and distances do not match" << std::endl; std::exit(1); @@ -75,6 +145,9 @@ void read_knng(const std::vector &knng_files, } const id_t src = ids[0]; + if (!keep_ids.contains(src)) { + continue; // Source vertex not of interest; skip this row + } for (std::size_t i = 1; i < ids.size(); ++i) { graph.add(src, std::make_pair(ids[i], dists[i])); } From b612336c787f3d41f495c16e26f708b729cb43cf Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Tue, 15 Sep 2026 23:22:56 -0700 Subject: [PATCH 31/37] Change the header format of clustering output file --- src/clustering/run_distributed_hdbscan_clustering.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/clustering/run_distributed_hdbscan_clustering.cpp b/src/clustering/run_distributed_hdbscan_clustering.cpp index 78943a9..42935a8 100644 --- a/src/clustering/run_distributed_hdbscan_clustering.cpp +++ b/src/clustering/run_distributed_hdbscan_clustering.cpp @@ -1180,7 +1180,7 @@ int main(int argc, char *argv[]) { std::filesystem::path label_file_path = opt.cluster_ids_out_path / label_file_name; std::ofstream labels_ofs(label_file_path); - labels_ofs << "point_id\tcluster_id\n"; + labels_ofs << "# point_id\tcluster_id\n"; for (auto &[point_id, cluster_id] : point_to_cluster_id_map) { labels_ofs << point_id << "\t" << cluster_id << "\n"; From 9d88646248d770e03270877c932914198d68f593 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Wed, 16 Sep 2026 01:20:37 -0700 Subject: [PATCH 32/37] Change noise clustering algorithm to consider distance --- script/benchmark/run_clams_bench.py | 60 +++++++++++--------- src/clustering/cluster_noise_points.cpp | 74 +++++++++++++++---------- 2 files changed, 79 insertions(+), 55 deletions(-) diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py index edf2228..66b33a3 100644 --- a/script/benchmark/run_clams_bench.py +++ b/script/benchmark/run_clams_bench.py @@ -96,6 +96,8 @@ def parse_options(): help='Use YGM cluster evaluation calculation.') parser.add_argument('-z', '--dummy_cluster_id', action='store_true', help='Assign a singleton cluster to each noise point in the clustering result.') + parser.add_argument('--cluster_noise_points', action='store_true', + help='Enable clustering of noise points after the main clustering process.') # For output parser.add_argument('-o', '--output_root_dir', @@ -142,11 +144,11 @@ def parse_options(): parser.add_argument('-Y', '--ygm_evaluator_exe', default=f'{cwd}/tpls/clams-cc/build/src/clustering_metrics', help='Path to the YGM clustering evaluation executable.') - parser.add_argument('--noise_point_assigner_exe', + parser.add_argument('--noise_point_clustering_exe', default=f'{cwd}/src/clustering/cluster_noise_points', help='Path to the noise point assigner executable.') parser.add_argument('--pm_datastore_copy_exe', - default=f'{cwd}/src/pm_datastore/copy_pm_datastore', + default=f'{cwd}/src/tools/copy_pm_datastore', help='Path to the PM datastore copy executable.') # Etc @@ -205,7 +207,8 @@ def run_clustering(job_script, set_cmd, work_dir, amst_approx_bound, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, ground_truth_path, singleton_cluster_to_noise_points, - noise_point_assigner_exe): + noise_point_clustering_exe, + cluster_noise_points): # Set the min cluster size environment variable for this run add_cmd(set_cmd, job_script, False, False) @@ -221,6 +224,16 @@ def run_clustering(job_script, set_cmd, work_dir, amst_approx_bound, add_srun_cmd(num_tasks_per_node, hpc_clustering_command, job_script) finish_stage(job_script) + if cluster_noise_points: + start_stage(job_script, "Assign clusters to noise points") + no_noise_points_cluster_label_file = f"{work_dir}/cluster_labels_a{amst_approx_bound}_m${{MIN_CLUSTER_SIZE}}_no_noise.txt" + cluster_assign_command = (f"{noise_point_clustering_exe} -M " + f"-m {amst_ds_path} " + f"-c {cluster_label_dir} " + f"-o {no_noise_points_cluster_label_file}") + add_cmd(cluster_assign_command, job_script) + finish_stage(job_script) + if ground_truth_path: start_stage(job_script, "Evaluating Clustering Results") job_script.write(f"echo \"MST-based cluster guess: False\"\n") @@ -230,22 +243,14 @@ def run_clustering(job_script, set_cmd, work_dir, amst_approx_bound, singleton_cluster_to_noise_points) finish_stage(job_script) - start_stage(job_script, "Assign clusters to noise points") - cluster_label_file_no_noise = f"{cluster_label_dir[:-4]}-noise_assigned.txt" - cluster_assign_command = (f"{noise_point_assigner_exe} -M " - f"-m {amst_ds_path} " - f"-c {cluster_label_dir} " - f"-o {cluster_label_file_no_noise}") - add_cmd(cluster_assign_command, job_script) - finish_stage(job_script) - - start_stage(job_script, "Evaluate clustering results") - job_script.write(f"echo \"MST-based cluster guess: True\"\n") - add_clustering_evaluation( - job_script, cluster_label_file_no_noise, amst_ds_path, - ground_truth_path, evaluator, ygm_cluster_eval, - num_tasks_per_node, verbose, singleton_cluster_to_noise_points) - finish_stage(job_script) + if cluster_noise_points: + start_stage(job_script, "Evaluate clustering results") + job_script.write(f"echo \"MST-based cluster guess: True\"\n") + add_clustering_evaluation( + job_script, no_noise_points_cluster_label_file, amst_ds_path, + ground_truth_path, evaluator, ygm_cluster_eval, + num_tasks_per_node, verbose, singleton_cluster_to_noise_points) + finish_stage(job_script) job_script.write("echo \"\" \n") @@ -267,13 +272,15 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, amst_exe, amst_approx_bound_list, clustering_exe, evaluator, - ygm_cluster_eval, verbose, + ygm_cluster_eval, ground_truth_path, singleton_cluster_to_noise_points, min_cluster_size_set_cmnds, min_samples, - noise_point_assigner_exe, - input_dnnd_ds_path='',): + noise_point_clustering_exe, + cluster_noise_points, + input_dnnd_ds_path, + verbose): create_dir(job_dir) job_script_path = f'{job_dir}/job.sh' @@ -339,7 +346,7 @@ def gen_clams_bench_script(job_name, job_dir, work_dir, amst_ds_path, clustering_exe, evaluator, ygm_cluster_eval, num_tasks_per_node, verbose, ground_truth_path, singleton_cluster_to_noise_points, - noise_point_assigner_exe) + noise_point_clustering_exe, cluster_noise_points) try_no += 1 finish_stage(job_script) @@ -413,13 +420,14 @@ def main(): opts.clustering_exe, evaluator, opts.ygm_cluster_eval, - opts.verbose, opts.ground_truth_path, opts.dummy_cluster_id, min_cluster_size_set_cmnds, opts.min_samples, - opts.noise_point_assigner_exe, - opts.input_dnnd_ds_path) + opts.noise_point_clustering_exe, + opts.cluster_noise_points, + opts.input_dnnd_ds_path, + opts.verbose) print(f"Generated batch script: {job_script}") job_submission_cmd = f"sbatch {opts.sbatch_opts} {job_script}" diff --git a/src/clustering/cluster_noise_points.cpp b/src/clustering/cluster_noise_points.cpp index 43ab552..e2c37cd 100644 --- a/src/clustering/cluster_noise_points.cpp +++ b/src/clustering/cluster_noise_points.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include "../common.hpp" @@ -30,6 +31,9 @@ using namespace clams; template using map_t = boost::unordered::unordered_flat_map; +template +using set_t = boost::unordered::unordered_flat_set; + struct option { std::filesystem::path mst_edges_path; bool metall_mst{false}; @@ -80,7 +84,7 @@ int main(int argc, char* argv[]) { option opt; parse_option(argc, argv, opt); - map_t> mst_graph; + map_t>> mst_graph; if (opt.metall_mst) { spdlog::info("Attaching MST in Metall datastore"); metall::manager metall_manager(metall::open_read_only, opt.mst_edges_path); @@ -95,8 +99,8 @@ int main(int argc, char* argv[]) { spdlog::info("#of MST edges: {}", input_mst_edges->size()); spdlog::info("Copying MST edges from Metall datastore"); for (const auto& edge : *input_mst_edges) { - mst_graph[edge.ids[0]].push_back(edge.ids[1]); - mst_graph[edge.ids[1]].push_back(edge.ids[0]); + mst_graph[edge.ids[0]].push_back({edge.ids[1], edge.distance}); + mst_graph[edge.ids[1]].push_back({edge.ids[0], edge.distance}); } } else { spdlog::info("Reading MST edges"); @@ -104,8 +108,8 @@ int main(int argc, char* argv[]) { read_edges(opt.mst_edges_path, input_mst_edges); spdlog::info("#of MST edges: {}", input_mst_edges.size()); for (const auto& edge : input_mst_edges) { - mst_graph[edge.ids[0]].push_back(edge.ids[1]); - mst_graph[edge.ids[1]].push_back(edge.ids[0]); + mst_graph[edge.ids[0]].push_back({edge.ids[1], edge.distance}); + mst_graph[edge.ids[1]].push_back({edge.ids[0], edge.distance}); } } @@ -135,9 +139,9 @@ int main(int argc, char* argv[]) { std::size_t n_noise_points = 0; std::size_t n_assigned_points = 0; - // NOTE: this algorithm is not determinstic because threads update the shared - // point_cluster_map concurrently. - // We employ this algorithm because it is simple and fast. + // NOTE: point_cluster_map may be updated concurrently by multiple threads, + // which causes data races, but it is acceptable for this use case. We employ + // this algorithm because it is simple and fast. OMP_DIRECTIVE(parallel for reduction(+ : n_noise_points, n_assigned_points)) for (size_t i = 0; i < point_ids.size(); ++i) { const auto point_id = point_ids.at(i); @@ -146,32 +150,44 @@ int main(int argc, char* argv[]) { } ++n_noise_points; - std::deque bfs_queue; - map_t visited; - bfs_queue.push_back(point_id); - visited[point_id] = true; - bool found_cluster = false; + std::vector> bfs_front; - while (!bfs_queue.empty() && !found_cluster) { - const auto current_point_id = bfs_queue.front(); - bfs_queue.pop_front(); + set_t visited; + bfs_front.push_back({point_id, 0}); + visited.insert(point_id); + bool found_cluster = false; - for (const auto neighbor_id : mst_graph.at(current_point_id)) { - if (visited.find(neighbor_id) != visited.end()) { - continue; // Already visited, e.g., the node we came from - } - if (point_cluster_map.at(neighbor_id) != k_noise_cluster_id) { - // Found a neighbor that belongs to a cluster - point_cluster_map[point_id] = point_cluster_map.at(neighbor_id); - found_cluster = true; - ++n_assigned_points; - break; - } else { - visited[neighbor_id] = true; - bfs_queue.push_back(neighbor_id); + // Level-order traversal (BFS) to find the nearest cluster for the noise + // point. Within each level, visit closer neighbors first + while (!bfs_front.empty()) { + std::vector> next; + for (const auto& [pid, _] : bfs_front) { + // Traverse the neighbors of the current point in the MST + for (const auto& [nid, ndist] : mst_graph.at(pid)) { + if (visited.count(nid) > 0) { + continue; // Already visited, e.g., the node we came from + } + if (point_cluster_map.at(nid) != k_noise_cluster_id) { + // Found a neighbor that belongs to a cluster + // Update the table without locking, which may cause data races but + // is acceptable for this use case + point_cluster_map.at(point_id) = point_cluster_map.at(nid); + found_cluster = true; + ++n_assigned_points; + goto BFS_COMPLETE; + } else { + visited.insert(nid); + next.push_back({nid, ndist}); + } } } + bfs_front.swap(next); + // Sort the BFS front by distance to prioritize closer neighbors + std::sort( + bfs_front.begin(), bfs_front.end(), + [](const auto& a, const auto& b) { return a.second < b.second; }); } + BFS_COMPLETE: if (!found_cluster) { spdlog::warn("Point {} could not be assigned to any cluster.", point_id); From bb50b72c66962daa4554d22e30b13c25708955b9 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Wed, 16 Sep 2026 21:14:15 -0700 Subject: [PATCH 33/37] Optimize noise point clustering code --- src/clustering/cluster_noise_points.cpp | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/src/clustering/cluster_noise_points.cpp b/src/clustering/cluster_noise_points.cpp index e2c37cd..8279046 100644 --- a/src/clustering/cluster_noise_points.cpp +++ b/src/clustering/cluster_noise_points.cpp @@ -25,6 +25,7 @@ #include #include "../common.hpp" +#include "../details/multithread_adjacency_list.hpp" using namespace clams; @@ -84,7 +85,7 @@ int main(int argc, char* argv[]) { option opt; parse_option(argc, argv, opt); - map_t>> mst_graph; + multithread_adjacency_list> mst_graph; if (opt.metall_mst) { spdlog::info("Attaching MST in Metall datastore"); metall::manager metall_manager(metall::open_read_only, opt.mst_edges_path); @@ -98,18 +99,22 @@ int main(int argc, char* argv[]) { } spdlog::info("#of MST edges: {}", input_mst_edges->size()); spdlog::info("Copying MST edges from Metall datastore"); - for (const auto& edge : *input_mst_edges) { - mst_graph[edge.ids[0]].push_back({edge.ids[1], edge.distance}); - mst_graph[edge.ids[1]].push_back({edge.ids[0], edge.distance}); + OMP_DIRECTIVE(parallel for) + for (size_t i = 0; i < input_mst_edges->size(); ++i) { + const auto& edge = input_mst_edges->at(i); + mst_graph.add(edge.ids[0], {edge.ids[1], edge.distance}); + mst_graph.add(edge.ids[1], {edge.ids[0], edge.distance}); } } else { spdlog::info("Reading MST edges"); weighted_edge_list_t input_mst_edges; read_edges(opt.mst_edges_path, input_mst_edges); spdlog::info("#of MST edges: {}", input_mst_edges.size()); - for (const auto& edge : input_mst_edges) { - mst_graph[edge.ids[0]].push_back({edge.ids[1], edge.distance}); - mst_graph[edge.ids[1]].push_back({edge.ids[0], edge.distance}); + OMP_DIRECTIVE(parallel for) + for (size_t i = 0; i < input_mst_edges.size(); ++i) { + const auto& edge = input_mst_edges.at(i); + mst_graph.add(edge.ids[0], {edge.ids[1], edge.distance}); + mst_graph.add(edge.ids[1], {edge.ids[0], edge.distance}); } } @@ -163,7 +168,9 @@ int main(int argc, char* argv[]) { std::vector> next; for (const auto& [pid, _] : bfs_front) { // Traverse the neighbors of the current point in the MST - for (const auto& [nid, ndist] : mst_graph.at(pid)) { + for (auto nitr = mst_graph.values_begin(pid); + nitr != mst_graph.values_end(pid); ++nitr) { + const auto& [nid, ndist] = *nitr; if (visited.count(nid) > 0) { continue; // Already visited, e.g., the node we came from } From 6f0efb86aec68be853ae2e2ddaa77e1d2de1dc13 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Wed, 16 Sep 2026 21:14:58 -0700 Subject: [PATCH 34/37] Separate point reading time in build-knng-neo-dnnd --- src/knng/build_knng_neo.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/knng/build_knng_neo.cpp b/src/knng/build_knng_neo.cpp index 6ea98ee..0dbd9dd 100644 --- a/src/knng/build_knng_neo.cpp +++ b/src/knng/build_knng_neo.cpp @@ -69,16 +69,20 @@ int main(int argc, char **argv) { return 1; } - ygm::utility::timer neo_dnnd_const_timer; // Build kNNG using neo_dnnd neo_dnnd_t neo_dnnd( saltatlas::distance::distance_function(opt.distance_name), comm, opt.verbose); + + ygm::utility::timer read_points_timer; comm.cout0() << "\n<>" << std::endl; neo_dnnd.load_points(point_files.begin(), point_files.end(), opt.point_file_format); + comm.cout0() << "\nReading points took (s)\t" << read_points_timer.elapsed() + << std::endl; + ygm::utility::timer neo_dnnd_const_timer; comm.cout0() << "\n<>" << std::endl; auto knng = neo_dnnd.build(opt.index_k, opt.r, opt.delta, opt.replicate_rate, opt.batch_size); From 02237d8831a68c6c40a502c33b94df334d1678ff Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Fri, 2 Oct 2026 13:55:29 -0700 Subject: [PATCH 35/37] Minor bugfix in benchmark utilities --- script/benchmark/clustering_utilities.py | 2 +- script/benchmark/extract_clams_clustering_qualities.sh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/script/benchmark/clustering_utilities.py b/script/benchmark/clustering_utilities.py index 81646e9..93777b4 100644 --- a/script/benchmark/clustering_utilities.py +++ b/script/benchmark/clustering_utilities.py @@ -220,7 +220,7 @@ def read_label_data(data_path): label = int(line) if pid in labels_dict: - print(f"Duplicate ID in the ground truth file: {pid}") + print(f"Duplicate ID in the input files: {pid}") exit(1) labels_dict[pid] = label diff --git a/script/benchmark/extract_clams_clustering_qualities.sh b/script/benchmark/extract_clams_clustering_qualities.sh index 74ff9d9..661ddbb 100755 --- a/script/benchmark/extract_clams_clustering_qualities.sh +++ b/script/benchmark/extract_clams_clustering_qualities.sh @@ -34,7 +34,7 @@ for file in "$@"; do /kNNG k:/ || /^[[:space:]]*k:[[:space:]]*[0-9]+[[:space:]]*$/ { k = lastfield($0); next } /Min cluster size/ { mcs = lastfield($0); next } /MST-based cluster guess/ { mst_clustering = lastfield($0); next } - /#of final clusters/ { fc = lastfield($0); next } + /Number of flat clusters:/ { fc = lastfield($0); next } /Cluster coverage \(%\)/ { coverage = lastfield($0); next } /ARI/ { ari = lastfield($0); have_ari = 1; next } /AMI/ { From 905cd5237df7503058c3f7087833289d49562348 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Sun, 4 Oct 2026 12:57:03 -0700 Subject: [PATCH 36/37] Reorganizing python scripts --- script/benchmark/bench_utilities.py | 64 +---- script/benchmark/clustering_utilities.py | 162 +---------- .../benchmark/evaluate_clustering_quality.py | 6 +- .../evaluate_correlation_knng_and_clusters.py | 105 +------- script/benchmark/run_hdbscan_bench.py | 6 +- script/data_readers.py | 255 ++++++++++++++++++ script/etc/reassign_sequential_cluster_ids.py | 5 +- script/utilities.py | 178 ++++++++++++ 8 files changed, 450 insertions(+), 331 deletions(-) create mode 100644 script/data_readers.py create mode 100644 script/utilities.py diff --git a/script/benchmark/bench_utilities.py b/script/benchmark/bench_utilities.py index 4d16045..4435f02 100644 --- a/script/benchmark/bench_utilities.py +++ b/script/benchmark/bench_utilities.py @@ -1,44 +1,10 @@ # Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS # Project Developers. See the top-level COPYRIGHT file for details. - import subprocess import os - -# Execute a shell command -# out_file: output file to write the stdout -# err_file: output file to write the stderr -# command: command to execute -# cwd: working directory -def execute_cmd(log_file, err_log_file, command, cwd='./'): - print(f'In {os.path.abspath(cwd)}') - print(f'Command: {command}') - result = subprocess.run(f'{command}', shell=True, - cwd=cwd, - stdout=subprocess.PIPE, - stderr=subprocess.PIPE) - - with open(log_file, 'a') as f: - f.write(f'Command: {command}') - f.write(result.stdout.decode("utf8")) - with open(err_log_file, 'a') as f: - f.write(result.stderr.decode("utf8")) - - if result.returncode == 0: - print("Command executed successfully.\n") - else: - print(f"Error executing command. See {err_log_file}.\n") - exit(1) - - -def create_dir(dir_path): - os.makedirs(dir_path, exist_ok=True) - # Check if the directory was created - if not os.path.exists(dir_path): - print(f"Error: Could not create directory {dir_path}") - exit(1) - +from script.utilities import * def set_up_batch_header(script_file, job_name, job_dir, num_nodes): # Make dir @@ -76,31 +42,3 @@ def add_srun_cmd(num_tasks_per_node, command, script_file, echo=True, check_return_code=True): add_cmd(f"srun --ntasks-per-node={num_tasks_per_node} {command}", script_file, echo, check_return_code) - - -def grep_file(file, text_pattern): - # Return the lines that contain the pattern - with open(file, 'r') as f: - lines = f.readlines() - return [line for line in lines if text_pattern in line] - -# Parse a range string and return a list of integers -# e.g., "1-5" -> [1, 2, 3, 4, 5] -# e.g., "1,3,5" -> [1, 3, 5] -def parse_range(range_str): - if range_str is None: - return [None] - - if '-' in range_str: - start, end = range_str.split('-') - return list(range(int(start), int(end) + 1)) - else: - return [int(x) for x in range_str.split(',')] - -# Parse a list of floats from a string and return a list of floats -# e.g., "1.0,2.5,3.14" -> [1.0, 2.5, 3.14] -def parse_float_list(float_list_str): - if float_list_str is None: - return [None] - - return [float(x) for x in float_list_str.split(',')] \ No newline at end of file diff --git a/script/benchmark/clustering_utilities.py b/script/benchmark/clustering_utilities.py index 93777b4..2d280bb 100644 --- a/script/benchmark/clustering_utilities.py +++ b/script/benchmark/clustering_utilities.py @@ -6,7 +6,8 @@ from sklearn.metrics import adjusted_rand_score, adjusted_mutual_info_score import datetime import os - +import sys +from script.utilities import * # Modified from: # https://gist.github.com/lmcinnes/24ed5c22c80125be5133811d677eae7b @@ -75,162 +76,3 @@ def assign_singleton_cluster_to_noise_points(cluster_labels, noise_id_offset): print(f"Assigned singleton clusters to {cnt_noise} noise points") return new_labels - - -def show_time_now(): - now = datetime.datetime.now() - print(now.strftime("%Y-%m-%d %H:%M:%S")) - - -def find_files_in_dir(dir_path, ext=''): - if not os.path.isdir(dir_path): - return [dir_path] - - files = [] - for file in os.listdir(dir_path): - if len(ext) == 0 or file.endswith(ext): - files.append(os.path.join(dir_path, file)) - return files - - -# Read point(feature) data from file(s) -# If data_path is a directory, load all files in the directory. -# Return a dense numpy array of feature vectors. -# Each row is a feature vector. -# Input files can have point IDs in the first column. -# Point IDs are expected to be 0, 1, 2, ..., N-1, -# where N is the number of points. -# If the first column is not a point ID, set has_ids to False. -# If there are multiple files in a directory, the point IDs must be present. -def read_point_data(data_path, has_ids=True): - print(f"Loading data from {data_path}", flush=True) - files = find_files_in_dir(data_path) - - if len(files) == 0: - print(f"No files found in {data_path}") - exit(1) - - if len(files) > 1 and not has_ids: - print("Multiple files are provided," - "but has_ids is false.") - - points_table = {} - dimensions = 0 - for file in files: - for line in open(file, 'r'): - if len(line.strip()) == 0: - print(f"Empty line found in {file}") - continue - items = line.split() - if has_ids: - pid = int(items[0]) - items = items[1:] - else: - pid = len(points_table) - - if pid in points_table: - print(f"Duplicate ID in the feature file: {pid}") - exit(1) - - if dimensions == 0: - dimensions = len(items) - elif dimensions != len(items): - print( - f"Dimension mismatch in {file}: {len(items)} != {dimensions}") - exit(1) - - points_table[pid] = list(map(float, items)) - - print(f"Loaded {len(points_table)} items from {len(files)} files", flush=True) - - # numpy array of feature vectors - # if the IDs are not continuous, fill the missing IDs with -1 - max_id = max(points_table.keys()) - points = np.full((max_id + 1, len(points_table[0])), 0.0) - for id, feature in points_table.items(): - points[id] = feature - - return points - - -# Load cluster label or ground truth label data -# If data_path is a directory, load all files in the directory -# Return a dense numpy array of labels -# Unclustered points are filled with -1 -# -# There are two accepted file types: -# -# 1. Only cluster IDs: -# Each line contains a cluster ID. -# Point IDs are assumed to be 0, 1, 2, ..., N-1, where N is the number of points. -# This mode is not acceptable if there are multiple files in a directory. -# -# 2. Point IDs and cluster IDs: -# Contains two columns: -# the first column is for point IDs and the second column is for cluster IDs. -# -# Both File types can also contain comment lines, which must start from #. -def read_label_data(data_path): - print(f"Loading data from {data_path}", flush=True) - labels_dict = {} - files = [] - if os.path.isdir(data_path): - files = [os.path.join(data_path, f) for f in os.listdir(data_path)] - else: - files.append(data_path) - - if len(files) == 0: - print(f"No files found in {data_path}") - exit(1) - - contains_ids = False - for line in open(files[0], 'r'): - if len(line.split()) >= 2: - contains_ids = True - break - - if contains_ids: - print("Loading point IDs and labels", flush=True) - else: - print("Loading only labels", flush=True) - - if len(files) > 1 and not contains_ids: - print("Multiple files are provided," - "but no IDs are found in the first file.") - - max_id = 0 - for file in files: - for line in open(file, 'r'): - - # Skip comments - if line.startswith("#"): - continue - - if contains_ids: - items = line.split() - if len(items) < 2: - print(f"Invalid line: {line}") - exit(1) - pid = int(items[0]) - label = int(items[1]) - max_id = max(max_id, pid) - else: - pid = len(labels_dict) - max_id = pid - label = int(line) - - if pid in labels_dict: - print(f"Duplicate ID in the input files: {pid}") - exit(1) - labels_dict[pid] = label - - print(f"Loaded {len(labels_dict)} items from {len(files)} files", flush=True) - print(f"Max ID: {max_id}", flush=True) - - # numpy array of labels - # if the IDs are not continuous, fill the missing IDs with -1 - labels = np.full(max_id + 1, -1) - for id, label in labels_dict.items(): - labels[id] = label - - return labels diff --git a/script/benchmark/evaluate_clustering_quality.py b/script/benchmark/evaluate_clustering_quality.py index 6dc700a..0396516 100644 --- a/script/benchmark/evaluate_clustering_quality.py +++ b/script/benchmark/evaluate_clustering_quality.py @@ -37,7 +37,7 @@ import argparse import numpy as np from clustering_utilities import * - +from script.utilities import * def parse_options(): parser = argparse.ArgumentParser( @@ -63,8 +63,8 @@ def parse_options(): def main(): opt = parse_options() - cluster_labels = read_label_data(opt.cluster_path) - true_labels = read_label_data(opt.gt_path) + cluster_labels = read_cluster_labels(opt.cluster_path, "*") + true_labels = read_cluster_labels(opt.gt_path, "*") try: eval_clusters(cluster_labels, true_labels, opt.singleton_noise_points, opt.ignore_true_noise_points) diff --git a/script/benchmark/evaluate_correlation_knng_and_clusters.py b/script/benchmark/evaluate_correlation_knng_and_clusters.py index 7ddc7e4..392c13e 100644 --- a/script/benchmark/evaluate_correlation_knng_and_clusters.py +++ b/script/benchmark/evaluate_correlation_knng_and_clusters.py @@ -35,27 +35,8 @@ import glob from pathlib import Path from typing import Iterable, Iterator, Sequence - - -def resolve_files(paths: Sequence[str]) -> list[Path]: - """Resolve file, directory, and glob inputs to a unique sorted file list.""" - resolved: set[Path] = set() - for value in paths: - path = Path(value) - if path.is_dir(): - matches = set(path.glob("*.txt")) | set(path.glob("knng-*.txt")) - elif path.is_file(): - matches = {path} - else: - matches = {Path(match) for match in glob.glob(value)} - - resolved.update(match.resolve() for match in matches if match.is_file()) - - if not resolved: - raise FileNotFoundError(f"No files matched: {paths}") - - return sorted(resolved) - +from script.utilities import find_files_in_dir +from script.data_readers import read_cluster_labels def iter_knng_records(path: Path, contains_distance: bool) -> Iterator[tuple[int, list[int]]]: """Yield (source_id, neighbor_ids) pairs from one KNNG dump file.""" @@ -79,62 +60,16 @@ def iter_knng_records(path: Path, contains_distance: bool) -> Iterator[tuple[int yield src_id, neighbor_ids -def read_cluster_labels(cluster_path: str | Path) -> dict[int, int]: - """Read cluster labels from a file or directory into a point->cluster map.""" - labels: dict[int, int] = {} - files = resolve_files([str(cluster_path)]) - - if not files: - raise FileNotFoundError(f"No cluster files found in {cluster_path}") - - # Heuristic: determine whether the file format contains point IDs by scanning - # the first non-comment, non-empty line. - first_file = files[0] - contains_ids = False - with first_file.open("r") as handle: - for line in handle: - stripped = line.strip() - if not stripped or stripped.startswith("#"): - continue - if len(stripped.split()) >= 2: - contains_ids = True - break - - for file in files: - with file.open("r") as handle: - for line in handle: - stripped = line.strip() - if not stripped or stripped.startswith("#"): - continue - - items = stripped.split() - if contains_ids: - if len(items) < 2: - raise ValueError(f"{file}: invalid cluster row: {line.rstrip()}") - point_id, cluster_id = int(items[0]), int(items[1]) - else: - point_id = len(labels) - cluster_id = int(items[0]) - - if point_id in labels: - raise ValueError(f"Duplicate point ID {point_id} in cluster file {file}") - labels[point_id] = cluster_id - - if not labels: - raise ValueError(f"No cluster labels were loaded from {cluster_path}") - - return labels - def evaluate_knng_clusters(graph_path: str | Path, cluster_path: str | Path, k: int = 10, no_distance: bool = False): """Return per-point and aggregate same-cluster neighbor statistics.""" if k <= 0: raise ValueError("-k/--neighbors must be positive") - labels = read_cluster_labels(cluster_path) + labels = read_cluster_labels(cluster_path, "*") print(f"Loaded {len(labels)} cluster labels from {cluster_path}") - files = resolve_files([str(graph_path)]) + files = find_files_in_dir([str(graph_path)]) per_point: list[tuple[int, int, int, int, float]] = [] total_same = 0 @@ -179,36 +114,6 @@ def evaluate_knng_clusters(graph_path: str | Path, cluster_path: str | Path, k: "num_clustered_points": len(labels), } - -def _parse_k_values(raw_values: Sequence[Sequence[str]] | None) -> list[int]: - """Normalize one or more k values from CLI input (supports repeated or comma-separated lists).""" - values: list[int] = [] - raw_entries = raw_values or [] - - for group in raw_entries: - if isinstance(group, str): - group_tokens = [group] - else: - group_tokens = list(group) - - for token in group_tokens: - for item in str(token).split(","): - item = item.strip() - if not item: - continue - try: - value = int(item) - except ValueError as exc: - raise argparse.ArgumentTypeError(f"invalid neighbor count: {item!r}") from exc - if value <= 0: - raise argparse.ArgumentTypeError(f"neighbor count must be positive: {value}") - values.append(value) - - if not values: - return [10] - return values - - def _parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser( description=( @@ -253,7 +158,7 @@ def _parse_args() -> argparse.Namespace: help="Optional output file path for the summary table", ) args = parser.parse_args() - args.neighbors = _parse_k_values(args.neighbors or [[10]]) + args.neighbors = parse_range(args.neighbors or [[10]]) return args diff --git a/script/benchmark/run_hdbscan_bench.py b/script/benchmark/run_hdbscan_bench.py index 94b8a6c..3670d2b 100644 --- a/script/benchmark/run_hdbscan_bench.py +++ b/script/benchmark/run_hdbscan_bench.py @@ -18,7 +18,7 @@ from clustering_utilities import * from bench_utilities import * - +from script.utilities import * def parse_options(): parser = argparse.ArgumentParser( @@ -101,7 +101,7 @@ def run_hdbscan(points, min_cluster_size, min_samples, # Evaluate the clustering quality if gt_file: print('\nLoading ground truth data') - gt_labels = read_label_data(gt_file) + gt_labels = read_cluster_labels(gt_file, "*") print('\nEvaluating clustering quality') eval_clusters(clusters.labels_, gt_labels) @@ -152,7 +152,7 @@ def main(): opts = parse_options() - points = read_point_data(opts.point_data_path, opts.has_ids) + points = read_points(opts.point_data_path, "*", opts.has_ids) min_cluster_size_list = parse_range(opts.min_cluster_size_range) min_samples_list = parse_range(opts.min_samples_range) diff --git a/script/data_readers.py b/script/data_readers.py new file mode 100644 index 0000000..3a9393a --- /dev/null +++ b/script/data_readers.py @@ -0,0 +1,255 @@ +# Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS +# Project Developers. See the top-level COPYRIGHT file for details. + +from __future__ import annotations + +import math +from pathlib import Path + +import numpy as np +from .utilities import find_files_in_dir + +def read_points( + data_path: str | Path, + pattern: str | None = None, + has_ids: bool = True, +) -> np.ndarray: + """Load whitespace-separated feature vectors into a dense NumPy array. + + Each non-empty row contains feature values, optionally preceded by a point + ID. When IDs are present, missing IDs are represented by 0-based row numbers. + Multiple files require explicit IDs so records can be combined unambiguously. + """ + print(f"Loading data from {data_path}", flush=True) + files = find_files_in_dir([data_path], pattern) + if len(files) > 1 and not has_ids: + raise ValueError("Multiple point files require point IDs") + + points_table: dict[int, list[float]] = {} + dimensions: int | None = None + for file in files: + with file.open("r", encoding="utf-8") as input_file: + for line_number, line in enumerate(input_file, start=1): + items = line.split() + if not items: + continue + + if has_ids: + try: + point_id = int(items.pop(0)) + except ValueError as exc: + raise ValueError( + f"{file}:{line_number}: invalid point ID" + ) from exc + if point_id < 0: + raise ValueError(f"{file}:{line_number}: point IDs must be non-negative") + else: + point_id = len(points_table) + + if not items: + raise ValueError(f"{file}:{line_number}: feature vector is empty") + if point_id in points_table: + raise ValueError(f"{file}:{line_number}: duplicate point ID {point_id}") + + try: + feature = [float(value) for value in items] + except ValueError as exc: + raise ValueError(f"{file}:{line_number}: invalid feature value") from exc + + if dimensions is None: + dimensions = len(feature) + elif len(feature) != dimensions: + raise ValueError( + f"{file}:{line_number}: dimension mismatch: " + f"{len(feature)} != {dimensions}" + ) + points_table[point_id] = feature + + if not points_table or dimensions is None: + raise ValueError(f"No feature vectors were loaded from {data_path}") + + print(f"Loaded {len(points_table)} items from {len(files)} files", flush=True) + + points = np.zeros((max(points_table) + 1, dimensions), dtype=float) + for point_id, feature in points_table.items(): + points[point_id] = feature + + return points + + +def read_cluster_labels( + data_path: str | Path, + pattern: str | None = None, +) -> np.ndarray: + """Load cluster labels from one or more files into a dense integer array. + + Rows may contain only a cluster ID (point ID is the 0-based row number) or a + point ID followed by a cluster ID. Blank lines and comments are ignored. + Multiple files require explicit IDs. + """ + files = find_files_in_dir([data_path], pattern) + + print(f"Found {len(files)} files in {data_path}", flush=True) + first_data_row: list[str] | None = None + for file in files: + with file.open("r", encoding="utf-8") as input_file: + for line in input_file: + stripped = line.strip() + if stripped and not stripped.startswith("#"): + first_data_row = stripped.split() + break + if first_data_row is not None: + break + + if first_data_row is None: + raise ValueError(f"No cluster labels were found in {data_path}") + + contains_ids = len(first_data_row) >= 2 + + if contains_ids: + print("Loading point IDs and labels", flush=True) + else: + print("Loading only labels", flush=True) + + if len(files) > 1 and not contains_ids: + raise ValueError("Multiple cluster files require point IDs") + + labels_dict: dict[int, int] = {} + for file in files: + with file.open("r", encoding="utf-8") as input_file: + for line_number, line in enumerate(input_file, start=1): + stripped = line.strip() + if not stripped or stripped.startswith("#"): + continue + items = stripped.split() + + try: + if contains_ids: + if len(items) < 2: + raise ValueError("expected point ID and cluster ID") + point_id, label = int(items[0]), int(items[1]) + else: + point_id, label = len(labels_dict), int(items[0]) + except ValueError as exc: + raise ValueError(f"{file}:{line_number}: invalid cluster row") from exc + + if point_id < 0: + raise ValueError(f"{file}:{line_number}: point IDs must be non-negative") + if point_id in labels_dict: + raise ValueError(f"{file}:{line_number}: duplicate point ID {point_id}") + labels_dict[point_id] = label + + if not labels_dict: + raise ValueError(f"No cluster labels were found in {data_path}") + + print(f"Loaded {len(labels_dict)} items from {len(files)} files", flush=True) + max_id = max(labels_dict) + print(f"Max ID: {max_id}", flush=True) + + labels = np.full(max_id + 1, -1) + for point_id, label in labels_dict.items(): + labels[point_id] = label + + return labels + +Graph: TypeAlias = dict[int, list[tuple[int, float]]] + +def read_knng( + graph_path: str | Path, + knng_format: str, + path_pattern: str | None = None, +) -> Graph: + """Load a kNN graph as source IDs mapped to neighbor/distance pairs. + + ``knng_format`` selects ``N`` (neighbor IDs only), ``B`` (blocked ID rows + followed by their distance rows), or ``M`` (alternating ID and distance + rows). Append ``I`` to any format when each ID row starts with its source + point ID; otherwise source IDs are assigned sequentially across the files. + Distances are ``None`` for the ID-only format. With explicit source IDs, + the corresponding first distance is treated as the dummy source distance. + """ + format_code = knng_format.upper() + base_formats = [code for code in format_code if code in "NBM"] + if len(base_formats) != 1 or any(code not in "NBMI" for code in format_code): + raise ValueError(f"Unsupported kNNG format: {knng_format!r}") + base_format = base_formats[0] + contains_ids = "I" in format_code + + files = find_files_in_dir([graph_path], path_pattern) + graph: dict[int, list[tuple[int, float | None]]] = {} + next_source_id = 0 + + def parse_id_row(file: Path, line_number: int, line: str) -> tuple[int, list[int]]: + nonlocal next_source_id + try: + values = [int(value) for value in line.split()] + except ValueError as exc: + raise ValueError(f"{file}:{line_number}: invalid neighbor ID row") from exc + if contains_ids: + if not values: + raise ValueError(f"{file}:{line_number}: missing source point ID") + source_id, neighbor_ids = values[0], values[1:] + else: + source_id, neighbor_ids = next_source_id, values + next_source_id += 1 + if source_id < 0 or any(neighbor_id < 0 for neighbor_id in neighbor_ids): + raise ValueError(f"{file}:{line_number}: point IDs must be non-negative") + if source_id in graph: + raise ValueError(f"{file}:{line_number}: duplicate source point ID {source_id}") + return source_id, neighbor_ids + + def parse_distances(file: Path, line_number: int, line: str) -> list[float]: + try: + distances = [float(value) for value in line.split()] + except ValueError as exc: + raise ValueError(f"{file}:{line_number}: invalid distance row") from exc + if any(not math.isfinite(distance) or distance < 0 for distance in distances): + raise ValueError(f"{file}:{line_number}: distances must be finite and non-negative") + return distances + + for file in files: + with file.open("r", encoding="utf-8") as input_file: + lines = [ + (line_number, line.strip()) + for line_number, line in enumerate(input_file, start=1) + if line.strip() and not line.lstrip().startswith("#") + ] + + if base_format == "N": + for line_number, line in lines: + source_id, neighbor_ids = parse_id_row(file, line_number, line) + graph[source_id] = [(neighbor_id, None) for neighbor_id in neighbor_ids] + continue + + if len(lines) % 2 != 0: + raise ValueError(f"{file}: expected an even number of kNNG lines") + half = len(lines) // 2 if base_format == "B" else None + if base_format == "B": + id_rows = lines[:half] + distance_rows = lines[half:] + else: + id_rows = lines[::2] + distance_rows = lines[1::2] + + for (id_line_number, id_line), (distance_line_number, distance_line) in zip( + id_rows, distance_rows + ): + source_id, neighbor_ids = parse_id_row(file, id_line_number, id_line) + distances = parse_distances(file, distance_line_number, distance_line) + if contains_ids: + if len(distances) != len(neighbor_ids) + 1: + raise ValueError( + f"{file}:{distance_line_number}: expected a dummy distance and " + f"{len(neighbor_ids)} neighbor distances" + ) + distances = distances[1:] + elif len(distances) != len(neighbor_ids): + raise ValueError( + f"{file}:{distance_line_number}: {len(neighbor_ids)} neighbors but " + f"{len(distances)} distances" + ) + graph[source_id] = list(zip(neighbor_ids, distances)) + + if not graph: + raise ValueError(f"No kNNG rows were loaded from {graph_path}") + return graph diff --git a/script/etc/reassign_sequential_cluster_ids.py b/script/etc/reassign_sequential_cluster_ids.py index 12ed8f3..d0aaac8 100644 --- a/script/etc/reassign_sequential_cluster_ids.py +++ b/script/etc/reassign_sequential_cluster_ids.py @@ -28,7 +28,8 @@ import argparse import numpy as np from script.benchmark.clustering_utilities import * - +from script.utilities import * +from script.data_readers import * def main(): parser = argparse.ArgumentParser( @@ -43,7 +44,7 @@ def main(): help='Path to the output file') args = parser.parse_args() - cluster_ids = read_label_data(args.input_cluster_ids_path) + cluster_ids = read_cluster_labels(args.input_cluster_ids_path, "*") # Make new cluster IDs new_ids = dict() diff --git a/script/utilities.py b/script/utilities.py new file mode 100644 index 0000000..457b297 --- /dev/null +++ b/script/utilities.py @@ -0,0 +1,178 @@ +# Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS +# Project Developers. See the top-level COPYRIGHT file for details. + +from __future__ import annotations +import datetime +from collections.abc import Sequence +import glob +from pathlib import Path +import subprocess + +import numpy as np +from numpy.typing import ArrayLike, NDArray + + +def _as_vector_pair( + vec1: ArrayLike, vec2: ArrayLike +) -> tuple[NDArray[np.float64], NDArray[np.float64]]: + """Convert and validate two finite, non-empty one-dimensional vectors.""" + first = np.asarray(vec1, dtype=np.float64) + second = np.asarray(vec2, dtype=np.float64) + if first.ndim != 1 or second.ndim != 1: + raise ValueError("Vectors must be one-dimensional") + if first.size != second.size: + raise ValueError("Vectors must have the same length") + if first.size == 0: + raise ValueError("Vectors must not be empty") + if not np.isfinite(first).all() or not np.isfinite(second).all(): + raise ValueError("Vectors must contain only finite values") + return first, second + + +def show_time_now(): + now = datetime.datetime.now() + print(now.strftime("%Y-%m-%d %H:%M:%S")) + + +def find_files_in_dir( + paths: Sequence[str | Path], pattern: str | None = None +) -> list[Path]: + """Resolve files, directories, and glob patterns to unique sorted paths. + + ``pattern`` filters entries inside directory inputs. Nonexistent path inputs + are treated as glob patterns. Directory searches are not recursive. + """ + resolved: set[Path] = set() + for value in paths: + path = Path(value) + if path.is_dir(): + matches = path.glob(pattern or "*") + elif path.is_file(): + matches = (path,) + else: + # Nonexistent inputs may still be glob expressions, such as *.txt. + matches = (Path(match) for match in glob.glob(str(value))) + + resolved.update(match.resolve() for match in matches if match.is_file()) + + if not resolved: + formatted_paths = ", ".join(map(str, paths)) + raise FileNotFoundError(f"No files matched: {formatted_paths}") + + return sorted(resolved) + + +def create_dir(dir_path: str | Path) -> None: + """Create a directory and any missing parent directories if needed.""" + Path(dir_path).mkdir(parents=True, exist_ok=True) + + +def grep_file(file: str | Path, text_pattern: str) -> list[str]: + """Return the lines in ``file`` that contain the given substring.""" + with Path(file).open("r", encoding="utf-8") as input_file: + return [line for line in input_file if text_pattern in line] + + +def parse_range(range_str: str | None) -> list[int | None]: + """Parse an inclusive integer range or comma-separated integer list. + + Examples: ``"1-3"`` becomes ``[1, 2, 3]`` and ``"1,3"`` becomes + ``[1, 3]``. ``None`` is preserved as ``[None]`` for optional CLI values. + """ + if range_str is None: + return [None] + + value = range_str.strip() + if not value: + raise ValueError("range must not be empty") + + if "-" in value: + parts = value.split("-") + if len(parts) != 2: + raise ValueError(f"invalid integer range: {range_str!r}") + start, end = (int(part.strip()) for part in parts) + if start > end: + raise ValueError(f"range start must not exceed its end: {range_str!r}") + return list(range(start, end + 1)) + + return [int(part.strip()) for part in value.split(",")] + + +def parse_float_list(float_list_str: str | None) -> list[float | None]: + """Parse comma-separated floating-point values into a list. + + ``None`` is preserved as ``[None]`` for optional CLI values. + """ + if float_list_str is None: + return [None] + + value = float_list_str.strip() + if not value: + raise ValueError("float list must not be empty") + return [float(part.strip()) for part in value.split(",")] + + +def execute_cmd( + log_file: str | Path, + err_log_file: str | Path, + command: str, + cwd: str | Path = "./", +) -> None: + """Run a shell command, append its output to logs, and exit on failure. + + The command is interpreted by the system shell. Standard output and + standard error are appended to their respective log files. + """ + print(f"In {Path(cwd).resolve()}") + print(f"Command: {command}") + result = subprocess.run( + command, + shell=True, + cwd=cwd, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + text=True, + encoding="utf-8", + errors="replace", + check=False, + ) + + with Path(log_file).open("a", encoding="utf-8") as output_file: + output_file.write(f"Command: {command}\n") + output_file.write(result.stdout) + with Path(err_log_file).open("a", encoding="utf-8") as error_file: + error_file.write(result.stderr) + + if result.returncode == 0: + print("Command executed successfully.\n") + else: + print(f"Error executing command. See {err_log_file}.\n") + raise SystemExit(result.returncode) + + +# Distance (similarity) functions +# The smaller the distance value, the more similar the vectors are. +# Inputs are numpy arrays or lists of floats. +def l2_distance(vec1: ArrayLike, vec2: ArrayLike) -> float: + """Compute the Euclidean (L2) distance between two vectors.""" + first, second = _as_vector_pair(vec1, vec2) + return float(np.linalg.norm(first - second)) + + +def l1_distance(vec1: ArrayLike, vec2: ArrayLike) -> float: + """Compute the Manhattan (L1) distance between two vectors.""" + first, second = _as_vector_pair(vec1, vec2) + return float(np.sum(np.abs(first - second))) + + +def cosine_distance(vec1: ArrayLike, vec2: ArrayLike) -> float: + """Compute the Cosine distance between two vectors (1 - cosine similarity).""" + first, second = _as_vector_pair(vec1, vec2) + norm1 = np.linalg.norm(first) + norm2 = np.linalg.norm(second) + if norm1 == 0 or norm2 == 0: + raise ValueError("Vectors must not be zero vectors") + + cosine_similarity = np.dot(first / norm1, second / norm2) + cosine_similarity = np.clip(cosine_similarity, -1.0, 1.0) + return float(1.0 - cosine_similarity) From 77f9cfca298ad80c1b49e209ee1fd3fc452f2978 Mon Sep 17 00:00:00 2001 From: Keita Iwabuchi Date: Sun, 4 Oct 2026 12:58:14 -0700 Subject: [PATCH 37/37] (WIP) Implementing unsupervised clustering quality evaluation --- ...valuate_clustering_quality_unsupervised.py | 391 ++++++++++++++++++ 1 file changed, 391 insertions(+) create mode 100644 script/benchmark/evaluate_clustering_quality_unsupervised.py diff --git a/script/benchmark/evaluate_clustering_quality_unsupervised.py b/script/benchmark/evaluate_clustering_quality_unsupervised.py new file mode 100644 index 0000000..b9e0936 --- /dev/null +++ b/script/benchmark/evaluate_clustering_quality_unsupervised.py @@ -0,0 +1,391 @@ +#!/usr/bin/env python3 +"""Evaluate clustering separation without ground-truth labels. + +For every clustered point, traverse the directed kNNG to find reachable +same-cluster and different-cluster points with the smallest direct feature +distance from the source. The graph defines reachability only; stored edge +distances are ignored. + +kNNG records use two lines per point: + point_id neighbor_id_1 neighbor_id_2 ... + dummy_distance distance_1 distance_2 ... + +Cluster assignments contain two columns: point_id cluster_id. Cluster ID -1 +denotes noise and is excluded as a source or candidate. +""" + +from __future__ import annotations +from script.utilities import cosine_distance, l1_distance, l2_distance +from script.data_readers import read_cluster_labels, read_knng, read_points + +import argparse +import csv +import heapq +import math +from collections import deque +from pathlib import Path +import sys +from collections.abc import Callable +from typing import TypeAlias + +import numpy as np +from numpy.typing import ArrayLike + +PROJECT_ROOT = Path(__file__).resolve().parents[2] +if str(PROJECT_ROOT) not in sys.path: + sys.path.insert(0, str(PROJECT_ROOT)) + + +Graph: TypeAlias = dict[int, list[tuple[int, float | None]]] +DistanceFunction: TypeAlias = Callable[[ArrayLike, ArrayLike], float] +ResultRow: TypeAlias = dict[str, int | float | None] + + +def find_closest_cluster_point_dijkstra( + source_id: int, + labels: dict[int, int], + graph: Graph, + points: np.ndarray, + dist_func: DistanceFunction, +) -> dict[str, tuple[int, int, float] | None]: + """Find reachable cluster points nearest by direct feature distance. + + The graph determines reachability only. The frontier is prioritized by + direct source-to-node distance. Returned values are + ``(point_id, hops, direct_distance)``; equal candidate distances prefer + fewer distances. + """ + + if source_id < 0 or source_id >= len(points): + raise ValueError( + "Source point ID is outside the loaded feature vectors") + if source_id not in labels: + raise ValueError(f"Source point ID {source_id} has no cluster label") + + source_cluster = labels[source_id] + if source_cluster < 0: + # Skip noise points as source + return {"in": None, "out": None} + + # found contains: + # "in" or "out": (closest_point_id, hops, distance) or None + found: dict[str, tuple[int, int, float] | None] = {"in": None, "out": None} + + # Frontier priority queue contains tuples of (direct_distance, hops, point_id) + frontier = [(0.0, 0, source_id)] + + visited: set[int] = set() + visited.add(source_id) + + while frontier: + source_distance, hops, point_id = heapq.heappop(frontier) + + neighbor_hops = hops + 1 + for neighbor_id, _stored_distance in graph.get(point_id, ()): + if neighbor_id in visited: + continue + visited.add(neighbor_id) + + neighbor_distance = float( + dist_func(points[source_id], points[neighbor_id])) + if not math.isfinite(neighbor_distance) or neighbor_distance < 0: + raise ValueError( + "Distance function must return a finite non-negative value" + ) + heapq.heappush( + frontier, + (neighbor_distance, neighbor_hops, neighbor_id), + ) + + category = "in" if labels[neighbor_id] == source_cluster else "out" + current = found[category] + if current is None or neighbor_distance < current[2]: + found[category] = ( + neighbor_id, neighbor_hops, neighbor_distance) + + return found + + +def find_closest_cluster_point_bfs( + source_id: int, + labels: dict[int, int], + graph: Graph, + points: np.ndarray, + dist_func: DistanceFunction, +) -> dict[str, tuple[int, int, float] | None]: + """Find closest cluster points by BFS depth, ranking same-depth ties by distance. + + The first reachable graph layer containing a point in each category is + considered. Within that layer, direct feature distance selects the point. + """ + if source_id < 0 or source_id >= len(points): + raise ValueError( + "Source point ID is outside the loaded feature vectors") + if source_id not in labels: + raise ValueError(f"Source point ID {source_id} has no cluster label") + + source_cluster = labels[source_id] + if source_cluster < 0: + return {"in": None, "out": None} + + found: dict[str, tuple[int, int, float] | None] = {"in": None, "out": None} + visited = {source_id} + frontier = deque([source_id]) + hops = 0 + + while frontier and any(result is None for result in found.values()): + next_frontier: list[int] = [] + layer_matches: dict[str, tuple[int, float] + | None] = {"in": None, "out": None} + + while frontier: + point_id = frontier.popleft() + if point_id != source_id and labels.get(point_id, -1) >= 0: + category = "in" if labels[point_id] == source_cluster else "out" + if found[category] is None: + distance = float( + dist_func(points[source_id], points[point_id])) + if not math.isfinite(distance) or distance < 0: + raise ValueError( + "Distance function must return a finite non-negative value" + ) + current = layer_matches[category] + if current is None or (distance, point_id) < current: + layer_matches[category] = (distance, point_id) + + for neighbor_id, _stored_distance in graph.get(point_id, ()): + if neighbor_id in visited: + continue + if neighbor_id < 0 or neighbor_id >= len(points): + raise ValueError( + "Neighbor point ID is out of the expected range") + visited.add(neighbor_id) + next_frontier.append(neighbor_id) + + for category, match in layer_matches.items(): + if match is not None: + distance, point_id = match + found[category] = (point_id, hops, distance) + + frontier = deque(next_frontier) + hops += 1 + + return found + + +def compute_score(distance: float, epsilon: float) -> float: + """Return the base-10 log score for a non-negative distance.""" + if not math.isfinite(distance) or distance < 0: + raise ValueError("distance must be finite and non-negative") + if not math.isfinite(epsilon) or epsilon <= 0: + raise ValueError("epsilon must be a finite positive number") + return math.log10(distance + epsilon) + + +def find_closest_reachable_distances( + points: np.ndarray, + graph: Graph, + dist_func: DistanceFunction, + labels: dict[int, int], + algorithm: str = "bfs", +) -> list[ResultRow]: + """Find closest reachable in- and out-of-cluster distances per point.""" + results: list[ResultRow] = [] + for source_id in sorted(labels): + cluster_id = labels[source_id] + if cluster_id < 0: + continue + if algorithm == "bfs": + closest = find_closest_cluster_point_bfs( + source_id, labels, graph, points, dist_func + ) + elif algorithm == "dijkstra": + closest = find_closest_cluster_point_dijkstra( + source_id, labels, graph, points, dist_func + ) + else: + raise ValueError(f"Unsupported algorithm: {algorithm}") + row: ResultRow = { + "point_id": source_id, + "cluster_id": cluster_id, + } + for category in ("in", "out"): + match = closest[category] + row[f"closest_{category}_point_id"] = match[0] if match else None + row[f"closest_{category}_hops"] = match[1] if match else None + row[f"closest_{category}_distance"] = match[2] if match else None + results.append(row) + return results + + +# def write_csv(path: str | Path, results: list[ResultRow]) -> None: +# """Write per-point results as CSV.""" +# fields = [ +# "point_id", +# "cluster_id", +# "closest_in_point_id", +# "closest_in_hops", +# "closest_in_distance", +# "closest_out_point_id", +# "closest_out_hops", +# "closest_out_distance", +# ] +# output_path = Path(path) +# output_path.parent.mkdir(parents=True, exist_ok=True) +# with output_path.open("w", encoding="utf-8", newline="") as output_file: +# writer = csv.DictWriter(output_file, fieldnames=fields) +# writer.writeheader() +# writer.writerows(results) + + +def write_txt(path: str | Path, results: list[ResultRow]) -> None: + """Write per-point results as plain text.""" + output_path = Path(path) + output_path.parent.mkdir(parents=True, exist_ok=True) + # Dump keys in the first line + keys = list(results[0].keys()) if results else [] + with output_path.open("w", encoding="utf-8") as output_file: + output_file.write("\t".join(keys) + "\n") + for row in results: + output_file.write("\t".join(str(row[key]) for key in keys) + "\n") + + +def compute_score_hist(results: list[ResultRow], epsilon: float, num_bins: int) -> dict[str, list[float]]: + """Compute histograms of log scores for in-cluster and out-of-cluster distances.""" + hist = {"in": [0.0] * num_bins, "out": [0.0] * num_bins} + + # First, compute the min and max score + min_scores = {"in": float("inf"), "out": float("inf")} + max_scores = {"in": float("-inf"), "out": float("-inf")} + for row in results: + if row["closest_in_distance"] is not None: + score = compute_score(float(row["closest_in_distance"]), epsilon) + min_scores["in"] = min(min_scores["in"], score) + max_scores["in"] = max(max_scores["in"], score) + if row["closest_out_distance"] is not None: + score = compute_score(float(row["closest_out_distance"]), epsilon) + min_scores["out"] = min(min_scores["out"], score) + max_scores["out"] = max(max_scores["out"], score) + min_score = min(min_scores["in"], min_scores["out"]) + max_score = max(max_scores["in"], max_scores["out"]) + + bin_size = max((max_score - min_score) / num_bins, 1e-12) + + for row in results: + if row["closest_in_distance"] is not None: + score = compute_score(float(row["closest_in_distance"]), epsilon) + bin_index = min(int((score - min_score) / bin_size), num_bins - 1) + hist["in"][bin_index] += 1 + if row["closest_out_distance"] is not None: + score = compute_score(float(row["closest_out_distance"]), epsilon) + bin_index = min(int((score - min_score) / bin_size), num_bins - 1) + hist["out"][bin_index] += 1 + return hist + + +def print_summary(results: list[ResultRow]) -> None: + """Print counts and means of distances by category, and the global score.""" + print(f"Evaluated {len(results)} clustered points") + for category, description in (("in", "in-cluster"), ("out", "out-of-cluster")): + values = [ + float(row[f"closest_{category}_distance"]) + for row in results + if row[f"closest_{category}_distance"] is not None + ] + if values: + print( + f"Closest {description}: {len(values)} points; " + f"mean distance = {sum(values) / len(values):.6f}" + ) + else: + print(f"Closest {description}: no reachable candidates") + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("-p", "--points", required=True, + help="Feature vectors file or directory") + parser.add_argument( + "-I", "--ids", action="store_true", + help="Feature-vector rows include point IDs", + ) + parser.add_argument("-g", "--graph", required=True, + help="kNNG file or directory of files") + parser.add_argument( + "--knng-format", + choices=("M", "MI", "B", "BI"), + help="Distance-bearing kNNG format (default: MI)", + required=True, + ) + parser.add_argument( + "-f", "--dist_func", choices=("l2", "l1", "cosine"), + help="Distance function to use", + required=True, + ) + parser.add_argument("-c", "--cluster", required=True, + help="Point ID and cluster ID file") + parser.add_argument( + "-e", "--epsilon", type=float, default=0.001, + help="Positive value added before taking log10 (default: 0.001)", + ) + parser.add_argument( + "-a", "--algorithm", action="store", + help="Traversal algorithm used (BFS or Dijkstra)", choices=("bfs", "dijkstra"), + default="bfs", + ) + parser.add_argument( + "-v", "--verbose", action="store_true", + help="Enable verbose output", + ) + parser.add_argument( + "-o", "--output", help="Optional per-point CSV output path") + return parser.parse_args() + + +def main() -> None: + args = parse_args() + try: + print("Reading points...") + points = read_points(args.points, pattern=None, has_ids=args.ids) + print(f"Loaded {len(points)} points") + + print("Reading graph...") + graph = read_knng(args.graph, args.knng_format) + print(f"Loaded {len(graph)} graph points") + + print("Reading cluster labels...") + labels_array = read_cluster_labels(args.cluster) + labels = { + point_id: int(label) + for point_id, label in enumerate(labels_array) + if label >= 0 + } + print(f"Loaded {len(labels)} assigned cluster labels") + + distance_functions: dict[str, DistanceFunction] = { + "l2": l2_distance, + "l1": l1_distance, + "cosine": cosine_distance, + } + results = find_closest_reachable_distances( + points, graph, distance_functions[args.dist_func], labels, algorithm=args.algorithm + ) + + print_summary(results) + + hist = compute_score_hist(results, epsilon=args.epsilon, num_bins=100) + print("Histogram of closest reachable distances:") + for category, counts in hist.items(): + print(f"{category}: {counts}") + + if args.output: + file_name = args.output + "_scores.txt" + print(f"Writing per-point results to {file_name}...") + write_txt(file_name, results) + print(f"Per-point results written to {file_name}") + except (OSError, ValueError) as exc: + raise SystemExit(f"Error: {exc}") from exc + + +if __name__ == "__main__": + main()