diff --git a/script/benchmark/clustering_utilities.py b/script/benchmark/clustering_utilities.py index e72f8be..6cb80ae 100644 --- a/script/benchmark/clustering_utilities.py +++ b/script/benchmark/clustering_utilities.py @@ -10,23 +10,30 @@ # 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): + 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): 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 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}") - + print(f"Ground truth cluster coverage (%): {pct_clustered_gt * 100:.2f}%") + + if ignore_true_noise_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)} 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}") + + 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") @@ -34,16 +41,17 @@ 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") + 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}") @@ -58,6 +66,9 @@ def assign_singleton_cluster_to_noise_points(cluster_labels, noise_id_offset): if label == -1: cnt_noise += 1 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/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/extract_clams_clustering_qualities.sh b/script/benchmark/extract_clams_clustering_qualities.sh new file mode 100755 index 0000000..74ff9d9 --- /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,Found Clusters,MST clustering,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"; 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 } + /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, mst_clustering, 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 19ee390..66696fd 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,18 +161,25 @@ 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(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 + +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, 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 "" @@ -188,6 +197,66 @@ 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) + finish_stage(job_script) + + if ground_truth_path: + 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, + 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" + 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) + finish_stage(job_script) + + 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, + num_tasks_per_node, verbose, singleton_cluster_to_noise_points) + finish_stage(job_script) + + job_script.write("echo \"\" \n") + # Function to generate the batch script for running a benchmark # @@ -220,25 +289,23 @@ 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: - 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 '' 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\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,104 +315,42 @@ 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) + finish_stage(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}" add_cmd(conv2coredist_cmd, job_script) + finish_stage(job_script) 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) + finish_stage(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") - 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: - 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") - # 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, - 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 + 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}") @@ -384,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/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 6b4f929..316268a 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. @@ -27,8 +30,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,54 +76,11 @@ 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); - 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); @@ -134,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; @@ -159,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; @@ -170,8 +128,15 @@ 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 + // 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); @@ -190,11 +155,12 @@ 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 } 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; @@ -210,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: {}", diff --git a/src/common.hpp b/src/common.hpp index b75efc6..705e429 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; @@ -77,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) { @@ -159,4 +165,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/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/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..35178a4 --- /dev/null +++ b/src/tools/evaluate_correlation_knng_and_clusters.cpp @@ -0,0 +1,162 @@ +// Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS +// Project Developers. See the top-level COPYRIGHT file for details. + +// 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 + +#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; +}