diff --git a/script/benchmark/bench_utilities.py b/script/benchmark/bench_utilities.py
index 4d16045..4435f02 100644
--- a/script/benchmark/bench_utilities.py
+++ b/script/benchmark/bench_utilities.py
@@ -1,44 +1,10 @@
 # Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS
 # Project Developers. See the top-level COPYRIGHT file for details.
 
-
 import subprocess
 import os
 
-
-# Execute a shell command
-# out_file: output file to write the stdout
-# err_file: output file to write the stderr
-# command: command to execute
-# cwd: working directory
-def execute_cmd(log_file, err_log_file, command, cwd='./'):
-    print(f'In {os.path.abspath(cwd)}')
-    print(f'Command: {command}')
-    result = subprocess.run(f'{command}', shell=True,
-                            cwd=cwd,
-                            stdout=subprocess.PIPE,
-                            stderr=subprocess.PIPE)
-
-    with open(log_file, 'a') as f:
-        f.write(f'Command: {command}')
-        f.write(result.stdout.decode("utf8"))
-    with open(err_log_file, 'a') as f:
-        f.write(result.stderr.decode("utf8"))
-
-    if result.returncode == 0:
-        print("Command executed successfully.\n")
-    else:
-        print(f"Error executing command. See {err_log_file}.\n")
-        exit(1)
-
-
-def create_dir(dir_path):
-    os.makedirs(dir_path, exist_ok=True)
-    # Check if the directory was created
-    if not os.path.exists(dir_path):
-        print(f"Error: Could not create directory {dir_path}")
-        exit(1)
-
+from script.utilities import *
 
 def set_up_batch_header(script_file, job_name, job_dir, num_nodes):
     # Make dir
@@ -76,31 +42,3 @@ def add_srun_cmd(num_tasks_per_node, command, script_file, echo=True,
                  check_return_code=True):
     add_cmd(f"srun --ntasks-per-node={num_tasks_per_node} {command}",
             script_file, echo, check_return_code)
-
-
-def grep_file(file, text_pattern):
-    # Return the lines that contain the pattern
-    with open(file, 'r') as f:
-        lines = f.readlines()
-        return [line for line in lines if text_pattern in line]
-
-# Parse a range string and return a list of integers
-# e.g., "1-5" -> [1, 2, 3, 4, 5]
-# e.g., "1,3,5" -> [1, 3, 5]
-def parse_range(range_str):
-    if range_str is None:
-        return [None]
-
-    if '-' in range_str:
-        start, end = range_str.split('-')
-        return list(range(int(start), int(end) + 1))
-    else:
-        return [int(x) for x in range_str.split(',')]
-
-# Parse a list of floats from a string and return a list of floats
-# e.g., "1.0,2.5,3.14" -> [1.0, 2.5, 3.14]
-def parse_float_list(float_list_str):
-    if float_list_str is None:
-        return [None]
-
-    return [float(x) for x in float_list_str.split(',')]
\ No newline at end of file
diff --git a/script/benchmark/clustering_utilities.py b/script/benchmark/clustering_utilities.py
index 81646e9..2d280bb 100644
--- a/script/benchmark/clustering_utilities.py
+++ b/script/benchmark/clustering_utilities.py
@@ -6,7 +6,8 @@
 from sklearn.metrics import adjusted_rand_score, adjusted_mutual_info_score
 import datetime
 import os
-
+import sys
+from script.utilities import *
 
 # Modified from:
 # https://gist.github.com/lmcinnes/24ed5c22c80125be5133811d677eae7b
@@ -75,162 +76,3 @@ def assign_singleton_cluster_to_noise_points(cluster_labels, noise_id_offset):
     print(f"Assigned singleton clusters to {cnt_noise} noise points")
 
     return new_labels
