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()
 
 #
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..7ddc7e4
--- /dev/null
+++ b/script/benchmark/evaluate_correlation_knng_and_clusters.py
@@ -0,0 +1,305 @@
+#!/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.
+    -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:
+        <source_id> <neighbor_id_1> <neighbor_id_2> ...
+        <dummy_distance> <distance_1> <distance_2> ...
+    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.
+    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)
+            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
+
+        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()
+
diff --git a/script/benchmark/run_clams_bench.py b/script/benchmark/run_clams_bench.py
index 4aec0eb..66b33a3 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',
@@ -99,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',
@@ -137,7 +136,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',
@@ -145,9 +144,12 @@ 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/tools/copy_pm_datastore',
+                        help='Path to the PM datastore copy executable.')
 
     # Etc
     parser.add_argument('-v', '--verbose', action='store_true',
@@ -162,20 +164,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,
@@ -185,74 +190,67 @@ 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,
-                       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)
 
-    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 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:
-        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"
-        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)
+        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")
 
@@ -268,18 +266,21 @@ 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,
-                                  clustering_exe, distributed_hdbscan,
+                                  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'
 
@@ -296,7 +297,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 ''
@@ -306,15 +307,15 @@ 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")
             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)
@@ -322,7 +323,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}"
@@ -332,7 +333,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)
@@ -342,10 +343,10 @@ 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)
+                    noise_point_clustering_exe, cluster_noise_points)
             try_no += 1
 
         finish_stage(job_script)
@@ -390,15 +391,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:
@@ -420,21 +412,22 @@ 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,
                                                amst_approx_bound_list,
-                                               clustering_exe,
-                                               opts.distributed_hdbscan,
+                                               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 316268a..8279046 100644
--- a/src/clustering/cluster_noise_points.cpp
+++ b/src/clustering/cluster_noise_points.cpp
@@ -22,14 +22,19 @@
 #include <spdlog/spdlog.h>
 #include <spdlog/stopwatch.h>
 #include <boost/unordered/unordered_flat_map.hpp>
+#include <boost/unordered/unordered_flat_set.hpp>
 
 #include "../common.hpp"
+#include "../details/multithread_adjacency_list.hpp"
 
 using namespace clams;
 
 template <typename K, typename V>
 using map_t = boost::unordered::unordered_flat_map<K, V>;
 
+template <typename K>
+using set_t = boost::unordered::unordered_flat_set<K>;
+
 struct option {
   std::filesystem::path mst_edges_path;
   bool                  metall_mst{false};
@@ -80,7 +85,7 @@ int main(int argc, char* argv[]) {
   option opt;
   parse_option(argc, argv, opt);
 
-  map_t<id_t, std::vector<id_t>> mst_graph;
+  multithread_adjacency_list<id_t, std::pair<id_t, distance_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);
@@ -94,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]);
-      mst_graph[edge.ids[1]].push_back(edge.ids[0]);
+    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]);
-      mst_graph[edge.ids[1]].push_back(edge.ids[0]);
+    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});
     }
   }
 
@@ -116,7 +125,8 @@ int main(int argc, char* argv[]) {
   }
 
   map_t<id_t, id_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());
 
@@ -134,9 +144,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);
@@ -145,32 +155,46 @@ int main(int argc, char* argv[]) {
     }
     ++n_noise_points;
 
-    std::deque<id_t>  bfs_queue;
-    map_t<id_t, bool> visited;
-    bfs_queue.push_back(point_id);
-    visited[point_id]  = true;
-    bool found_cluster = false;
+    std::vector<std::pair<id_t, distance_t>> bfs_front;
 