-
-
-def show_time_now():
-    now = datetime.datetime.now()
-    print(now.strftime("%Y-%m-%d %H:%M:%S"))
-
-
-def find_files_in_dir(dir_path, ext=''):
-    if not os.path.isdir(dir_path):
-        return [dir_path]
-
-    files = []
-    for file in os.listdir(dir_path):
-        if len(ext) == 0 or file.endswith(ext):
-            files.append(os.path.join(dir_path, file))
-    return files
-
-
-# Read point(feature) data from file(s)
-# If data_path is a directory, load all files in the directory.
-# Return a dense numpy array of feature vectors.
-# Each row is a feature vector.
-# Input files can have point IDs in the first column.
-# Point IDs are expected to be 0, 1, 2, ..., N-1,
-# where N is the number of points.
-# If the first column is not a point ID, set has_ids to False.
-# If there are multiple files in a directory, the point IDs must be present.
-def read_point_data(data_path, has_ids=True):
-    print(f"Loading data from {data_path}", flush=True)
-    files = find_files_in_dir(data_path)
-
-    if len(files) == 0:
-        print(f"No files found in {data_path}")
-        exit(1)
-
-    if len(files) > 1 and not has_ids:
-        print("Multiple files are provided,"
-              "but has_ids is false.")
-
-    points_table = {}
-    dimensions = 0
-    for file in files:
-        for line in open(file, 'r'):
-            if len(line.strip()) == 0:
-                print(f"Empty line found in {file}")
-                continue
-            items = line.split()
-            if has_ids:
-                pid = int(items[0])
-                items = items[1:]
-            else:
-                pid = len(points_table)
-
-            if pid in points_table:
-                print(f"Duplicate ID in the feature file: {pid}")
-                exit(1)
-
-            if dimensions == 0:
-                dimensions = len(items)
-            elif dimensions != len(items):
-                print(
-                    f"Dimension mismatch in {file}: {len(items)} != {dimensions}")
-                exit(1)
-
-            points_table[pid] = list(map(float, items))
-
-    print(f"Loaded {len(points_table)} items from {len(files)} files", flush=True)
-
-    # numpy array of feature vectors
-    # if the IDs are not continuous, fill the missing IDs with -1
-    max_id = max(points_table.keys())
-    points = np.full((max_id + 1, len(points_table[0])), 0.0)
-    for id, feature in points_table.items():
-        points[id] = feature
-
-    return points
-
-
-# Load cluster label or ground truth label data
-# If data_path is a directory, load all files in the directory
-# Return a dense numpy array of labels
-# Unclustered points are filled with -1
-#
-# There are two accepted file types:
-#
-# 1. Only cluster IDs:
-# Each line contains a cluster ID.
-# Point IDs are assumed to be 0, 1, 2, ..., N-1, where N is the number of points.
-# This mode is not acceptable if there are multiple files in a directory.
-#
-# 2. Point IDs and cluster IDs:
-# Contains two columns:
-# the first column is for point IDs and the second column is for cluster IDs.
-#
-# Both File types can also contain comment lines, which must start from #.
-def read_label_data(data_path):
-    print(f"Loading data from {data_path}", flush=True)
-    labels_dict = {}
-    files = []
-    if os.path.isdir(data_path):
-        files = [os.path.join(data_path, f) for f in os.listdir(data_path)]
-    else:
-        files.append(data_path)
-
-    if len(files) == 0:
-        print(f"No files found in {data_path}")
-        exit(1)
-
-    contains_ids = False
-    for line in open(files[0], 'r'):
-        if len(line.split()) >= 2:
-            contains_ids = True
-            break
-
-    if contains_ids:
-        print("Loading point IDs and labels", flush=True)
-    else:
-        print("Loading only labels", flush=True)
-
-    if len(files) > 1 and not contains_ids:
-        print("Multiple files are provided,"
-              "but no IDs are found in the first file.")
-
-    max_id = 0
-    for file in files:
-        for line in open(file, 'r'):
-
-            # Skip comments
-            if line.startswith("#"):
-                continue
-
-            if contains_ids:
-                items = line.split()
-                if len(items) < 2:
-                    print(f"Invalid line: {line}")
-                    exit(1)
-                pid = int(items[0])
-                label = int(items[1])
-                max_id = max(max_id, pid)
-            else:
-                pid = len(labels_dict)
-                max_id = pid
-                label = int(line)
-
-            if pid in labels_dict:
-                print(f"Duplicate ID in the ground truth file: {pid}")
-                exit(1)
-            labels_dict[pid] = label
-
-    print(f"Loaded {len(labels_dict)} items from {len(files)} files", flush=True)
-    print(f"Max ID: {max_id}", flush=True)
-
-    # numpy array of labels
-    # if the IDs are not continuous, fill the missing IDs with -1
-    labels = np.full(max_id + 1, -1)
-    for id, label in labels_dict.items():
-        labels[id] = label
-
-    return labels
diff --git a/script/benchmark/evaluate_clustering_quality.py b/script/benchmark/evaluate_clustering_quality.py
index 6dc700a..0396516 100644
--- a/script/benchmark/evaluate_clustering_quality.py
+++ b/script/benchmark/evaluate_clustering_quality.py
@@ -37,7 +37,7 @@
 import argparse
 import numpy as np
 from clustering_utilities import *