-    while (!bfs_queue.empty() && !found_cluster) {
-      const auto current_point_id = bfs_queue.front();
-      bfs_queue.pop_front();
+    set_t<id_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<std::pair<id_t, distance_t>> next;
+      for (const auto& [pid, _] : bfs_front) {
+        // Traverse the neighbors of the current point in the MST
+        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
+          }
+          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);
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 <unistd.h>
+#include <algorithm>
+#include <cstdlib>
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+#include <limits>
+#include <queue>
+#include <stack>
+#include <string>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+#include <spdlog/spdlog.h>
+#include <spdlog/stopwatch.h>
+#include <boost/unordered/unordered_flat_map.hpp>
+#include <metall/metall.hpp>
+
+#include "../common.hpp"
+#include "../details/multithread_adjacency_list.hpp"
+
+using namespace clams;
+
+template <typename K, typename V>
+using map_t   = boost::unordered::unordered_flat_map<K, V>;
+using graph_t = multithread_adjacency_list<id_t, std::pair<id_t, distance_t>>;
+
+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
+      << "<<Usage>>\n"
+         "Required arguments:\n"
+         "  -m <path> path to a file or a directory that contains input MST.\n"
+         "  -c <path> path to input cluster IDs.\n"
+         "  -o <path> 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<weighted_edge_list_t>(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<id_t, id_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<distance_t, id_t, id_t>;
+  std::priority_queue<heap_edge_t, std::vector<heap_edge_t>, 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/clustering/run_distributed_hdbscan_clustering.cpp b/src/clustering/run_distributed_hdbscan_clustering.cpp
index 0b14b5c..5d16cea 100644
--- a/src/clustering/run_distributed_hdbscan_clustering.cpp
+++ b/src/clustering/run_distributed_hdbscan_clustering.cpp
@@ -1195,7 +1195,7 @@ int main(int argc, char *argv[]) {
       spdlog::warn("Unable to open label file: {}", label_file_path.string());
     }
 
-    label_ofs << "point_id\tcluster_id\n";
+    label_ofs << "# point_id\tcluster_id\n";
     for (const auto &[point_id, cluster_id] : point_to_cluster_id_map) {
       label_ofs << point_id << "\t" << cluster_id << "\n";
     }
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 <typename cluster_id_table_t>
-void read_cluster_ids(const std::filesystem::path &input_path,
-                      cluster_id_table_t          &point_cluster_map) {
+void read_cluster_ids(const std::vector<std::filesystem::path> &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;
+    }
   }
 }
 
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 <cassert>
-#include <mutex>
 #include <algorithm>
+#include <cassert>
 #include <cstddef>
 #include <functional>
+#include <mutex>
 
-#define CLAMS_USE_STL_CONTAINERS_IN_ADJLIST 0
-#if CLAMS_USE_STL_CONTAINERS_IN_ADJLIST
-#include <vector>
-#include <unordered_map>
-#include <scoped_allocator>
-#else
 #include <metall/container/scoped_allocator.hpp>
-#include <metall/container/unordered_map.hpp>
+#include <metall/container/unordered_node_map.hpp>
 #include <metall/container/vector.hpp>
-#endif
 
 #include <metall/utility/hash.hpp>
 #include <metall/utility/mutex.hpp>
@@ -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 <typename _key_type, typename _value_type,
           typename _base_allocator_type = std::allocator<std::byte>>
 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<std::pair<const key_type, list_type>>>;
-  using key_table_type =
-      container::unordered_map<key_type, list_type, metall::utility::hash<>,
-                               std::equal_to<key_type>,
-                               key_table_allocator_type>;
+  using key_table_type = container::unordered_node_map<
+      key_type, list_type, metall::utility::hash<>, std::equal_to<key_type>,
+      key_table_allocator_type>;
 
   using bank_table_allocator_type =
       container::scoped_allocator_adaptor<other_allocator_type<key_table_type>>;
@@ -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<k_num_banks>(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
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 <utility>
 #include <vector>
 
+#include <boost/unordered/unordered_flat_set.hpp>
+
 #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<std::filesystem::path> &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<std::filesystem::path> &knng_files,
       }
 
       std::vector<distance_t> 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<std::filesystem::path>              &knng_files,
+    const boost::unordered::unordered_flat_set<id_t>       &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<id_t> ids;
       {
+        std::getline(ifs, line);
+        std::istringstream iss(line);
+        id_t               buf;
+        while (iss >> buf) {
+          ids.push_back(buf);
+        }
+      }
+
+      std::vector<distance_t> 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<std::filesystem::path> &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]));
       }
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<typename neo_dnnd_t::point_type,
                                                dist_t>(opt.distance_name),
         comm, opt.verbose);
+
+    ygm::utility::timer read_points_timer;
     comm.cout0() << "\n<<Read Points>>" << 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<<kNNG Construction>>" << std::endl;
     auto knng = neo_dnnd.build(opt.index_k, opt.r, opt.delta,
                                opt.replicate_rate, opt.batch_size);
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 <algorithm>
 #include <cstdlib>
 #include <filesystem>
-#include <iostream>
-#include <vector>
 #include <fstream>
+#include <iostream>
 #include <string>
 #include <utility>
+#include <vector>
 
 #include <metall/metall.hpp>
 
@@ -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 <pm_knng_path> -o <output_path>" << 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<id_t, id_t>;
 // 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;