-
+from script.utilities import *
 
 def parse_options():
     parser = argparse.ArgumentParser(
@@ -63,8 +63,8 @@ def parse_options():
 def main():
     opt = parse_options()
 
-    cluster_labels = read_label_data(opt.cluster_path)
-    true_labels = read_label_data(opt.gt_path)
+    cluster_labels = read_cluster_labels(opt.cluster_path, "*")
+    true_labels = read_cluster_labels(opt.gt_path, "*")
 
     try:
         eval_clusters(cluster_labels, true_labels, opt.singleton_noise_points, opt.ignore_true_noise_points)
diff --git a/script/benchmark/evaluate_clustering_quality_unsupervised.py b/script/benchmark/evaluate_clustering_quality_unsupervised.py
new file mode 100644
index 0000000..b9e0936
--- /dev/null
+++ b/script/benchmark/evaluate_clustering_quality_unsupervised.py
@@ -0,0 +1,391 @@
+#!/usr/bin/env python3
+"""Evaluate clustering separation without ground-truth labels.
+
+For every clustered point, traverse the directed kNNG to find reachable
+same-cluster and different-cluster points with the smallest direct feature
+distance from the source. The graph defines reachability only; stored edge
+distances are ignored.
+
+kNNG records use two lines per point:
+	point_id neighbor_id_1 neighbor_id_2 ...
+	dummy_distance distance_1 distance_2 ...
+
+Cluster assignments contain two columns: point_id cluster_id. Cluster ID -1
+denotes noise and is excluded as a source or candidate.
+"""
+
+from __future__ import annotations
+from script.utilities import cosine_distance, l1_distance, l2_distance
+from script.data_readers import read_cluster_labels, read_knng, read_points
+
+import argparse
+import csv
+import heapq
+import math
+from collections import deque
+from pathlib import Path
+import sys
+from collections.abc import Callable
+from typing import TypeAlias
+
+import numpy as np
+from numpy.typing import ArrayLike
+
+PROJECT_ROOT = Path(__file__).resolve().parents[2]
+if str(PROJECT_ROOT) not in sys.path:
+    sys.path.insert(0, str(PROJECT_ROOT))
+
+
+Graph: TypeAlias = dict[int, list[tuple[int, float | None]]]
+DistanceFunction: TypeAlias = Callable[[ArrayLike, ArrayLike], float]
+ResultRow: TypeAlias = dict[str, int | float | None]
+
+
+def find_closest_cluster_point_dijkstra(
+    source_id: int,
+    labels: dict[int, int],
+    graph: Graph,
+    points: np.ndarray,
+    dist_func: DistanceFunction,
+) -> dict[str, tuple[int, int, float] | None]:
+    """Find reachable cluster points nearest by direct feature distance.
+
+    The graph determines reachability only. The frontier is prioritized by
+    direct source-to-node distance. Returned values are
+    ``(point_id, hops, direct_distance)``; equal candidate distances prefer
+    fewer distances.
+    """
+
+    if source_id < 0 or source_id >= len(points):
+        raise ValueError(
+            "Source point ID is outside the loaded feature vectors")
+    if source_id not in labels:
+        raise ValueError(f"Source point ID {source_id} has no cluster label")
+
+    source_cluster = labels[source_id]
+    if source_cluster < 0:
+        # Skip noise points as source
+        return {"in": None, "out": None}
+
+    # found contains:
+    # "in" or "out": (closest_point_id, hops, distance) or None
+    found: dict[str, tuple[int, int, float] | None] = {"in": None, "out": None}
+
+    # Frontier priority queue contains tuples of (direct_distance, hops, point_id)
+    frontier = [(0.0, 0, source_id)]
+
+    visited: set[int] = set()
+    visited.add(source_id)
+
+    while frontier:
+        source_distance, hops, point_id = heapq.heappop(frontier)
+
+        neighbor_hops = hops + 1
+        for neighbor_id, _stored_distance in graph.get(point_id, ()):
+            if neighbor_id in visited:
+                continue
+            visited.add(neighbor_id)
+
+            neighbor_distance = float(
+                dist_func(points[source_id], points[neighbor_id]))
+            if not math.isfinite(neighbor_distance) or neighbor_distance < 0:
+                raise ValueError(
+                    "Distance function must return a finite non-negative value"
+                )
+            heapq.heappush(
+                frontier,
+                (neighbor_distance, neighbor_hops, neighbor_id),
+            )
+
+            category = "in" if labels[neighbor_id] == source_cluster else "out"
+            current = found[category]
+            if current is None or neighbor_distance < current[2]:
+                found[category] = (
+                    neighbor_id, neighbor_hops, neighbor_distance)
+
+    return found
+
+
+def find_closest_cluster_point_bfs(
+    source_id: int,
+    labels: dict[int, int],
+    graph: Graph,
+    points: np.ndarray,
+    dist_func: DistanceFunction,
+) -> dict[str, tuple[int, int, float] | None]:
+    """Find closest cluster points by BFS depth, ranking same-depth ties by distance.
+
+    The first reachable graph layer containing a point in each category is
+    considered. Within that layer, direct feature distance selects the point.
+    """
+    if source_id < 0 or source_id >= len(points):
+        raise ValueError(
+            "Source point ID is outside the loaded feature vectors")
+    if source_id not in labels:
+        raise ValueError(f"Source point ID {source_id} has no cluster label")
+
+    source_cluster = labels[source_id]
+    if source_cluster < 0:
+        return {"in": None, "out": None}
+
+    found: dict[str, tuple[int, int, float] | None] = {"in": None, "out": None}
+    visited = {source_id}
+    frontier = deque([source_id])
+    hops = 0
+
+    while frontier and any(result is None for result in found.values()):
+        next_frontier: list[int] = []
+        layer_matches: dict[str, tuple[int, float]
+                            | None] = {"in": None, "out": None}
+
+        while frontier:
+            point_id = frontier.popleft()
+            if point_id != source_id and labels.get(point_id, -1) >= 0:
+                category = "in" if labels[point_id] == source_cluster else "out"
+                if found[category] is None:
+                    distance = float(
+                        dist_func(points[source_id], points[point_id]))
+                    if not math.isfinite(distance) or distance < 0:
+                        raise ValueError(
+                            "Distance function must return a finite non-negative value"
+                        )
+                    current = layer_matches[category]
+                    if current is None or (distance, point_id) < current:
+                        layer_matches[category] = (distance, point_id)
+
+            for neighbor_id, _stored_distance in graph.get(point_id, ()):
+                if neighbor_id in visited:
+                    continue
+                if neighbor_id < 0 or neighbor_id >= len(points):
+                    raise ValueError(
+                        "Neighbor point ID is out of the expected range")
+                visited.add(neighbor_id)
+                next_frontier.append(neighbor_id)
+
+        for category, match in layer_matches.items():
+            if match is not None:
+                distance, point_id = match
+                found[category] = (point_id, hops, distance)
+
+        frontier = deque(next_frontier)
+        hops += 1
+
+    return found
+
+
+def compute_score(distance: float, epsilon: float) -> float:
+    """Return the base-10 log score for a non-negative distance."""
+    if not math.isfinite(distance) or distance < 0:
+        raise ValueError("distance must be finite and non-negative")
+    if not math.isfinite(epsilon) or epsilon <= 0:
+        raise ValueError("epsilon must be a finite positive number")
+    return math.log10(distance + epsilon)
+
+
+def find_closest_reachable_distances(
+    points: np.ndarray,
+    graph: Graph,
+    dist_func: DistanceFunction,
+    labels: dict[int, int],
+    algorithm: str = "bfs",
+) -> list[ResultRow]:
+    """Find closest reachable in- and out-of-cluster distances per point."""
+    results: list[ResultRow] = []
+    for source_id in sorted(labels):
+        cluster_id = labels[source_id]
+        if cluster_id < 0:
+            continue
+        if algorithm == "bfs":
+            closest = find_closest_cluster_point_bfs(
+                source_id, labels, graph, points, dist_func
+            )
+        elif algorithm == "dijkstra":
+            closest = find_closest_cluster_point_dijkstra(
+                source_id, labels, graph, points, dist_func
+            )
+        else:
+            raise ValueError(f"Unsupported algorithm: {algorithm}")
+        row: ResultRow = {
+            "point_id": source_id,
+            "cluster_id": cluster_id,
+        }
+        for category in ("in", "out"):
+            match = closest[category]
+            row[f"closest_{category}_point_id"] = match[0] if match else None
+            row[f"closest_{category}_hops"] = match[1] if match else None
+            row[f"closest_{category}_distance"] = match[2] if match else None
+        results.append(row)
+    return results
+
+
+# def write_csv(path: str | Path, results: list[ResultRow]) -> None:
+#     """Write per-point results as CSV."""
+#     fields = [
+#         "point_id",
+#         "cluster_id",
+#         "closest_in_point_id",
+#         "closest_in_hops",
+#         "closest_in_distance",
+#         "closest_out_point_id",
+#         "closest_out_hops",
+#         "closest_out_distance",
+#     ]
+#     output_path = Path(path)
+#     output_path.parent.mkdir(parents=True, exist_ok=True)
+#     with output_path.open("w", encoding="utf-8", newline="") as output_file:
+#         writer = csv.DictWriter(output_file, fieldnames=fields)
+#         writer.writeheader()
+#         writer.writerows(results)
+
+
+def write_txt(path: str | Path, results: list[ResultRow]) -> None:
+    """Write per-point results as plain text."""
+    output_path = Path(path)
+    output_path.parent.mkdir(parents=True, exist_ok=True)
+    # Dump keys in the first line
+    keys = list(results[0].keys()) if results else []
+    with output_path.open("w", encoding="utf-8") as output_file:
+        output_file.write("\t".join(keys) + "\n")
+        for row in results:
+            output_file.write("\t".join(str(row[key]) for key in keys) + "\n")
+
+
+def compute_score_hist(results: list[ResultRow], epsilon: float, num_bins: int) -> dict[str, list[float]]:
+    """Compute histograms of log scores for in-cluster and out-of-cluster distances."""
+    hist = {"in": [0.0] * num_bins, "out": [0.0] * num_bins}
+
+    # First, compute the min and max score
+    min_scores = {"in": float("inf"), "out": float("inf")}
+    max_scores = {"in": float("-inf"), "out": float("-inf")}
+    for row in results:
+        if row["closest_in_distance"] is not None:
+            score = compute_score(float(row["closest_in_distance"]), epsilon)
+            min_scores["in"] = min(min_scores["in"], score)
+            max_scores["in"] = max(max_scores["in"], score)
+        if row["closest_out_distance"] is not None:
+            score = compute_score(float(row["closest_out_distance"]), epsilon)
+            min_scores["out"] = min(min_scores["out"], score)
+            max_scores["out"] = max(max_scores["out"], score)
+    min_score = min(min_scores["in"], min_scores["out"])
+    max_score = max(max_scores["in"], max_scores["out"])
+
+    bin_size = max((max_score - min_score) / num_bins, 1e-12)
+
+    for row in results:
+        if row["closest_in_distance"] is not None:
+            score = compute_score(float(row["closest_in_distance"]), epsilon)
+            bin_index = min(int((score - min_score) / bin_size), num_bins - 1)
+            hist["in"][bin_index] += 1
+        if row["closest_out_distance"] is not None:
+            score = compute_score(float(row["closest_out_distance"]), epsilon)
+            bin_index = min(int((score - min_score) / bin_size), num_bins - 1)
+            hist["out"][bin_index] += 1
+    return hist
+
+
+def print_summary(results: list[ResultRow]) -> None:
+    """Print counts and means of distances by category, and the global score."""
+    print(f"Evaluated {len(results)} clustered points")
+    for category, description in (("in", "in-cluster"), ("out", "out-of-cluster")):
+        values = [
+            float(row[f"closest_{category}_distance"])
+            for row in results
+            if row[f"closest_{category}_distance"] is not None
+        ]
+        if values:
+            print(
+                f"Closest {description}: {len(values)} points; "
+                f"mean distance = {sum(values) / len(values):.6f}"
+            )
+        else:
+            print(f"Closest {description}: no reachable candidates")
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("-p", "--points", required=True,
+                        help="Feature vectors file or directory")
+    parser.add_argument(
+        "-I", "--ids", action="store_true",
+        help="Feature-vector rows include point IDs",
+    )
+    parser.add_argument("-g", "--graph", required=True,
+                        help="kNNG file or directory of files")
+    parser.add_argument(
+        "--knng-format",
+        choices=("M", "MI", "B", "BI"),
+        help="Distance-bearing kNNG format (default: MI)",
+        required=True,
+    )
+    parser.add_argument(
+        "-f", "--dist_func", choices=("l2", "l1", "cosine"),
+        help="Distance function to use",
+        required=True,
+    )
+    parser.add_argument("-c", "--cluster", required=True,
+                        help="Point ID and cluster ID file")
+    parser.add_argument(
+        "-e", "--epsilon", type=float, default=0.001,
+        help="Positive value added before taking log10 (default: 0.001)",
+    )
+    parser.add_argument(
+        "-a", "--algorithm", action="store",
+        help="Traversal algorithm used (BFS or Dijkstra)", choices=("bfs", "dijkstra"),
+        default="bfs",
+    )
+    parser.add_argument(
+        "-v", "--verbose", action="store_true",
+        help="Enable verbose output",
+    )
+    parser.add_argument(
+        "-o", "--output", help="Optional per-point CSV output path")
+    return parser.parse_args()
+
+
+def main() -> None:
+    args = parse_args()
+    try:
+        print("Reading points...")
+        points = read_points(args.points, pattern=None, has_ids=args.ids)
+        print(f"Loaded {len(points)} points")
+
+        print("Reading graph...")
+        graph = read_knng(args.graph, args.knng_format)
+        print(f"Loaded {len(graph)} graph points")
+
+        print("Reading cluster labels...")
+        labels_array = read_cluster_labels(args.cluster)
+        labels = {
+            point_id: int(label)
+            for point_id, label in enumerate(labels_array)
+            if label >= 0
+        }
+        print(f"Loaded {len(labels)} assigned cluster labels")
+
+        distance_functions: dict[str, DistanceFunction] = {
+            "l2": l2_distance,
+            "l1": l1_distance,
+            "cosine": cosine_distance,
+        }
+        results = find_closest_reachable_distances(
+            points, graph, distance_functions[args.dist_func], labels, algorithm=args.algorithm
+        )
+
+        print_summary(results)
+
+        hist = compute_score_hist(results, epsilon=args.epsilon, num_bins=100)
+        print("Histogram of closest reachable distances:")
+        for category, counts in hist.items():
+            print(f"{category}: {counts}")
+
+        if args.output:
+            file_name = args.output + "_scores.txt"
+            print(f"Writing per-point results to {file_name}...")
+            write_txt(file_name, results)
+            print(f"Per-point results written to {file_name}")
+    except (OSError, ValueError) as exc:
+        raise SystemExit(f"Error: {exc}") from exc
+
+
+if __name__ == "__main__":
+    main()
diff --git a/script/benchmark/evaluate_correlation_knng_and_clusters.py b/script/benchmark/evaluate_correlation_knng_and_clusters.py
index 7ddc7e4..392c13e 100644
--- a/script/benchmark/evaluate_correlation_knng_and_clusters.py
+++ b/script/benchmark/evaluate_correlation_knng_and_clusters.py
@@ -35,27 +35,8 @@
 import glob
 from pathlib import Path
 from typing import Iterable, Iterator, Sequence
-
-
-def resolve_files(paths: Sequence[str]) -> list[Path]:
-    """Resolve file, directory, and glob inputs to a unique sorted file list."""
-    resolved: set[Path] = set()
-    for value in paths:
-        path = Path(value)
-        if path.is_dir():
-            matches = set(path.glob("*.txt")) | set(path.glob("knng-*.txt"))
-        elif path.is_file():
-            matches = {path}
-        else:
-            matches = {Path(match) for match in glob.glob(value)}
-
-        resolved.update(match.resolve() for match in matches if match.is_file())
-
-    if not resolved:
-        raise FileNotFoundError(f"No files matched: {paths}")
-
-    return sorted(resolved)
-
+from script.utilities import find_files_in_dir
+from script.data_readers import read_cluster_labels
 
 def iter_knng_records(path: Path, contains_distance: bool) -> Iterator[tuple[int, list[int]]]:
     """Yield (source_id, neighbor_ids) pairs from one KNNG dump file."""
@@ -79,62 +60,16 @@ def iter_knng_records(path: Path, contains_distance: bool) -> Iterator[tuple[int
 
             yield src_id, neighbor_ids
 
-def read_cluster_labels(cluster_path: str | Path) -> dict[int, int]:
-    """Read cluster labels from a file or directory into a point->cluster map."""
-    labels: dict[int, int] = {}
-    files = resolve_files([str(cluster_path)])
-
-    if not files:
-        raise FileNotFoundError(f"No cluster files found in {cluster_path}")
-
-    # Heuristic: determine whether the file format contains point IDs by scanning
-    # the first non-comment, non-empty line.
-    first_file = files[0]
-    contains_ids = False
-    with first_file.open("r") as handle:
-        for line in handle:
-            stripped = line.strip()
-            if not stripped or stripped.startswith("#"):
-                continue
-            if len(stripped.split()) >= 2:
-                contains_ids = True
-            break
-
-    for file in files:
-        with file.open("r") as handle:
-            for line in handle:
-                stripped = line.strip()
-                if not stripped or stripped.startswith("#"):
-                    continue
-
-                items = stripped.split()
-                if contains_ids:
-                    if len(items) < 2:
-                        raise ValueError(f"{file}: invalid cluster row: {line.rstrip()}")
-                    point_id, cluster_id = int(items[0]), int(items[1])
-                else:
-                    point_id = len(labels)
-                    cluster_id = int(items[0])
-
-                if point_id in labels:
-                    raise ValueError(f"Duplicate point ID {point_id} in cluster file {file}")
-                labels[point_id] = cluster_id
-
-    if not labels:
-        raise ValueError(f"No cluster labels were loaded from {cluster_path}")
-
-    return labels
-
 
 def evaluate_knng_clusters(graph_path: str | Path, cluster_path: str | Path, k: int = 10, no_distance: bool = False):
     """Return per-point and aggregate same-cluster neighbor statistics."""
     if k <= 0:
         raise ValueError("-k/--neighbors must be positive")
 
-    labels = read_cluster_labels(cluster_path)
+    labels = read_cluster_labels(cluster_path, "*")
     print(f"Loaded {len(labels)} cluster labels from {cluster_path}")
 
-    files = resolve_files([str(graph_path)])
+    files = find_files_in_dir([str(graph_path)])
 
     per_point: list[tuple[int, int, int, int, float]] = []
     total_same = 0
@@ -179,36 +114,6 @@ def evaluate_knng_clusters(graph_path: str | Path, cluster_path: str | Path, k:
         "num_clustered_points": len(labels),
     }
 
-
-def _parse_k_values(raw_values: Sequence[Sequence[str]] | None) -> list[int]:
-    """Normalize one or more k values from CLI input (supports repeated or comma-separated lists)."""
-    values: list[int] = []
-    raw_entries = raw_values or []
-
-    for group in raw_entries:
-        if isinstance(group, str):
-            group_tokens = [group]
-        else:
-            group_tokens = list(group)
-
-        for token in group_tokens:
-            for item in str(token).split(","):
-                item = item.strip()
-                if not item:
-                    continue
-                try:
-                    value = int(item)
-                except ValueError as exc:
-                    raise argparse.ArgumentTypeError(f"invalid neighbor count: {item!r}") from exc
-                if value <= 0:
-                    raise argparse.ArgumentTypeError(f"neighbor count must be positive: {value}")
-                values.append(value)
-
-    if not values:
-        return [10]
-    return values
-
-
 def _parse_args() -> argparse.Namespace:
     parser = argparse.ArgumentParser(
         description=(
@@ -253,7 +158,7 @@ def _parse_args() -> argparse.Namespace:
         help="Optional output file path for the summary table",
     )
     args = parser.parse_args()
-    args.neighbors = _parse_k_values(args.neighbors or [[10]])
+    args.neighbors = parse_range(args.neighbors or [[10]])
     return args
 
 
diff --git a/script/benchmark/extract_clams_clustering_qualities.sh b/script/benchmark/extract_clams_clustering_qualities.sh
index 74ff9d9..661ddbb 100755
--- a/script/benchmark/extract_clams_clustering_qualities.sh
+++ b/script/benchmark/extract_clams_clustering_qualities.sh
@@ -34,7 +34,7 @@ for file in "$@"; do
     /kNNG k:/ || /^[[:space:]]*k:[[:space:]]*[0-9]+[[:space:]]*$/ { k = lastfield($0); next }
     /Min cluster size/       { mcs = lastfield($0); next }
     /MST-based cluster guess/ { mst_clustering = lastfield($0); next }
-    /#of final clusters/     { fc = lastfield($0); next }
+    /Number of flat clusters:/     { fc = lastfield($0); next }
     /Cluster coverage \(%\)/ { coverage = lastfield($0); next }
     /ARI/                    { ari = lastfield($0); have_ari = 1; next }
     /AMI/ {
diff --git a/script/benchmark/run_hdbscan_bench.py b/script/benchmark/run_hdbscan_bench.py
index 94b8a6c..3670d2b 100644
--- a/script/benchmark/run_hdbscan_bench.py
+++ b/script/benchmark/run_hdbscan_bench.py
@@ -18,7 +18,7 @@
 
 from clustering_utilities import *
 from bench_utilities import *
-
+from script.utilities import *
 
 def parse_options():
     parser = argparse.ArgumentParser(
@@ -101,7 +101,7 @@ def run_hdbscan(points, min_cluster_size, min_samples,
     # Evaluate the clustering quality
     if gt_file:
         print('\nLoading ground truth data')
-        gt_labels = read_label_data(gt_file)
+        gt_labels = read_cluster_labels(gt_file, "*")
 
         print('\nEvaluating clustering quality')
         eval_clusters(clusters.labels_, gt_labels)
@@ -152,7 +152,7 @@ def main():
     opts = parse_options()
 
 
-    points = read_point_data(opts.point_data_path, opts.has_ids)
+    points = read_points(opts.point_data_path, "*", opts.has_ids)
 
     min_cluster_size_list = parse_range(opts.min_cluster_size_range)
     min_samples_list = parse_range(opts.min_samples_range)
diff --git a/script/data_readers.py b/script/data_readers.py
new file mode 100644
index 0000000..3a9393a
--- /dev/null
+++ b/script/data_readers.py
@@ -0,0 +1,255 @@
+# Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS
+# Project Developers. See the top-level COPYRIGHT file for details.
+
+from __future__ import annotations
+
+import math
+from pathlib import Path
+
+import numpy as np
+from .utilities import find_files_in_dir
+
+def read_points(
+    data_path: str | Path,
+    pattern: str | None = None,
+    has_ids: bool = True,
+) -> np.ndarray:
+    """Load whitespace-separated feature vectors into a dense NumPy array.
+
+    Each non-empty row contains feature values, optionally preceded by a point
+    ID. When IDs are present, missing IDs are represented by 0-based row numbers.
+    Multiple files require explicit IDs so records can be combined unambiguously.
+    """
+    print(f"Loading data from {data_path}", flush=True)
+    files = find_files_in_dir([data_path], pattern)
+    if len(files) > 1 and not has_ids:
+        raise ValueError("Multiple point files require point IDs")
+
+    points_table: dict[int, list[float]] = {}
+    dimensions: int | None = None
+    for file in files:
+        with file.open("r", encoding="utf-8") as input_file:
+            for line_number, line in enumerate(input_file, start=1):
+                items = line.split()
+                if not items:
+                    continue
+
+                if has_ids:
+                    try:
+                        point_id = int(items.pop(0))
+                    except ValueError as exc:
+                        raise ValueError(
+                            f"{file}:{line_number}: invalid point ID"
+                        ) from exc
+                    if point_id < 0:
+                        raise ValueError(f"{file}:{line_number}: point IDs must be non-negative")
+                else:
+                    point_id = len(points_table)
+
+                if not items:
+                    raise ValueError(f"{file}:{line_number}: feature vector is empty")
+                if point_id in points_table:
+                    raise ValueError(f"{file}:{line_number}: duplicate point ID {point_id}")
+
+                try:
+                    feature = [float(value) for value in items]
+                except ValueError as exc:
+                    raise ValueError(f"{file}:{line_number}: invalid feature value") from exc
+
+                if dimensions is None:
+                    dimensions = len(feature)
+                elif len(feature) != dimensions:
+                    raise ValueError(
+                        f"{file}:{line_number}: dimension mismatch: "
+                        f"{len(feature)} != {dimensions}"
+                    )
+                points_table[point_id] = feature
+
+    if not points_table or dimensions is None:
+        raise ValueError(f"No feature vectors were loaded from {data_path}")
+
+    print(f"Loaded {len(points_table)} items from {len(files)} files", flush=True)
+
+    points = np.zeros((max(points_table) + 1, dimensions), dtype=float)
+    for point_id, feature in points_table.items():
+        points[point_id] = feature
+
+    return points
+
+
+def read_cluster_labels(
+    data_path: str | Path,
+    pattern: str | None = None,
+) -> np.ndarray:
+    """Load cluster labels from one or more files into a dense integer array.
+
+    Rows may contain only a cluster ID (point ID is the 0-based row number) or a
+    point ID followed by a cluster ID. Blank lines and comments are ignored.
+    Multiple files require explicit IDs.
+    """
+    files = find_files_in_dir([data_path], pattern)
+
+    print(f"Found {len(files)} files in {data_path}", flush=True)
+    first_data_row: list[str] | None = None
+    for file in files:
+        with file.open("r", encoding="utf-8") as input_file:
+            for line in input_file:
+                stripped = line.strip()
+                if stripped and not stripped.startswith("#"):
+                    first_data_row = stripped.split()
+                    break
+        if first_data_row is not None:
+            break
+
+    if first_data_row is None:
+        raise ValueError(f"No cluster labels were found in {data_path}")
+
+    contains_ids = len(first_data_row) >= 2
+
+    if contains_ids:
+        print("Loading point IDs and labels", flush=True)
+    else:
+        print("Loading only labels", flush=True)
+
+    if len(files) > 1 and not contains_ids:
+        raise ValueError("Multiple cluster files require point IDs")
+
+    labels_dict: dict[int, int] = {}
+    for file in files:
+        with file.open("r", encoding="utf-8") as input_file:
+            for line_number, line in enumerate(input_file, start=1):
+                stripped = line.strip()
+                if not stripped or stripped.startswith("#"):
+                    continue
+                items = stripped.split()
+
+                try:
+                    if contains_ids:
+                        if len(items) < 2:
+                            raise ValueError("expected point ID and cluster ID")
+                        point_id, label = int(items[0]), int(items[1])
+                    else:
+                        point_id, label = len(labels_dict), int(items[0])
+                except ValueError as exc:
+                    raise ValueError(f"{file}:{line_number}: invalid cluster row") from exc
+
+                if point_id < 0:
+                    raise ValueError(f"{file}:{line_number}: point IDs must be non-negative")
+                if point_id in labels_dict:
+                    raise ValueError(f"{file}:{line_number}: duplicate point ID {point_id}")
+                labels_dict[point_id] = label
+
+    if not labels_dict:
+        raise ValueError(f"No cluster labels were found in {data_path}")
+
+    print(f"Loaded {len(labels_dict)} items from {len(files)} files", flush=True)
+    max_id = max(labels_dict)
+    print(f"Max ID: {max_id}", flush=True)
+
+    labels = np.full(max_id + 1, -1)
+    for point_id, label in labels_dict.items():
+        labels[point_id] = label
+
+    return labels
+
+Graph: TypeAlias = dict[int, list[tuple[int, float]]]
+
+def read_knng(
+    graph_path: str | Path,
+    knng_format: str,
+    path_pattern: str | None = None,
+) -> Graph:
+    """Load a kNN graph as source IDs mapped to neighbor/distance pairs.
+
+    ``knng_format`` selects ``N`` (neighbor IDs only), ``B`` (blocked ID rows
+    followed by their distance rows), or ``M`` (alternating ID and distance
+    rows). Append ``I`` to any format when each ID row starts with its source
+    point ID; otherwise source IDs are assigned sequentially across the files.
+    Distances are ``None`` for the ID-only format. With explicit source IDs,
+    the corresponding first distance is treated as the dummy source distance.
+    """
+    format_code = knng_format.upper()
+    base_formats = [code for code in format_code if code in "NBM"]
+    if len(base_formats) != 1 or any(code not in "NBMI" for code in format_code):
+        raise ValueError(f"Unsupported kNNG format: {knng_format!r}")
+    base_format = base_formats[0]
+    contains_ids = "I" in format_code
+
+    files = find_files_in_dir([graph_path], path_pattern)
+    graph: dict[int, list[tuple[int, float | None]]] = {}
+    next_source_id = 0
+
+    def parse_id_row(file: Path, line_number: int, line: str) -> tuple[int, list[int]]:
+        nonlocal next_source_id
+        try:
+            values = [int(value) for value in line.split()]
+        except ValueError as exc:
+            raise ValueError(f"{file}:{line_number}: invalid neighbor ID row") from exc
+        if contains_ids:
+            if not values:
+                raise ValueError(f"{file}:{line_number}: missing source point ID")
+            source_id, neighbor_ids = values[0], values[1:]
+        else:
+            source_id, neighbor_ids = next_source_id, values
+            next_source_id += 1
+        if source_id < 0 or any(neighbor_id < 0 for neighbor_id in neighbor_ids):
+            raise ValueError(f"{file}:{line_number}: point IDs must be non-negative")
+        if source_id in graph:
+            raise ValueError(f"{file}:{line_number}: duplicate source point ID {source_id}")
+        return source_id, neighbor_ids
+
+    def parse_distances(file: Path, line_number: int, line: str) -> list[float]:
+        try:
+            distances = [float(value) for value in line.split()]
+        except ValueError as exc:
+            raise ValueError(f"{file}:{line_number}: invalid distance row") from exc
+        if any(not math.isfinite(distance) or distance < 0 for distance in distances):
+            raise ValueError(f"{file}:{line_number}: distances must be finite and non-negative")
+        return distances
+
+    for file in files:
+        with file.open("r", encoding="utf-8") as input_file:
+            lines = [
+                (line_number, line.strip())
+                for line_number, line in enumerate(input_file, start=1)
+                if line.strip() and not line.lstrip().startswith("#")
+            ]
+
+        if base_format == "N":
+            for line_number, line in lines:
+                source_id, neighbor_ids = parse_id_row(file, line_number, line)
+                graph[source_id] = [(neighbor_id, None) for neighbor_id in neighbor_ids]
+            continue
+
+        if len(lines) % 2 != 0:
+            raise ValueError(f"{file}: expected an even number of kNNG lines")
+        half = len(lines) // 2 if base_format == "B" else None
+        if base_format == "B":
+            id_rows = lines[:half]
+            distance_rows = lines[half:]
+        else:
+            id_rows = lines[::2]
+            distance_rows = lines[1::2]
+
+        for (id_line_number, id_line), (distance_line_number, distance_line) in zip(
+            id_rows, distance_rows
+        ):
+            source_id, neighbor_ids = parse_id_row(file, id_line_number, id_line)
+            distances = parse_distances(file, distance_line_number, distance_line)
+            if contains_ids:
+                if len(distances) != len(neighbor_ids) + 1:
+                    raise ValueError(
+                        f"{file}:{distance_line_number}: expected a dummy distance and "
+                        f"{len(neighbor_ids)} neighbor distances"
+                    )
+                distances = distances[1:]
+            elif len(distances) != len(neighbor_ids):
+                raise ValueError(
+                    f"{file}:{distance_line_number}: {len(neighbor_ids)} neighbors but "
+                    f"{len(distances)} distances"
+                )
+            graph[source_id] = list(zip(neighbor_ids, distances))
+
+    if not graph:
+        raise ValueError(f"No kNNG rows were loaded from {graph_path}")
+    return graph
diff --git a/script/etc/reassign_sequential_cluster_ids.py b/script/etc/reassign_sequential_cluster_ids.py
index 12ed8f3..d0aaac8 100644
--- a/script/etc/reassign_sequential_cluster_ids.py
+++ b/script/etc/reassign_sequential_cluster_ids.py
@@ -28,7 +28,8 @@
 import argparse
 import numpy as np
 from script.benchmark.clustering_utilities import *
-
+from script.utilities import *
+from script.data_readers import *
 
 def main():
     parser = argparse.ArgumentParser(
@@ -43,7 +44,7 @@ def main():
                         help='Path to the output file')
     args = parser.parse_args()
 
-    cluster_ids = read_label_data(args.input_cluster_ids_path)
+    cluster_ids = read_cluster_labels(args.input_cluster_ids_path, "*")
 
     # Make new cluster IDs
     new_ids = dict()
diff --git a/script/utilities.py b/script/utilities.py
new file mode 100644
index 0000000..457b297
--- /dev/null
+++ b/script/utilities.py
@@ -0,0 +1,178 @@
+# Copyright 2023-2026 Lawrence Livermore National Security, LLC and other ClaMS
+# Project Developers. See the top-level COPYRIGHT file for details.
+
+from __future__ import annotations
+import datetime
+from collections.abc import Sequence
+import glob
+from pathlib import Path
+import subprocess
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+
+
+def _as_vector_pair(
+    vec1: ArrayLike, vec2: ArrayLike
+) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
+    """Convert and validate two finite, non-empty one-dimensional vectors."""
+    first = np.asarray(vec1, dtype=np.float64)
+    second = np.asarray(vec2, dtype=np.float64)
+    if first.ndim != 1 or second.ndim != 1:
+        raise ValueError("Vectors must be one-dimensional")
+    if first.size != second.size:
+        raise ValueError("Vectors must have the same length")
+    if first.size == 0:
+        raise ValueError("Vectors must not be empty")
+    if not np.isfinite(first).all() or not np.isfinite(second).all():
+        raise ValueError("Vectors must contain only finite values")
+    return first, second
+
+
+def show_time_now():
+    now = datetime.datetime.now()
+    print(now.strftime("%Y-%m-%d %H:%M:%S"))
+
+
+def find_files_in_dir(
+    paths: Sequence[str | Path], pattern: str | None = None
+) -> list[Path]:
+    """Resolve files, directories, and glob patterns to unique sorted paths.
+
+    ``pattern`` filters entries inside directory inputs. Nonexistent path inputs
+    are treated as glob patterns. Directory searches are not recursive.
+    """
+    resolved: set[Path] = set()
+    for value in paths:
+        path = Path(value)
+        if path.is_dir():
+            matches = path.glob(pattern or "*")
+        elif path.is_file():
+            matches = (path,)
+        else:
+            # Nonexistent inputs may still be glob expressions, such as *.txt.
+            matches = (Path(match) for match in glob.glob(str(value)))
+
+        resolved.update(match.resolve() for match in matches if match.is_file())
+
+    if not resolved:
+        formatted_paths = ", ".join(map(str, paths))
+        raise FileNotFoundError(f"No files matched: {formatted_paths}")
+
+    return sorted(resolved)
+
+
+def create_dir(dir_path: str | Path) -> None:
+    """Create a directory and any missing parent directories if needed."""
+    Path(dir_path).mkdir(parents=True, exist_ok=True)
+
+
+def grep_file(file: str | Path, text_pattern: str) -> list[str]:
+    """Return the lines in ``file`` that contain the given substring."""
+    with Path(file).open("r", encoding="utf-8") as input_file:
+        return [line for line in input_file if text_pattern in line]
+
+
+def parse_range(range_str: str | None) -> list[int | None]:
+    """Parse an inclusive integer range or comma-separated integer list.
+
+    Examples: ``"1-3"`` becomes ``[1, 2, 3]`` and ``"1,3"`` becomes
+    ``[1, 3]``. ``None`` is preserved as ``[None]`` for optional CLI values.
+    """
+    if range_str is None:
+        return [None]
+
+    value = range_str.strip()
+    if not value:
+        raise ValueError("range must not be empty")
+
+    if "-" in value:
+        parts = value.split("-")
+        if len(parts) != 2:
+            raise ValueError(f"invalid integer range: {range_str!r}")
+        start, end = (int(part.strip()) for part in parts)
+        if start > end:
+            raise ValueError(f"range start must not exceed its end: {range_str!r}")
+        return list(range(start, end + 1))
+
+    return [int(part.strip()) for part in value.split(",")]
+
+
+def parse_float_list(float_list_str: str | None) -> list[float | None]:
+    """Parse comma-separated floating-point values into a list.
+
+    ``None`` is preserved as ``[None]`` for optional CLI values.
+    """
+    if float_list_str is None:
+        return [None]
+
+    value = float_list_str.strip()
+    if not value:
+        raise ValueError("float list must not be empty")
+    return [float(part.strip()) for part in value.split(",")]
+
+
+def execute_cmd(
+    log_file: str | Path,
+    err_log_file: str | Path,
+    command: str,
+    cwd: str | Path = "./",
+) -> None:
+    """Run a shell command, append its output to logs, and exit on failure.
+
+    The command is interpreted by the system shell. Standard output and
+    standard error are appended to their respective log files.
+    """
+    print(f"In {Path(cwd).resolve()}")
+    print(f"Command: {command}")
+    result = subprocess.run(
+        command,
+        shell=True,
+        cwd=cwd,
+        stdout=subprocess.PIPE,
+        stderr=subprocess.PIPE,
+        text=True,
+        encoding="utf-8",
+        errors="replace",
+        check=False,
+    )
+
+    with Path(log_file).open("a", encoding="utf-8") as output_file:
+        output_file.write(f"Command: {command}\n")
+        output_file.write(result.stdout)
+    with Path(err_log_file).open("a", encoding="utf-8") as error_file:
+        error_file.write(result.stderr)
+
+    if result.returncode == 0:
+        print("Command executed successfully.\n")
+    else:
+        print(f"Error executing command. See {err_log_file}.\n")
+        raise SystemExit(result.returncode)
+
+
+# Distance (similarity) functions
+# The smaller the distance value, the more similar the vectors are.
+# Inputs are numpy arrays or lists of floats.
+def l2_distance(vec1: ArrayLike, vec2: ArrayLike) -> float:
+    """Compute the Euclidean (L2) distance between two vectors."""
+    first, second = _as_vector_pair(vec1, vec2)
+    return float(np.linalg.norm(first - second))
+
+
+def l1_distance(vec1: ArrayLike, vec2: ArrayLike) -> float:
+    """Compute the Manhattan (L1) distance between two vectors."""
+    first, second = _as_vector_pair(vec1, vec2)
+    return float(np.sum(np.abs(first - second)))
+
+
+def cosine_distance(vec1: ArrayLike, vec2: ArrayLike) -> float:
+    """Compute the Cosine distance between two vectors (1 - cosine similarity)."""
+    first, second = _as_vector_pair(vec1, vec2)
+    norm1 = np.linalg.norm(first)
+    norm2 = np.linalg.norm(second)
+    if norm1 == 0 or norm2 == 0:
+        raise ValueError("Vectors must not be zero vectors")
+
+    cosine_similarity = np.dot(first / norm1, second / norm2)
+    cosine_similarity = np.clip(cosine_similarity, -1.0, 1.0)
+    return float(1.0 - cosine_similarity)