From 5e93be28feecfbce117b48160468fe045ffbdf09 Mon Sep 17 00:00:00 2001
From: vc1492a <vc1492a@gmail.com>
Date: Fri, 20 Mar 2026 11:12:01 -0700
Subject: [PATCH 1/5] feat: add parallel distance computation and vectorized
 pipeline

Rewrite the distance computation engine from scratch on top of v0.3.5:

- Vectorized kNN distances using NumPy broadcasting with chunked
  processing for memory efficiency and progress bar support
- Add n_jobs parameter for cross-cluster multiprocessing via
  concurrent.futures (n_jobs=-1 uses all cores)
- Restructure Numba path with non-generator kernels that support
  numba.prange for thread-level parallelism
- Optional scipy.spatial.distance.cdist and scipy.special.erf
  acceleration when scipy is available
- Vectorize _standard_distances, _prob_distances, and
  _norm_prob_outlier_factor pipeline methods
- Fully backward-compatible: all existing API calls work unchanged

Closes #36

Made-with: Cursor
---
 PyNomaly/loop.py   | 291 +++++++++++++++++++++++++++++++++++++++------
 README.md          |  42 ++++++-
 tests/test_loop.py |  50 ++++++++
 3 files changed, 343 insertions(+), 40 deletions(-)

diff --git a/PyNomaly/loop.py b/PyNomaly/loop.py
index 48f499b..46bdf31 100644
--- a/PyNomaly/loop.py
+++ b/PyNomaly/loop.py
@@ -1,12 +1,79 @@
+from concurrent.futures import ProcessPoolExecutor, as_completed
 from math import erf, sqrt
 import numpy as np
+import os
 from python_utils.terminal import get_terminal_size
 import sys
 from typing import Tuple, Union
 import warnings
 
+try:
+    from scipy.spatial.distance import cdist as _scipy_cdist
+except ImportError:
+    _scipy_cdist = None
+
+try:
+    from scipy.special import erf as _scipy_erf
+except ImportError:
+    _scipy_erf = None
+
 try:
     import numba
+
+    def _numba_distance_kernel_parallel(clust_points_vector, n_neighbors):
+        n = clust_points_vector.shape[0]
+        d_features = clust_points_vector.shape[1]
+        local_distances = np.full((n, n_neighbors), 9e10, dtype=np.float64)
+        local_indexes = np.full((n, n_neighbors), 0, dtype=np.int64)
+        for i in numba.prange(n):
+            for j in range(n):
+                if i == j:
+                    continue
+                d = 0.0
+                for f in range(d_features):
+                    diff_val = clust_points_vector[i, f] - clust_points_vector[j, f]
+                    d += diff_val * diff_val
+                d = d ** 0.5
+                idx_max = 0
+                for k in range(1, n_neighbors):
+                    if local_distances[i, k] > local_distances[i, idx_max]:
+                        idx_max = k
+                if d < local_distances[i, idx_max]:
+                    local_distances[i, idx_max] = d
+                    local_indexes[i, idx_max] = j
+        return local_distances, local_indexes
+
+    _numba_kernel_parallel = numba.jit(
+        _numba_distance_kernel_parallel, nopython=True, parallel=True, cache=True
+    )
+
+    def _numba_distance_kernel_sequential(clust_points_vector, n_neighbors):
+        n = clust_points_vector.shape[0]
+        d_features = clust_points_vector.shape[1]
+        local_distances = np.full((n, n_neighbors), 9e10, dtype=np.float64)
+        local_indexes = np.full((n, n_neighbors), 0, dtype=np.int64)
+        for i in range(n):
+            for j in range(n):
+                if i == j:
+                    continue
+                d = 0.0
+                for f in range(d_features):
+                    diff_val = clust_points_vector[i, f] - clust_points_vector[j, f]
+                    d += diff_val * diff_val
+                d = d ** 0.5
+                idx_max = 0
+                for k in range(1, n_neighbors):
+                    if local_distances[i, k] > local_distances[i, idx_max]:
+                        idx_max = k
+                if d < local_distances[i, idx_max]:
+                    local_distances[i, idx_max] = d
+                    local_indexes[i, idx_max] = j
+        return local_distances, local_indexes
+
+    _numba_kernel_sequential = numba.jit(
+        _numba_distance_kernel_sequential, nopython=True, cache=True
+    )
+
 except ImportError:
     pass
 
@@ -15,6 +82,40 @@
 __license__ = "Apache License, Version 2.0"
 
 
+def _cluster_distances_worker(args):
+    """Top-level worker for ProcessPoolExecutor (must be picklable)."""
+    clust_points_vector, global_indices, n_neighbors = args
+    n = clust_points_vector.shape[0]
+
+    if clust_points_vector.ndim == 1:
+        clust_points_vector = clust_points_vector.reshape(-1, 1)
+
+    local_distances = np.full((n, n_neighbors), 9e10, dtype=float)
+    local_indexes = np.full((n, n_neighbors), 9e10, dtype=float)
+
+    chunk_size = min(256, n)
+    for chunk_start in range(0, n, chunk_size):
+        chunk_end = min(chunk_start + chunk_size, n)
+        chunk = clust_points_vector[chunk_start:chunk_end]
+
+        if _scipy_cdist is not None:
+            dist = _scipy_cdist(chunk, clust_points_vector, metric="euclidean")
+        else:
+            diff = chunk[:, np.newaxis, :] - clust_points_vector[np.newaxis, :, :]
+            dist = np.sqrt((diff ** 2).sum(axis=2))
+
+        row_idx = np.arange(chunk_end - chunk_start)
+        dist[row_idx, row_idx + chunk_start] = np.inf
+
+        knn_idx = np.argpartition(dist, n_neighbors, axis=1)[:, :n_neighbors]
+        knn_dists = np.take_along_axis(dist, knn_idx, axis=1)
+
+        local_distances[chunk_start:chunk_end] = knn_dists
+        local_indexes[chunk_start:chunk_end] = global_indices[knn_idx]
+
+    return local_distances, local_indexes, global_indices
+
+
 # Custom Exceptions
 class PyNomalyError(Exception):
     """Base exception for PyNomaly."""
@@ -74,6 +175,9 @@ class LocalOutlierProbability(object):
     sample (optional, default 10)
     :param cluster_labels: a numpy array of cluster assignments w.r.t. each 
     sample (optional, default None)
+    :param n_jobs: the number of parallel workers for distance computation.
+    Use -1 to use all available CPU cores, or 1 for sequential processing
+    (optional, default 1)
     :return:
     """ """
 
@@ -315,7 +419,8 @@ def new_f(*args, **kwds):
                     "n_neighbors": {"type": types[5]},
                     "cluster_labels": {"type": types[6]},
                     "use_numba": {"type": types[7]},
-                    "progress_bar": {"type": types[8]},
+                    "n_jobs": {"type": types[8]},
+                    "progress_bar": {"type": types[9]},
                 }
                 for x in kwds:
                     opt_types[x]["value"] = kwds[x]
@@ -348,6 +453,7 @@ def new_f(*args, **kwds):
         (int, np.integer),
         list,
         bool,
+        (int, np.integer),
         bool,
     )
     def __init__(
@@ -359,6 +465,7 @@ def __init__(
         n_neighbors=10,
         cluster_labels=None,
         use_numba=False,
+        n_jobs=1,
         progress_bar=False,
     ) -> None:
         self.data = data
@@ -368,6 +475,7 @@ def __init__(
         self.n_neighbors = n_neighbors
         self.cluster_labels = cluster_labels
         self.use_numba = use_numba
+        self.n_jobs = n_jobs
         self.points_vector = None
         self.prob_distances = None
         self.prob_distances_ev = None
@@ -383,6 +491,13 @@ def __init__(
                 "Numba is not available, falling back to pure python mode.", UserWarning
             )
 
+        if self.n_jobs < -1 or self.n_jobs == 0:
+            warnings.warn(
+                "n_jobs must be -1 or a positive integer. Defaulting to 1.",
+                UserWarning,
+            )
+            self.n_jobs = 1
+
         self._validate_inputs()
         self._check_extent()
 
@@ -444,10 +559,8 @@ def _norm_prob_outlier_factor(
         outlier factor of the input observation.
         :return: the normalized probabilistic outlier factor.
         """
-        npofs = []
-        for i in ev_probabilistic_outlier_factor:
-            npofs.append(extent * sqrt(i))
-        return npofs
+        ev_arr = np.array(ev_probabilistic_outlier_factor, dtype=float)
+        return (extent * np.sqrt(ev_arr)).tolist()
 
     @staticmethod
     def _local_outlier_probability(
@@ -461,11 +574,14 @@ def _local_outlier_probability(
         input observation.
         :return: the local outlier probability.
         """
-        erf_vec = np.vectorize(erf)
         if np.all(plof_val == nplof_val):
             return np.zeros(plof_val.shape)
-        else:
-            return np.maximum(0, erf_vec(plof_val / (nplof_val * np.sqrt(2.0))))
+        plof_f = np.asarray(plof_val, dtype=float)
+        nplof_f = np.asarray(nplof_val, dtype=float)
+        if _scipy_erf is not None:
+            return np.maximum(0, _scipy_erf(plof_f / (nplof_f * np.sqrt(2.0))))
+        erf_vec = np.vectorize(erf)
+        return np.maximum(0, erf_vec(plof_f / (nplof_f * np.sqrt(2.0))))
 
     def _n_observations(self) -> int:
         """
@@ -564,6 +680,102 @@ def _compute_distance_and_neighbor_matrix(
 
             yield distances, indexes, i
 
+    def _distances_vectorized(
+        self, clusters, distances, indexes, progress_bar
+    ) -> None:
+        """Vectorized kNN distance computation with chunked progress."""
+        progress = "="
+        total_points = sum(cv.shape[0] for cv, _ in clusters)
+        completed = 0
+
+        for clust_points_vector, global_indices in clusters:
+            n = clust_points_vector.shape[0]
+
+            if clust_points_vector.ndim == 1:
+                clust_points_vector = clust_points_vector.reshape(-1, 1)
+
+            chunk_size = min(256, n)
+            for chunk_start in range(0, n, chunk_size):
+                chunk_end = min(chunk_start + chunk_size, n)
+                chunk = clust_points_vector[chunk_start:chunk_end]
+
+                if _scipy_cdist is not None:
+                    dist = _scipy_cdist(
+                        chunk, clust_points_vector, metric="euclidean"
+                    )
+                else:
+                    diff = (
+                        chunk[:, np.newaxis, :]
+                        - clust_points_vector[np.newaxis, :, :]
+                    )
+                    dist = np.sqrt((diff ** 2).sum(axis=2))
+
+                row_idx = np.arange(chunk_end - chunk_start)
+                dist[row_idx, row_idx + chunk_start] = np.inf
+
+                knn_idx = np.argpartition(dist, self.n_neighbors, axis=1)[
+                    :, : self.n_neighbors
+                ]
+                knn_dists = np.take_along_axis(dist, knn_idx, axis=1)
+
+                chunk_global = global_indices[chunk_start:chunk_end]
+                distances[chunk_global] = knn_dists
+                indexes[chunk_global] = global_indices[knn_idx]
+
+                completed += chunk_end - chunk_start
+                if progress_bar:
+                    progress = Utils.emit_progress_bar(
+                        progress, completed, total_points
+                    )
+
+    def _distances_numba(
+        self, clusters, distances, indexes, progress_bar, parallel=False
+    ) -> None:
+        """Numba-accelerated distance computation."""
+        progress = "="
+        kernel = _numba_kernel_parallel if parallel else _numba_kernel_sequential
+
+        for idx, (clust_points_vector, global_indices) in enumerate(clusters):
+            if clust_points_vector.ndim == 1:
+                clust_points_vector = clust_points_vector.reshape(-1, 1)
+
+            local_dists, local_idxs = kernel(
+                clust_points_vector.astype(np.float64), self.n_neighbors
+            )
+
+            distances[global_indices] = local_dists
+            indexes[global_indices] = global_indices[local_idxs]
+
+            if progress_bar:
+                progress = Utils.emit_progress_bar(
+                    progress, idx + 1, len(clusters)
+                )
+
+    def _distances_parallel(
+        self, clusters, distances, indexes, n_jobs, progress_bar
+    ) -> None:
+        """Parallel distance computation across clusters via multiprocessing."""
+        progress = "="
+        worker_args = [
+            (cv, gi, self.n_neighbors) for cv, gi in clusters
+        ]
+
+        with ProcessPoolExecutor(max_workers=n_jobs) as executor:
+            futures = {
+                executor.submit(_cluster_distances_worker, args): idx
+                for idx, args in enumerate(worker_args)
+            }
+            completed_clusters = 0
+            for future in as_completed(futures):
+                local_dists, local_idxs, global_indices = future.result()
+                distances[global_indices] = local_dists
+                indexes[global_indices] = local_idxs
+                completed_clusters += 1
+                if progress_bar:
+                    progress = Utils.emit_progress_bar(
+                        progress, completed_clusters, len(clusters)
+                    )
+
     def _distances(self, progress_bar: bool = False) -> None:
         """
         Provides the distances between each observation and it's closest
@@ -576,27 +788,40 @@ def _distances(self, progress_bar: bool = False) -> None:
         distances = np.full(
             [self._n_observations(), self.n_neighbors], 9e10, dtype=float
         )
-        indexes = np.full([self._n_observations(), self.n_neighbors], 9e10, dtype=float)
-        self.points_vector = self._convert_to_array(self.data)
-        compute = (
-            numba.jit(self._compute_distance_and_neighbor_matrix, cache=True)
-            if self.use_numba
-            else self._compute_distance_and_neighbor_matrix
+        indexes = np.full(
+            [self._n_observations(), self.n_neighbors], 9e10, dtype=float
         )
-        progress = "="
-        for cluster_id in set(self._cluster_labels()):
-            indices = np.where(self._cluster_labels() == cluster_id)
+        self.points_vector = self._convert_to_array(self.data)
+
+        cluster_labels = self._cluster_labels()
+        cluster_ids = sorted(set(cluster_labels))
+
+        clusters = []
+        for cluster_id in cluster_ids:
+            indices = np.where(cluster_labels == cluster_id)
             clust_points_vector = np.array(
                 self.points_vector.take(indices, axis=0)[0], dtype=np.float64
             )
-            # a generator that yields an updated distance matrix on each loop
-            for c in compute(clust_points_vector, indices, distances, indexes):
-                distances, indexes, i = c
-                # update the progress bar
-                if progress_bar is True:
-                    progress = Utils.emit_progress_bar(
-                        progress, i + 1, clust_points_vector.shape[0]
-                    )
+            clusters.append((clust_points_vector, indices[0]))
+
+        n_jobs = self.n_jobs
+        if n_jobs == -1:
+            n_jobs = os.cpu_count() or 1
+        n_jobs = min(n_jobs, len(clusters))
+
+        if self.use_numba:
+            self._distances_numba(
+                clusters, distances, indexes, progress_bar,
+                parallel=(n_jobs > 1)
+            )
+        elif n_jobs > 1 and len(clusters) > 1:
+            self._distances_parallel(
+                clusters, distances, indexes, n_jobs, progress_bar
+            )
+        else:
+            self._distances_vectorized(
+                clusters, distances, indexes, progress_bar
+            )
 
         self.distance_matrix = distances
         self.neighbor_matrix = indexes
@@ -626,18 +851,14 @@ def _standard_distances(self, data_store: np.ndarray) -> np.ndarray:
         Calculated the standard distance for each observation in the input
         data. First calculates the cardinality and then calculates the standard
         distance with respect to each observation.
-        :param data_store:
         :param data_store: the storage matrix that collects information on
         each observation.
         :return: the updated storage matrix that collects information on
         each observation.
         """
-        cardinality = [self.n_neighbors] * self._n_observations()
-        vals = data_store[:, 3].tolist()
-        std_distances = []
-        for c, v in zip(cardinality, vals):
-            std_distances.append(self._standard_distance(c, v))
-        return np.hstack((data_store, np.array([std_distances]).T))
+        ssd_vals = data_store[:, 3].astype(float)
+        std_distances = np.sqrt(ssd_vals / self.n_neighbors)
+        return np.hstack((data_store, std_distances.reshape(-1, 1)))
 
     def _prob_distances(self, data_store: np.ndarray) -> np.ndarray:
         """
@@ -648,10 +869,8 @@ def _prob_distances(self, data_store: np.ndarray) -> np.ndarray:
         :return: the updated storage matrix that collects information on
         each observation.
         """
-        prob_distances = []
-        for i in range(data_store[:, 4].shape[0]):
-            prob_distances.append(self._prob_distance(self.extent, data_store[:, 4][i]))
-        return np.hstack((data_store, np.array([prob_distances]).T))
+        prob_distances = self.extent * data_store[:, 4].astype(float)
+        return np.hstack((data_store, prob_distances.reshape(-1, 1)))
 
     def _prob_distances_ev(self, data_store) -> np.ndarray:
         """
diff --git a/README.md b/README.md
index d639fd8..4faab28 100644
--- a/README.md
+++ b/README.md
@@ -43,12 +43,14 @@ to calculate the Local Outlier Probability of each sample.
 - numpy >= 1.16.3
 - python-utils >= 2.3.0
 - (optional) numba >= 0.45.1
+- (optional) scipy >= 1.3.0
 
-Numba just-in-time (JIT) compiles the function with calculates the Euclidean 
+Numba just-in-time (JIT) compiles the function which calculates the Euclidean 
 distance between observations, providing a reduction in computation time 
 (significantly when a large number of observations are scored). Numba is not a 
 requirement and PyNomaly may still be used solely with numpy if desired
-(details below). 
+(details below). When scipy is available, PyNomaly uses its optimized distance 
+computation and error function implementations for additional performance gains.
 
 ## Quick Start
 
@@ -101,6 +103,36 @@ normalization scheme prior to using LoOP, especially when working with multiple
 Users must also appropriately handle missing values prior to using LoOP, as LoOP does not support Pandas
 DataFrames or Numpy arrays with missing values.
 
+### Parallel Processing
+
+PyNomaly supports parallel distance computation via the `n_jobs` parameter.
+When data includes multiple clusters (via `cluster_labels`), each cluster's 
+distance computation is independent and can be processed across multiple CPU cores:
+
+```python
+from PyNomaly import loop
+from sklearn.cluster import DBSCAN
+db = DBSCAN(eps=0.6, min_samples=50).fit(data)
+m = loop.LocalOutlierProbability(
+    data, extent=2, n_neighbors=20,
+    cluster_labels=list(db.labels_), n_jobs=-1
+).fit()
+scores = m.local_outlier_probabilities
+print(scores)
+```
+
+Set `n_jobs=-1` to use all available CPU cores, or specify a positive integer 
+for a fixed number of workers. The default value of `n_jobs=1` runs sequentially. 
+Parallelism is most beneficial when the data contains multiple clusters. For 
+single-cluster data, `n_jobs` has no effect as there is only one unit of work.
+
+When using Numba (`use_numba=True`) with `n_jobs > 1`, PyNomaly uses Numba's 
+thread-level parallelism (`prange`) instead of multiprocessing.
+
+Note that parallel processing is only applicable when raw input data is provided. 
+If a pre-existing distance matrix is provided, the distance computation step is 
+skipped entirely and `n_jobs` has no effect.
+
 ### Utilizing Numba and Progress Bars
 
 It may be helpful to use just-in-time (JIT) compilation in the cases where a lot of 
@@ -117,10 +149,12 @@ print(scores)
 Numba must be installed if the above to use JIT compilation and improve the 
 speed of multiple calls to `LocalOutlierProbability()`, and PyNomaly has been 
 tested with Numba version 0.45.1. An example of the speed difference that can 
-be realized with using Numba is avaialble in `examples/numba_speed_diff.py`. 
+be realized with using Numba is available in `examples/numba_speed_diff.py`. 
 
-You may also choose to print progress bars _with our without_ the use of numba 
+You may also choose to print progress bars _with or without_ the use of numba 
 by passing `progress_bar=True` to the `LocalOutlierProbability()` method as above.
+Progress bars are supported across all execution modes (sequential, parallel, 
+and Numba).
 
 ### Choosing Parameters
 
diff --git a/tests/test_loop.py b/tests/test_loop.py
index b9f4aaf..4458e91 100644
--- a/tests/test_loop.py
+++ b/tests/test_loop.py
@@ -827,3 +827,53 @@ def test_distance_matrix_consistency(X_n120) -> None:
 
     # Compare scores allowing for minor floating-point differences
     assert_array_almost_equal(scores_data, scores_dist, decimal=10, err_msg="Inconsistent LoOP scores due to self-distances")
+
+
+def test_n_jobs_equivalence(X_n140_outliers) -> None:
+    """
+    Tests that n_jobs > 1 produces equivalent results to n_jobs=1 when
+    using cluster labels (multiple clusters processed in parallel).
+    """
+    a = [0] * 120
+    b = [1] * 20
+    cluster_labels = a + b
+
+    clf_seq = loop.LocalOutlierProbability(
+        X_n140_outliers, n_neighbors=10, cluster_labels=cluster_labels, n_jobs=1
+    )
+    scores_seq = clf_seq.fit().local_outlier_probabilities
+
+    clf_par = loop.LocalOutlierProbability(
+        X_n140_outliers, n_neighbors=10, cluster_labels=cluster_labels, n_jobs=2
+    )
+    scores_par = clf_par.fit().local_outlier_probabilities
+
+    assert_array_almost_equal(scores_seq, scores_par, decimal=10)
+
+
+def test_n_jobs_single_cluster(X_n120) -> None:
+    """
+    Tests that n_jobs=-1 works correctly with a single cluster (falls back
+    to sequential since there is only one cluster to process).
+    """
+    clf1 = loop.LocalOutlierProbability(X_n120, n_neighbors=10, n_jobs=1)
+    scores1 = clf1.fit().local_outlier_probabilities
+
+    clf2 = loop.LocalOutlierProbability(X_n120, n_neighbors=10, n_jobs=-1)
+    scores2 = clf2.fit().local_outlier_probabilities
+
+    assert_array_almost_equal(scores1, scores2, decimal=10)
+
+
+def test_n_jobs_invalid() -> None:
+    """
+    Tests that invalid n_jobs values produce a warning and default to 1.
+    """
+    X = np.array([[1, 2], [3, 4], [5, 6], [7, 8]])
+
+    with pytest.warns(UserWarning) as record:
+        clf = loop.LocalOutlierProbability(X, n_neighbors=2, n_jobs=0)
+
+    messages = [r.message.args[0] for r in record]
+    assert any("n_jobs must be -1 or a positive integer" in m for m in messages)
+    assert clf.n_jobs == 1

From 8dd865a0a9c6bc01785a47ba177778b8f08cc86c Mon Sep 17 00:00:00 2001
From: vc1492a <vc1492a@gmail.com>
Date: Fri, 20 Mar 2026 11:16:37 -0700
Subject: [PATCH 2/5] chore: bump version to 0.4.0 and update changelog

Update version across loop.py, setup.py, and README badge.
Add changelog entry documenting all new features and improvements.

Made-with: Cursor
---
 PyNomaly/loop.py |  2 +-
 README.md        |  2 +-
 changelog.md     | 21 +++++++++++++++++++++
 setup.py         |  4 ++--
 4 files changed, 25 insertions(+), 4 deletions(-)

diff --git a/PyNomaly/loop.py b/PyNomaly/loop.py
index 46bdf31..d3dcc84 100644
--- a/PyNomaly/loop.py
+++ b/PyNomaly/loop.py
@@ -78,7 +78,7 @@ def _numba_distance_kernel_sequential(clust_points_vector, n_neighbors):
     pass
 
 __author__ = "Valentino Constantinou"
-__version__ = "0.3.5"
+__version__ = "0.4.0"
 __license__ = "Apache License, Version 2.0"
 
 
diff --git a/README.md b/README.md
index 4faab28..261a6bd 100644
--- a/README.md
+++ b/README.md
@@ -7,7 +7,7 @@ scores in the range of [0,1] that are directly interpretable as the probability
 PyNomaly is a core library of [deepchecks](https://github.com/deepchecks/deepchecks), [OmniDocBench](https://github.com/opendatalab/OmniDocBench) and [pysad](https://github.com/selimfirat/pysad). 
 
 [![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](https://opensource.org/licenses/Apache-2.0)
-[![PyPi](https://img.shields.io/badge/pypi-0.3.5-blue.svg)](https://pypi.python.org/pypi/PyNomaly/0.3.5)
+[![PyPi](https://img.shields.io/badge/pypi-0.4.0-blue.svg)](https://pypi.python.org/pypi/PyNomaly/0.4.0)
 [![Total Downloads](https://static.pepy.tech/badge/pynomaly)](https://pepy.tech/projects/pynomaly)
 [![Monthly Downloads](https://static.pepy.tech/badge/pynomaly/month)](https://pepy.tech/projects/pynomaly)
 ![Tests](https://github.com/vc1492a/PyNomaly/actions/workflows/tests.yml/badge.svg)
diff --git a/changelog.md b/changelog.md
index ee92cdd..d845315 100644
--- a/changelog.md
+++ b/changelog.md
@@ -4,6 +4,27 @@ All notable changes to PyNomaly will be documented in this Changelog.
 The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/) 
 and adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html).
 
+## 0.4.0
+### Added
+- Parallel distance computation via the new `n_jobs` parameter. Set `n_jobs=-1` 
+to use all available CPU cores for cross-cluster multiprocessing, or specify a 
+positive integer for a fixed number of workers 
+([Issue #36](https://github.com/vc1492a/PyNomaly/issues/36)).
+- Numba `prange`-based parallel kernels for thread-level parallelism when 
+`use_numba=True` and `n_jobs > 1`.
+- Optional `scipy` acceleration: uses `scipy.spatial.distance.cdist` for 
+distance computation and `scipy.special.erf` for the error function when 
+scipy is available, with graceful fallback to pure NumPy.
+### Changed
+- Replaced the O(n^2) Python nested loop for distance computation with a 
+vectorized NumPy implementation using chunked broadcasting. Progress bar 
+support is preserved via chunk-level updates.
+- Restructured the Numba JIT path from a generator-based approach (incompatible 
+with Numba's parallel mode) to non-generator kernels that support `numba.prange`.
+- Vectorized `_standard_distances`, `_prob_distances`, and 
+`_norm_prob_outlier_factor` pipeline methods, replacing Python `for` loops 
+with NumPy array operations.
+
 ## 0.3.5
 ### Changed
 - Refactored the `Validate` class by dissolving it and moving validation methods 
diff --git a/setup.py b/setup.py
index d58b87e..a73829b 100644
--- a/setup.py
+++ b/setup.py
@@ -7,7 +7,7 @@
 setup(
     name='PyNomaly',
     packages=['PyNomaly'],
-    version='0.3.5',
+    version='0.4.0',
     description='A Python 3 implementation of LoOP: Local Outlier '
                 'Probabilities, a local density based outlier detection '
                 'method providing an outlier score in the range of [0,1].',
@@ -16,7 +16,7 @@
     long_description=long_description,
     long_description_content_type='text/markdown',
     url='https://github.com/vc1492a/PyNomaly',
-    download_url='https://github.com/vc1492a/PyNomaly/archive/0.3.5.tar.gz',
+    download_url='https://github.com/vc1492a/PyNomaly/archive/0.4.0.tar.gz',
     keywords=['outlier', 'anomaly', 'detection', 'machine', 'learning',
               'probability'],
     classifiers=[],

From 1de847fff4c5fd1d560ade2a72acf0f887fec16b Mon Sep 17 00:00:00 2001
From: vc1492a <vc1492a@gmail.com>
Date: Tue, 26 May 2026 11:12:17 -0700
Subject: [PATCH 3/5] test: improve coverage for parallel/vectorized paths and
 update benchmark

- Add tests for 1D data vectorized path, parallel progress bar,
  n_jobs=-2 validation, and single-cluster progress bar
- Update benchmark script with n_jobs parallel examples and
  __main__ guard for macOS multiprocessing compatibility
- Note n_jobs replacing num_threads in changelog

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 changelog.md                 |  6 +++
 examples/numba_speed_diff.py | 92 +++++++++++++++++++++++++-----------
 tests/test_loop.py           | 58 +++++++++++++++++++++++
 3 files changed, 129 insertions(+), 27 deletions(-)

diff --git a/changelog.md b/changelog.md
index d845315..ee0fe82 100644
--- a/changelog.md
+++ b/changelog.md
@@ -24,6 +24,12 @@ with Numba's parallel mode) to non-generator kernels that support `numba.prange`
 - Vectorized `_standard_distances`, `_prob_distances`, and 
 `_norm_prob_outlier_factor` pipeline methods, replacing Python `for` loops 
 with NumPy array operations.
+- Replaced the previous `num_threads` parameter (from the unreleased 
+`feature/numba_parallel` branch) with `n_jobs`, which provides a unified 
+interface for both process-level parallelism (via `concurrent.futures`) and 
+Numba thread-level parallelism (via `prange`). The `n_jobs` parameter follows 
+the scikit-learn convention (`-1` for all cores, positive integer for a fixed 
+number of workers).
 
 ## 0.3.5
 ### Changed
diff --git a/examples/numba_speed_diff.py b/examples/numba_speed_diff.py
index 2d4f622..cc36ba7 100644
--- a/examples/numba_speed_diff.py
+++ b/examples/numba_speed_diff.py
@@ -2,30 +2,68 @@
 from PyNomaly import loop
 import time
 
-# generate a large set of data
-data = np.ones(shape=(10000, 4))
-
-# first time the process without Numba
-# use the progress bar to track progress
-
-t1 = time.time()
-scores_numpy = loop.LocalOutlierProbability(
-    data,
-    n_neighbors=3,
-    use_numba=False,
-    progress_bar=True
-).fit().local_outlier_probabilities
-t2 = time.time()
-seconds_no_numba = t2 - t1
-print("\nComputation took " + str(seconds_no_numba) + " seconds without Numba JIT.")
-
-t3 = time.time()
-scores_numba = loop.LocalOutlierProbability(
-    data,
-    n_neighbors=3,
-    use_numba=True,
-    progress_bar=True
-).fit().local_outlier_probabilities
-t4 = time.time()
-seconds_numba = t4 - t3
-print("\nComputation took " + str(seconds_numba) + " seconds with Numba JIT.")
+
+def main():
+    # generate a large set of data
+    data = np.ones(shape=(10000, 4))
+
+    # Vectorized NumPy (default, no Numba)
+    t1 = time.time()
+    scores_numpy = loop.LocalOutlierProbability(
+        data,
+        n_neighbors=3,
+        use_numba=False,
+        progress_bar=True
+    ).fit().local_outlier_probabilities
+    t2 = time.time()
+    seconds_no_numba = t2 - t1
+    print("\nComputation took " + str(seconds_no_numba) + " seconds without Numba JIT.")
+
+    # Numba JIT (sequential)
+    t3 = time.time()
+    scores_numba = loop.LocalOutlierProbability(
+        data,
+        n_neighbors=3,
+        use_numba=True,
+        progress_bar=True
+    ).fit().local_outlier_probabilities
+    t4 = time.time()
+    seconds_numba = t4 - t3
+    print("\nComputation took " + str(seconds_numba) + " seconds with Numba JIT.")
+
+    # Multi-cluster parallel example
+    np.random.seed(42)
+    cluster_a = np.random.randn(5000, 4) + 0
+    cluster_b = np.random.randn(5000, 4) + 10
+    multi_data = np.vstack([cluster_a, cluster_b])
+    cluster_labels = [0] * 5000 + [1] * 5000
+
+    # Sequential (n_jobs=1) with clusters
+    t5 = time.time()
+    scores_seq = loop.LocalOutlierProbability(
+        multi_data,
+        n_neighbors=3,
+        cluster_labels=cluster_labels,
+        n_jobs=1,
+        progress_bar=True
+    ).fit().local_outlier_probabilities
+    t6 = time.time()
+    seconds_seq = t6 - t5
+    print("\nComputation took " + str(seconds_seq) + " seconds with n_jobs=1 (2 clusters).")
+
+    # Parallel (n_jobs=-1) with clusters
+    t7 = time.time()
+    scores_par = loop.LocalOutlierProbability(
+        multi_data,
+        n_neighbors=3,
+        cluster_labels=cluster_labels,
+        n_jobs=-1,
+        progress_bar=True
+    ).fit().local_outlier_probabilities
+    t8 = time.time()
+    seconds_par = t8 - t7
+    print("\nComputation took " + str(seconds_par) + " seconds with n_jobs=-1 (2 clusters).")
+
+
+if __name__ == '__main__':
+    main()
diff --git a/tests/test_loop.py b/tests/test_loop.py
index 4458e91..ec02bf4 100644
--- a/tests/test_loop.py
+++ b/tests/test_loop.py
@@ -829,6 +829,64 @@ def test_distance_matrix_consistency(X_n120) -> None:
     assert_array_almost_equal(scores_data, scores_dist, decimal=10, err_msg="Inconsistent LoOP scores due to self-distances")
 
 
+def test_vectorized_1d_data() -> None:
+    """
+    Tests the vectorized distance path with 1-dimensional data, exercising
+    the ndim == 1 reshape branch in _distances_vectorized.
+    """
+    X = np.array([1.0, 2.0, 3.0, 10.0, 11.0, 12.0, 50.0])
+    clf = loop.LocalOutlierProbability(X, n_neighbors=3)
+    scores = clf.fit().local_outlier_probabilities
+    assert scores is not None
+    assert len(scores) == len(X)
+    assert scores[-1] > 0
+
+
+def test_parallel_with_progress_bar(X_n140_outliers) -> None:
+    """
+    Tests that parallel processing works with the progress bar enabled,
+    covering the progress_bar branch inside _distances_parallel.
+    """
+    a = [0] * 120
+    b = [1] * 20
+    cluster_labels = a + b
+
+    clf = loop.LocalOutlierProbability(
+        X_n140_outliers, n_neighbors=10, cluster_labels=cluster_labels,
+        n_jobs=2, progress_bar=True
+    )
+    scores = clf.fit().local_outlier_probabilities
+    assert scores is not None
+    assert len(scores) == len(X_n140_outliers)
+
+
+def test_n_jobs_negative_two() -> None:
+    """
+    Tests that n_jobs=-2 (invalid) produces a warning and defaults to 1.
+    """
+    X = np.array([[1, 2], [3, 4], [5, 6], [7, 8]])
+
+    with pytest.warns(UserWarning) as record:
+        clf = loop.LocalOutlierProbability(X, n_neighbors=2, n_jobs=-2)
+
+    messages = [r.message.args[0] for r in record]
+    assert any("n_jobs must be -1 or a positive integer" in m for m in messages)
+    assert clf.n_jobs == 1
+
+
+def test_vectorized_progress_bar_single_cluster(X_n120) -> None:
+    """
+    Tests the progress bar within the vectorized (non-parallel) distance path
+    on single-cluster data.
+    """
+    clf = loop.LocalOutlierProbability(
+        X_n120, n_neighbors=10, n_jobs=1, progress_bar=True
+    )
+    scores = clf.fit().local_outlier_probabilities
+    assert scores is not None
+    assert len(scores) == len(X_n120)
+
+
 def test_n_jobs_equivalence(X_n140_outliers) -> None:
     """
     Tests that n_jobs > 1 produces equivalent results to n_jobs=1 when

From bf2efc7f732866e451cc852605c45fc137671cb0 Mon Sep 17 00:00:00 2001
From: vc1492a <vc1492a@gmail.com>
Date: Tue, 26 May 2026 11:30:26 -0700
Subject: [PATCH 4/5] fix: enable Numba prange for single-cluster data and add
 Numba tests

Numba prange dispatch was incorrectly capped to the number of clusters,
preventing intra-cluster thread parallelism for single-cluster data.
This fix provides 2-3x speedups on multi-core machines. Adds dedicated
Numba tests, a parallel benchmark script, and updated documentation.

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 PyNomaly/loop.py               |   5 +-
 README.md                      |  26 +++--
 changelog.md                   |   7 ++
 examples/parallel_benchmark.py | 191 +++++++++++++++++++++++++++++++++
 tests/test_loop.py             |  83 ++++++++++++++
 5 files changed, 303 insertions(+), 9 deletions(-)
 create mode 100644 examples/parallel_benchmark.py

diff --git a/PyNomaly/loop.py b/PyNomaly/loop.py
index d3dcc84..711fb10 100644
--- a/PyNomaly/loop.py
+++ b/PyNomaly/loop.py
@@ -807,12 +807,15 @@ def _distances(self, progress_bar: bool = False) -> None:
         n_jobs = self.n_jobs
         if n_jobs == -1:
             n_jobs = os.cpu_count() or 1
+        user_requested_parallel = n_jobs > 1
         n_jobs = min(n_jobs, len(clusters))
 
         if self.use_numba:
+            # Numba prange parallelizes within each cluster (over observations),
+            # so it benefits even with a single cluster.
             self._distances_numba(
                 clusters, distances, indexes, progress_bar,
-                parallel=(n_jobs > 1)
+                parallel=user_requested_parallel
             )
         elif n_jobs > 1 and len(clusters) > 1:
             self._distances_parallel(
diff --git a/README.md b/README.md
index 261a6bd..2753acc 100644
--- a/README.md
+++ b/README.md
@@ -123,11 +123,18 @@ print(scores)
 
 Set `n_jobs=-1` to use all available CPU cores, or specify a positive integer 
 for a fixed number of workers. The default value of `n_jobs=1` runs sequentially. 
-Parallelism is most beneficial when the data contains multiple clusters. For 
-single-cluster data, `n_jobs` has no effect as there is only one unit of work.
+
+Without Numba, `n_jobs` controls multiprocessing across clusters via 
+`ProcessPoolExecutor`. This is most beneficial for data with many clusters, 
+though process startup overhead means the speedup is modest and mainly appears 
+at larger scales (e.g. 20,000+ observations across 4+ clusters).
 
 When using Numba (`use_numba=True`) with `n_jobs > 1`, PyNomaly uses Numba's 
-thread-level parallelism (`prange`) instead of multiprocessing.
+thread-level parallelism (`prange`) instead of multiprocessing. This 
+parallelizes the distance computation *within* each cluster, providing 
+significant speedups (2-3x on 8 cores) even with a single cluster.
+For best performance with large datasets, use `use_numba=True` with 
+`n_jobs=-1`.
 
 Note that parallel processing is only applicable when raw input data is provided. 
 If a pre-existing distance matrix is provided, the distance computation step is 
@@ -146,10 +153,11 @@ scores = m.local_outlier_probabilities
 print(scores)
 ```
 
-Numba must be installed if the above to use JIT compilation and improve the 
-speed of multiple calls to `LocalOutlierProbability()`, and PyNomaly has been 
-tested with Numba version 0.45.1. An example of the speed difference that can 
-be realized with using Numba is available in `examples/numba_speed_diff.py`. 
+Numba must be installed to use JIT compilation and improve the 
+speed of multiple calls to `LocalOutlierProbability()`. PyNomaly has been 
+tested with Numba versions 0.45.1 through 0.65.1. An example of the speed 
+difference that can be realized with using Numba is available in 
+`examples/numba_speed_diff.py`. 
 
 You may also choose to print progress bars _with or without_ the use of numba 
 by passing `progress_bar=True` to the `LocalOutlierProbability()` method as above.
@@ -465,7 +473,9 @@ python3 -m pytest --cov=PyNomaly -s -v
 
 To run the tests with Numba enabled, simply set the flag `NUMBA` in `test_loop.py` 
 to `True`. Note that a drop in coverage is expected due to portions of the code 
-being compiled upon code execution. 
+being compiled upon code execution. Additionally, there are dedicated 
+Numba-specific tests (`test_numba_*`) that run automatically when Numba is 
+installed and are skipped otherwise.
 
 ## Versioning
 [Semantic versioning](http://semver.org/) is used for this project. If contributing, please conform to semantic
diff --git a/changelog.md b/changelog.md
index ee0fe82..69a2ad2 100644
--- a/changelog.md
+++ b/changelog.md
@@ -15,6 +15,13 @@ positive integer for a fixed number of workers
 - Optional `scipy` acceleration: uses `scipy.spatial.distance.cdist` for 
 distance computation and `scipy.special.erf` for the error function when 
 scipy is available, with graceful fallback to pure NumPy.
+- Dedicated Numba-specific tests (`test_numba_*`) that automatically run when 
+Numba is installed and are skipped otherwise, covering sequential equivalence, 
+parallel equivalence, and progress bar integration with Numba.
+### Fixed
+- Numba `prange` dispatch now correctly activates for single-cluster data 
+when `n_jobs > 1`, providing 2-3x speedups that were previously blocked by 
+incorrect capping of the parallelism flag to the number of clusters.
 ### Changed
 - Replaced the O(n^2) Python nested loop for distance computation with a 
 vectorized NumPy implementation using chunked broadcasting. Progress bar 
diff --git a/examples/parallel_benchmark.py b/examples/parallel_benchmark.py
new file mode 100644
index 0000000..b71a41d
--- /dev/null
+++ b/examples/parallel_benchmark.py
@@ -0,0 +1,191 @@
+"""
+Thorough benchmark of PyNomaly's parallel processing capabilities.
+
+Tests across multiple dimensions:
+  - Data sizes (total observations)
+  - Number of clusters (parallelism granularity)
+  - Execution modes: vectorized, multiprocessing (n_jobs), Numba sequential, Numba parallel
+  - Number of workers (1, 2, 4, all cores)
+
+Run:
+    python examples/parallel_benchmark.py
+"""
+
+import numpy as np
+import os
+import time
+from PyNomaly import loop
+
+N_FEATURES = 4
+N_NEIGHBORS = 10
+N_REPEATS = 3
+
+try:
+    import numba
+    HAS_NUMBA = True
+except ImportError:
+    HAS_NUMBA = False
+
+
+def generate_clustered_data(n_per_cluster, n_clusters, n_features=N_FEATURES, seed=42):
+    rng = np.random.RandomState(seed)
+    clusters = []
+    labels = []
+    for i in range(n_clusters):
+        center = rng.uniform(-20, 20, size=n_features)
+        points = rng.randn(n_per_cluster, n_features) + center
+        clusters.append(points)
+        labels.extend([i] * n_per_cluster)
+    data = np.vstack(clusters)
+    return data, labels
+
+
+def time_fit(data, cluster_labels, n_neighbors, use_numba, n_jobs, warmup=False):
+    """Time a single fit call. Returns elapsed seconds."""
+    clf = loop.LocalOutlierProbability(
+        data, n_neighbors=n_neighbors,
+        cluster_labels=cluster_labels,
+        use_numba=use_numba,
+        n_jobs=n_jobs
+    )
+    t0 = time.perf_counter()
+    clf.fit()
+    elapsed = time.perf_counter() - t0
+    return elapsed
+
+
+def run_benchmark(data, labels, n_neighbors, label, configs):
+    """Run a set of configurations and print results."""
+    n_total = len(data)
+    n_clusters = len(set(labels))
+
+    print(f"\n{'='*70}")
+    print(f"  {label}")
+    print(f"  {n_total:,} observations, {n_clusters} clusters, "
+          f"{n_total // n_clusters:,} per cluster, {data.shape[1]} features")
+    print(f"{'='*70}")
+    print(f"  {'Mode':<35} {'Best':>8} {'Mean':>8} {'Speedup':>8}")
+    print(f"  {'-'*35} {'-'*8} {'-'*8} {'-'*8}")
+
+    baseline_best = None
+    results = []
+
+    for mode_name, use_numba, n_jobs in configs:
+        if use_numba and not HAS_NUMBA:
+            continue
+
+        # Warmup for Numba JIT compilation
+        if use_numba:
+            small_data, small_labels = generate_clustered_data(50, n_clusters, seed=99)
+            time_fit(small_data, small_labels, min(n_neighbors, 10), use_numba, n_jobs, warmup=True)
+
+        times = []
+        for _ in range(N_REPEATS):
+            t = time_fit(data, labels, n_neighbors, use_numba, n_jobs)
+            times.append(t)
+
+        best = min(times)
+        mean = sum(times) / len(times)
+
+        if baseline_best is None:
+            baseline_best = best
+
+        speedup = baseline_best / best if best > 0 else float('inf')
+        results.append((mode_name, best, mean, speedup))
+        print(f"  {mode_name:<35} {best:>7.3f}s {mean:>7.3f}s {speedup:>7.2f}x")
+
+    return results
+
+
+def main():
+    cpu_count = os.cpu_count() or 1
+    print(f"PyNomaly Parallel Processing Benchmark")
+    print(f"CPU cores: {cpu_count}")
+    print(f"Numba available: {HAS_NUMBA}")
+    print(f"Repeats per config: {N_REPEATS}")
+
+    # Scenario 1: Small data, many clusters — tests multiprocessing dispatch overhead
+    data, labels = generate_clustered_data(500, 8)
+    run_benchmark(data, labels, N_NEIGHBORS,
+        "Scenario 1: Small clusters (500 pts x 8 clusters = 4,000 total)",
+        [
+            ("Vectorized (baseline)",       False, 1),
+            ("Multiprocessing n_jobs=2",     False, 2),
+            ("Multiprocessing n_jobs=4",     False, 4),
+            ("Multiprocessing n_jobs=-1",    False, -1),
+            ("Numba sequential",             True,  1),
+            ("Numba parallel (prange)",      True,  -1),
+        ]
+    )
+
+    # Scenario 2: Medium data, moderate clusters — practical use case
+    data, labels = generate_clustered_data(2000, 4)
+    run_benchmark(data, labels, N_NEIGHBORS,
+        "Scenario 2: Medium clusters (2,000 pts x 4 clusters = 8,000 total)",
+        [
+            ("Vectorized (baseline)",       False, 1),
+            ("Multiprocessing n_jobs=2",     False, 2),
+            ("Multiprocessing n_jobs=4",     False, 4),
+            ("Multiprocessing n_jobs=-1",    False, -1),
+            ("Numba sequential",             True,  1),
+            ("Numba parallel (prange)",      True,  -1),
+        ]
+    )
+
+    # Scenario 3: Large data, many clusters — where parallelism should help most
+    data, labels = generate_clustered_data(2000, 8)
+    run_benchmark(data, labels, N_NEIGHBORS,
+        "Scenario 3: Large + many clusters (2,000 pts x 8 clusters = 16,000 total)",
+        [
+            ("Vectorized (baseline)",       False, 1),
+            ("Multiprocessing n_jobs=2",     False, 2),
+            ("Multiprocessing n_jobs=4",     False, 4),
+            ("Multiprocessing n_jobs=8",     False, 8),
+            ("Multiprocessing n_jobs=-1",    False, -1),
+            ("Numba sequential",             True,  1),
+            ("Numba parallel (prange)",      True,  -1),
+        ]
+    )
+
+    # Scenario 4: Scaling test — fixed cluster size, increasing cluster count
+    print(f"\n{'='*70}")
+    print(f"  Scenario 4: Scaling test — 1,000 pts/cluster, increasing clusters")
+    print(f"{'='*70}")
+    print(f"  {'Clusters':>8} {'Total':>8} {'Seq':>8} {'Par':>8} {'Speedup':>8}")
+    print(f"  {'-'*8} {'-'*8} {'-'*8} {'-'*8} {'-'*8}")
+
+    for n_clusters in [2, 4, 8, 16]:
+        data, labels = generate_clustered_data(1000, n_clusters)
+
+        seq_times = []
+        par_times = []
+        for _ in range(N_REPEATS):
+            seq_times.append(time_fit(data, labels, N_NEIGHBORS, False, 1))
+            par_times.append(time_fit(data, labels, N_NEIGHBORS, False, -1))
+
+        seq_best = min(seq_times)
+        par_best = min(par_times)
+        speedup = seq_best / par_best if par_best > 0 else float('inf')
+        total = n_clusters * 1000
+
+        print(f"  {n_clusters:>8} {total:>7,} {seq_best:>7.3f}s {par_best:>7.3f}s {speedup:>7.2f}x")
+
+    # Scenario 5: Single large cluster — no parallelism benefit expected
+    data, labels = generate_clustered_data(10000, 1)
+    run_benchmark(data, labels, N_NEIGHBORS,
+        "Scenario 5: Single cluster (10,000 pts) — no multi-cluster parallelism",
+        [
+            ("Vectorized (baseline)",       False, 1),
+            ("Multiprocessing n_jobs=-1",    False, -1),
+            ("Numba sequential",             True,  1),
+            ("Numba parallel (prange)",      True,  -1),
+        ]
+    )
+
+    print(f"\n{'='*70}")
+    print("  Benchmark complete.")
+    print(f"{'='*70}")
+
+
+if __name__ == '__main__':
+    main()
diff --git a/tests/test_loop.py b/tests/test_loop.py
index ec02bf4..5a7106a 100644
--- a/tests/test_loop.py
+++ b/tests/test_loop.py
@@ -935,3 +935,86 @@ def test_n_jobs_invalid() -> None:
     messages = [r.message.args[0] for r in record]
     assert any("n_jobs must be -1 or a positive integer" in m for m in messages)
     assert clf.n_jobs == 1
+
+
+# --- Numba-specific tests ---
+# These test the Numba code paths directly, skipped if Numba is not installed.
+
+_has_numba = False
+try:
+    import numba
+    _has_numba = True
+except ImportError:
+    pass
+
+
+@pytest.mark.skipif(not _has_numba, reason="Numba not installed")
+def test_numba_sequential_equivalence(X_n8) -> None:
+    """
+    Tests that use_numba=True with n_jobs=1 produces equivalent results to
+    the default vectorized path.
+    """
+    clf_vec = loop.LocalOutlierProbability(X_n8, n_neighbors=5, use_numba=False)
+    scores_vec = clf_vec.fit().local_outlier_probabilities
+
+    clf_numba = loop.LocalOutlierProbability(X_n8, n_neighbors=5, use_numba=True)
+    scores_numba = clf_numba.fit().local_outlier_probabilities
+
+    assert_array_almost_equal(scores_vec, scores_numba, decimal=6)
+
+
+@pytest.mark.skipif(not _has_numba, reason="Numba not installed")
+def test_numba_parallel_equivalence(X_n140_outliers) -> None:
+    """
+    Tests that Numba parallel mode (use_numba=True, n_jobs > 1) produces
+    equivalent results to Numba sequential mode on clustered data.
+    """
+    a = [0] * 120
+    b = [1] * 20
+    cluster_labels = a + b
+
+    clf_seq = loop.LocalOutlierProbability(
+        X_n140_outliers, n_neighbors=10, cluster_labels=cluster_labels,
+        use_numba=True, n_jobs=1
+    )
+    scores_seq = clf_seq.fit().local_outlier_probabilities
+
+    clf_par = loop.LocalOutlierProbability(
+        X_n140_outliers, n_neighbors=10, cluster_labels=cluster_labels,
+        use_numba=True, n_jobs=2
+    )
+    scores_par = clf_par.fit().local_outlier_probabilities
+
+    assert_array_almost_equal(scores_seq, scores_par, decimal=10)
+
+
+@pytest.mark.skipif(not _has_numba, reason="Numba not installed")
+def test_numba_with_progress_bar(X_n120) -> None:
+    """
+    Tests that the Numba path works with the progress bar enabled.
+    """
+    clf = loop.LocalOutlierProbability(
+        X_n120, n_neighbors=10, use_numba=True, progress_bar=True
+    )
+    scores = clf.fit().local_outlier_probabilities
+    assert scores is not None
+    assert len(scores) == len(X_n120)
+
+
+@pytest.mark.skipif(not _has_numba, reason="Numba not installed")
+def test_numba_prange_single_cluster(X_n120) -> None:
+    """
+    Tests that Numba prange is used for single-cluster data when n_jobs > 1,
+    and produces equivalent results to the sequential path.
+    """
+    clf_seq = loop.LocalOutlierProbability(
+        X_n120, n_neighbors=10, use_numba=True, n_jobs=1
+    )
+    scores_seq = clf_seq.fit().local_outlier_probabilities
+
+    clf_par = loop.LocalOutlierProbability(
+        X_n120, n_neighbors=10, use_numba=True, n_jobs=-1
+    )
+    scores_par = clf_par.fit().local_outlier_probabilities
+
+    assert_array_almost_equal(scores_seq, scores_par, decimal=6)

From 60913098c1655df9604df24917210f2144834698 Mon Sep 17 00:00:00 2001
From: vc1492a <vc1492a@gmail.com>
Date: Tue, 26 May 2026 11:34:47 -0700
Subject: [PATCH 5/5] refactor: remove multiprocessing, use Numba prange for
 all parallelism

Multiprocessing via ProcessPoolExecutor was consistently slower than
sequential due to process startup and data serialization overhead.
Numba prange provides 2-3x speedups with zero overhead. n_jobs now
only takes effect with use_numba=True; a warning is issued otherwise.

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 PyNomaly/loop.py               |  83 +++-------------------
 README.md                      |  73 +++++++++-----------
 changelog.md                   |  26 +++----
 examples/numba_speed_diff.py   |  36 +++-------
 examples/parallel_benchmark.py | 122 ++++++++++++---------------------
 tests/test_loop.py             |  57 +++------------
 6 files changed, 113 insertions(+), 284 deletions(-)

diff --git a/PyNomaly/loop.py b/PyNomaly/loop.py
index 711fb10..31a32b1 100644
--- a/PyNomaly/loop.py
+++ b/PyNomaly/loop.py
@@ -1,4 +1,3 @@
-from concurrent.futures import ProcessPoolExecutor, as_completed
 from math import erf, sqrt
 import numpy as np
 import os
@@ -82,40 +81,6 @@ def _numba_distance_kernel_sequential(clust_points_vector, n_neighbors):
 __license__ = "Apache License, Version 2.0"
 
 
-def _cluster_distances_worker(args):
-    """Top-level worker for ProcessPoolExecutor (must be picklable)."""
-    clust_points_vector, global_indices, n_neighbors = args
-    n = clust_points_vector.shape[0]
-
-    if clust_points_vector.ndim == 1:
-        clust_points_vector = clust_points_vector.reshape(-1, 1)
-
-    local_distances = np.full((n, n_neighbors), 9e10, dtype=float)
-    local_indexes = np.full((n, n_neighbors), 9e10, dtype=float)
-
-    chunk_size = min(256, n)
-    for chunk_start in range(0, n, chunk_size):
-        chunk_end = min(chunk_start + chunk_size, n)
-        chunk = clust_points_vector[chunk_start:chunk_end]
-
-        if _scipy_cdist is not None:
-            dist = _scipy_cdist(chunk, clust_points_vector, metric="euclidean")
-        else:
-            diff = chunk[:, np.newaxis, :] - clust_points_vector[np.newaxis, :, :]
-            dist = np.sqrt((diff ** 2).sum(axis=2))
-
-        row_idx = np.arange(chunk_end - chunk_start)
-        dist[row_idx, row_idx + chunk_start] = np.inf
-
-        knn_idx = np.argpartition(dist, n_neighbors, axis=1)[:, :n_neighbors]
-        knn_dists = np.take_along_axis(dist, knn_idx, axis=1)
-
-        local_distances[chunk_start:chunk_end] = knn_dists
-        local_indexes[chunk_start:chunk_end] = global_indices[knn_idx]
-
-    return local_distances, local_indexes, global_indices
-
-
 # Custom Exceptions
 class PyNomalyError(Exception):
     """Base exception for PyNomaly."""
@@ -175,9 +140,9 @@ class LocalOutlierProbability(object):
     sample (optional, default 10)
     :param cluster_labels: a numpy array of cluster assignments w.r.t. each 
     sample (optional, default None)
-    :param n_jobs: the number of parallel workers for distance computation.
-    Use -1 to use all available CPU cores, or 1 for sequential processing
-    (optional, default 1)
+    :param n_jobs: controls Numba thread-level parallelism via prange.
+    Use -1 to use all available CPU cores, or 1 for sequential processing.
+    Only effective when use_numba=True (optional, default 1)
     :return:
     """ """
 
@@ -751,31 +716,6 @@ def _distances_numba(
                     progress, idx + 1, len(clusters)
                 )
 
-    def _distances_parallel(
-        self, clusters, distances, indexes, n_jobs, progress_bar
-    ) -> None:
-        """Parallel distance computation across clusters via multiprocessing."""
-        progress = "="
-        worker_args = [
-            (cv, gi, self.n_neighbors) for cv, gi in clusters
-        ]
-
-        with ProcessPoolExecutor(max_workers=n_jobs) as executor:
-            futures = {
-                executor.submit(_cluster_distances_worker, args): idx
-                for idx, args in enumerate(worker_args)
-            }
-            completed_clusters = 0
-            for future in as_completed(futures):
-                local_dists, local_idxs, global_indices = future.result()
-                distances[global_indices] = local_dists
-                indexes[global_indices] = local_idxs
-                completed_clusters += 1
-                if progress_bar:
-                    progress = Utils.emit_progress_bar(
-                        progress, completed_clusters, len(clusters)
-                    )
-
     def _distances(self, progress_bar: bool = False) -> None:
         """
         Provides the distances between each observation and it's closest
@@ -807,21 +747,20 @@ def _distances(self, progress_bar: bool = False) -> None:
         n_jobs = self.n_jobs
         if n_jobs == -1:
             n_jobs = os.cpu_count() or 1
-        user_requested_parallel = n_jobs > 1
-        n_jobs = min(n_jobs, len(clusters))
 
         if self.use_numba:
-            # Numba prange parallelizes within each cluster (over observations),
-            # so it benefits even with a single cluster.
             self._distances_numba(
                 clusters, distances, indexes, progress_bar,
-                parallel=user_requested_parallel
-            )
-        elif n_jobs > 1 and len(clusters) > 1:
-            self._distances_parallel(
-                clusters, distances, indexes, n_jobs, progress_bar
+                parallel=(n_jobs > 1)
             )
         else:
+            if n_jobs > 1:
+                warnings.warn(
+                    "n_jobs > 1 requires use_numba=True for parallel "
+                    "processing. Install Numba and set use_numba=True "
+                    "to enable parallelism. Falling back to sequential.",
+                    UserWarning,
+                )
             self._distances_vectorized(
                 clusters, distances, indexes, progress_bar
             )
diff --git a/README.md b/README.md
index 2753acc..8b46a45 100644
--- a/README.md
+++ b/README.md
@@ -103,66 +103,59 @@ normalization scheme prior to using LoOP, especially when working with multiple
 Users must also appropriately handle missing values prior to using LoOP, as LoOP does not support Pandas
 DataFrames or Numpy arrays with missing values.
 
-### Parallel Processing
+### Speeding Things Up with Numba
 
-PyNomaly supports parallel distance computation via the `n_jobs` parameter.
-When data includes multiple clusters (via `cluster_labels`), each cluster's 
-distance computation is independent and can be processed across multiple CPU cores:
+For large datasets, Numba's just-in-time (JIT) compilation can significantly 
+speed up distance computation. Enable it with `use_numba=True`:
+
+```python
+from PyNomaly import loop
+m = loop.LocalOutlierProbability(data, extent=2, n_neighbors=20, use_numba=True).fit()
+scores = m.local_outlier_probabilities
+print(scores)
+```
+
+To go further, set `n_jobs=-1` to enable Numba's thread-level parallelism 
+(`prange`), which distributes work across all available CPU cores:
 
 ```python
 from PyNomaly import loop
-from sklearn.cluster import DBSCAN
-db = DBSCAN(eps=0.6, min_samples=50).fit(data)
 m = loop.LocalOutlierProbability(
     data, extent=2, n_neighbors=20,
-    cluster_labels=list(db.labels_), n_jobs=-1
+    use_numba=True, n_jobs=-1
 ).fit()
 scores = m.local_outlier_probabilities
 print(scores)
 ```
 
-Set `n_jobs=-1` to use all available CPU cores, or specify a positive integer 
-for a fixed number of workers. The default value of `n_jobs=1` runs sequentially. 
-
-Without Numba, `n_jobs` controls multiprocessing across clusters via 
-`ProcessPoolExecutor`. This is most beneficial for data with many clusters, 
-though process startup overhead means the speedup is modest and mainly appears 
-at larger scales (e.g. 20,000+ observations across 4+ clusters).
+This provides **2-3x speedups** on multi-core machines (benchmarked on 8 cores). 
+The parallelism works within each cluster's distance computation, so it benefits 
+both single-cluster and multi-cluster data. Set `n_jobs=-1` to use all cores, 
+or specify a positive integer for a fixed number of threads. The default 
+`n_jobs=1` runs sequentially.
 
-When using Numba (`use_numba=True`) with `n_jobs > 1`, PyNomaly uses Numba's 
-thread-level parallelism (`prange`) instead of multiprocessing. This 
-parallelizes the distance computation *within* each cluster, providing 
-significant speedups (2-3x on 8 cores) even with a single cluster.
-For best performance with large datasets, use `use_numba=True` with 
-`n_jobs=-1`.
+Numba must be installed to use JIT compilation. PyNomaly has been tested with 
+Numba versions 0.45.1 through 0.65.1. An example of the speed difference is 
+available in `examples/numba_speed_diff.py`, and a comprehensive parallel 
+benchmark is in `examples/parallel_benchmark.py`.
 
-Note that parallel processing is only applicable when raw input data is provided. 
-If a pre-existing distance matrix is provided, the distance computation step is 
-skipped entirely and `n_jobs` has no effect.
+Note that `n_jobs` only takes effect when `use_numba=True`. If `n_jobs > 1` 
+is set without Numba, PyNomaly will warn and fall back to the sequential 
+vectorized path. Parallel processing is also only applicable when raw input 
+data is provided — if a pre-existing distance matrix is provided, the distance 
+computation step is skipped entirely.
 
-### Utilizing Numba and Progress Bars
+### Progress Bars
 
-It may be helpful to use just-in-time (JIT) compilation in the cases where a lot of 
-observations are scored. Numba, a JIT compiler for Python, may be used 
-with PyNomaly by setting `use_numba=True`:
+You may choose to print progress bars _with or without_ the use of Numba 
+by passing `progress_bar=True` to `LocalOutlierProbability()`:
 
 ```python
 from PyNomaly import loop
-m = loop.LocalOutlierProbability(data, extent=2, n_neighbors=20, use_numba=True, progress_bar=True).fit()
-scores = m.local_outlier_probabilities
-print(scores)
+m = loop.LocalOutlierProbability(data, use_numba=True, n_jobs=-1, progress_bar=True).fit()
 ```
 
-Numba must be installed to use JIT compilation and improve the 
-speed of multiple calls to `LocalOutlierProbability()`. PyNomaly has been 
-tested with Numba versions 0.45.1 through 0.65.1. An example of the speed 
-difference that can be realized with using Numba is available in 
-`examples/numba_speed_diff.py`. 
-
-You may also choose to print progress bars _with or without_ the use of numba 
-by passing `progress_bar=True` to the `LocalOutlierProbability()` method as above.
-Progress bars are supported across all execution modes (sequential, parallel, 
-and Numba).
+Progress bars are supported in both sequential and Numba execution modes.
 
 ### Choosing Parameters
 
diff --git a/changelog.md b/changelog.md
index 69a2ad2..9b41565 100644
--- a/changelog.md
+++ b/changelog.md
@@ -6,22 +6,16 @@ and adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html).
 
 ## 0.4.0
 ### Added
-- Parallel distance computation via the new `n_jobs` parameter. Set `n_jobs=-1` 
-to use all available CPU cores for cross-cluster multiprocessing, or specify a 
-positive integer for a fixed number of workers 
+- Parallel distance computation via Numba `prange` and the new `n_jobs` 
+parameter. Set `use_numba=True` with `n_jobs=-1` to use all available CPU 
+cores, providing 2-3x speedups on multi-core machines 
 ([Issue #36](https://github.com/vc1492a/PyNomaly/issues/36)).
-- Numba `prange`-based parallel kernels for thread-level parallelism when 
-`use_numba=True` and `n_jobs > 1`.
 - Optional `scipy` acceleration: uses `scipy.spatial.distance.cdist` for 
 distance computation and `scipy.special.erf` for the error function when 
 scipy is available, with graceful fallback to pure NumPy.
 - Dedicated Numba-specific tests (`test_numba_*`) that automatically run when 
 Numba is installed and are skipped otherwise, covering sequential equivalence, 
-parallel equivalence, and progress bar integration with Numba.
-### Fixed
-- Numba `prange` dispatch now correctly activates for single-cluster data 
-when `n_jobs > 1`, providing 2-3x speedups that were previously blocked by 
-incorrect capping of the parallelism flag to the number of clusters.
+parallel equivalence, single-cluster prange, and progress bar integration.
 ### Changed
 - Replaced the O(n^2) Python nested loop for distance computation with a 
 vectorized NumPy implementation using chunked broadcasting. Progress bar 
@@ -31,12 +25,12 @@ with Numba's parallel mode) to non-generator kernels that support `numba.prange`
 - Vectorized `_standard_distances`, `_prob_distances`, and 
 `_norm_prob_outlier_factor` pipeline methods, replacing Python `for` loops 
 with NumPy array operations.
-- Replaced the previous `num_threads` parameter (from the unreleased 
-`feature/numba_parallel` branch) with `n_jobs`, which provides a unified 
-interface for both process-level parallelism (via `concurrent.futures`) and 
-Numba thread-level parallelism (via `prange`). The `n_jobs` parameter follows 
-the scikit-learn convention (`-1` for all cores, positive integer for a fixed 
-number of workers).
+- The `n_jobs` parameter controls Numba thread-level parallelism via `prange`, 
+following the scikit-learn convention (`-1` for all cores, positive integer 
+for a fixed number of threads). Replaces the previous `num_threads` parameter 
+from the unreleased `feature/numba_parallel` branch. `n_jobs` only takes 
+effect when `use_numba=True`; without Numba a warning is issued and processing 
+falls back to the sequential vectorized path.
 
 ## 0.3.5
 ### Changed
diff --git a/examples/numba_speed_diff.py b/examples/numba_speed_diff.py
index cc36ba7..25260ee 100644
--- a/examples/numba_speed_diff.py
+++ b/examples/numba_speed_diff.py
@@ -29,40 +29,20 @@ def main():
     ).fit().local_outlier_probabilities
     t4 = time.time()
     seconds_numba = t4 - t3
-    print("\nComputation took " + str(seconds_numba) + " seconds with Numba JIT.")
+    print("\nComputation took " + str(seconds_numba) + " seconds with Numba JIT (sequential).")
 
-    # Multi-cluster parallel example
-    np.random.seed(42)
-    cluster_a = np.random.randn(5000, 4) + 0
-    cluster_b = np.random.randn(5000, 4) + 10
-    multi_data = np.vstack([cluster_a, cluster_b])
-    cluster_labels = [0] * 5000 + [1] * 5000
-
-    # Sequential (n_jobs=1) with clusters
+    # Numba JIT (parallel, n_jobs=-1)
     t5 = time.time()
-    scores_seq = loop.LocalOutlierProbability(
-        multi_data,
-        n_neighbors=3,
-        cluster_labels=cluster_labels,
-        n_jobs=1,
-        progress_bar=True
-    ).fit().local_outlier_probabilities
-    t6 = time.time()
-    seconds_seq = t6 - t5
-    print("\nComputation took " + str(seconds_seq) + " seconds with n_jobs=1 (2 clusters).")
-
-    # Parallel (n_jobs=-1) with clusters
-    t7 = time.time()
-    scores_par = loop.LocalOutlierProbability(
-        multi_data,
+    scores_numba_par = loop.LocalOutlierProbability(
+        data,
         n_neighbors=3,
-        cluster_labels=cluster_labels,
+        use_numba=True,
         n_jobs=-1,
         progress_bar=True
     ).fit().local_outlier_probabilities
-    t8 = time.time()
-    seconds_par = t8 - t7
-    print("\nComputation took " + str(seconds_par) + " seconds with n_jobs=-1 (2 clusters).")
+    t6 = time.time()
+    seconds_numba_par = t6 - t5
+    print("\nComputation took " + str(seconds_numba_par) + " seconds with Numba JIT (parallel, n_jobs=-1).")
 
 
 if __name__ == '__main__':
diff --git a/examples/parallel_benchmark.py b/examples/parallel_benchmark.py
index b71a41d..432a165 100644
--- a/examples/parallel_benchmark.py
+++ b/examples/parallel_benchmark.py
@@ -1,11 +1,11 @@
 """
-Thorough benchmark of PyNomaly's parallel processing capabilities.
+Benchmark of PyNomaly's Numba parallel processing (prange).
 
 Tests across multiple dimensions:
   - Data sizes (total observations)
-  - Number of clusters (parallelism granularity)
-  - Execution modes: vectorized, multiprocessing (n_jobs), Numba sequential, Numba parallel
-  - Number of workers (1, 2, 4, all cores)
+  - Number of clusters
+  - Sequential vs parallel Numba (n_jobs=1 vs n_jobs=-1)
+  - Comparison against the default vectorized NumPy path
 
 Run:
     python examples/parallel_benchmark.py
@@ -40,8 +40,7 @@ def generate_clustered_data(n_per_cluster, n_clusters, n_features=N_FEATURES, se
     return data, labels
 
 
-def time_fit(data, cluster_labels, n_neighbors, use_numba, n_jobs, warmup=False):
-    """Time a single fit call. Returns elapsed seconds."""
+def time_fit(data, cluster_labels, n_neighbors, use_numba, n_jobs):
     clf = loop.LocalOutlierProbability(
         data, n_neighbors=n_neighbors,
         cluster_labels=cluster_labels,
@@ -50,12 +49,10 @@ def time_fit(data, cluster_labels, n_neighbors, use_numba, n_jobs, warmup=False)
     )
     t0 = time.perf_counter()
     clf.fit()
-    elapsed = time.perf_counter() - t0
-    return elapsed
+    return time.perf_counter() - t0
 
 
 def run_benchmark(data, labels, n_neighbors, label, configs):
-    """Run a set of configurations and print results."""
     n_total = len(data)
     n_clusters = len(set(labels))
 
@@ -68,16 +65,14 @@ def run_benchmark(data, labels, n_neighbors, label, configs):
     print(f"  {'-'*35} {'-'*8} {'-'*8} {'-'*8}")
 
     baseline_best = None
-    results = []
 
     for mode_name, use_numba, n_jobs in configs:
         if use_numba and not HAS_NUMBA:
             continue
 
-        # Warmup for Numba JIT compilation
         if use_numba:
             small_data, small_labels = generate_clustered_data(50, n_clusters, seed=99)
-            time_fit(small_data, small_labels, min(n_neighbors, 10), use_numba, n_jobs, warmup=True)
+            time_fit(small_data, small_labels, min(n_neighbors, 10), use_numba, n_jobs)
 
         times = []
         for _ in range(N_REPEATS):
@@ -91,96 +86,63 @@ def run_benchmark(data, labels, n_neighbors, label, configs):
             baseline_best = best
 
         speedup = baseline_best / best if best > 0 else float('inf')
-        results.append((mode_name, best, mean, speedup))
         print(f"  {mode_name:<35} {best:>7.3f}s {mean:>7.3f}s {speedup:>7.2f}x")
 
-    return results
-
 
 def main():
     cpu_count = os.cpu_count() or 1
-    print(f"PyNomaly Parallel Processing Benchmark")
+    print(f"PyNomaly Parallel Benchmark (Numba prange)")
     print(f"CPU cores: {cpu_count}")
     print(f"Numba available: {HAS_NUMBA}")
     print(f"Repeats per config: {N_REPEATS}")
 
-    # Scenario 1: Small data, many clusters — tests multiprocessing dispatch overhead
+    if not HAS_NUMBA:
+        print("\nNumba is not installed. Install it with: pip install numba")
+        print("Only vectorized (baseline) results will be shown.\n")
+
+    configs = [
+        ("Vectorized NumPy (baseline)",     False, 1),
+        ("Numba sequential (n_jobs=1)",      True,  1),
+        ("Numba parallel (n_jobs=-1)",       True,  -1),
+    ]
+
     data, labels = generate_clustered_data(500, 8)
     run_benchmark(data, labels, N_NEIGHBORS,
         "Scenario 1: Small clusters (500 pts x 8 clusters = 4,000 total)",
-        [
-            ("Vectorized (baseline)",       False, 1),
-            ("Multiprocessing n_jobs=2",     False, 2),
-            ("Multiprocessing n_jobs=4",     False, 4),
-            ("Multiprocessing n_jobs=-1",    False, -1),
-            ("Numba sequential",             True,  1),
-            ("Numba parallel (prange)",      True,  -1),
-        ]
-    )
+        configs)
 
-    # Scenario 2: Medium data, moderate clusters — practical use case
     data, labels = generate_clustered_data(2000, 4)
     run_benchmark(data, labels, N_NEIGHBORS,
         "Scenario 2: Medium clusters (2,000 pts x 4 clusters = 8,000 total)",
-        [
-            ("Vectorized (baseline)",       False, 1),
-            ("Multiprocessing n_jobs=2",     False, 2),
-            ("Multiprocessing n_jobs=4",     False, 4),
-            ("Multiprocessing n_jobs=-1",    False, -1),
-            ("Numba sequential",             True,  1),
-            ("Numba parallel (prange)",      True,  -1),
-        ]
-    )
+        configs)
 
-    # Scenario 3: Large data, many clusters — where parallelism should help most
     data, labels = generate_clustered_data(2000, 8)
     run_benchmark(data, labels, N_NEIGHBORS,
-        "Scenario 3: Large + many clusters (2,000 pts x 8 clusters = 16,000 total)",
-        [
-            ("Vectorized (baseline)",       False, 1),
-            ("Multiprocessing n_jobs=2",     False, 2),
-            ("Multiprocessing n_jobs=4",     False, 4),
-            ("Multiprocessing n_jobs=8",     False, 8),
-            ("Multiprocessing n_jobs=-1",    False, -1),
-            ("Numba sequential",             True,  1),
-            ("Numba parallel (prange)",      True,  -1),
-        ]
-    )
+        "Scenario 3: Large (2,000 pts x 8 clusters = 16,000 total)",
+        configs)
 
-    # Scenario 4: Scaling test — fixed cluster size, increasing cluster count
+    # Single cluster scaling — shows prange benefit within a cluster
     print(f"\n{'='*70}")
-    print(f"  Scenario 4: Scaling test — 1,000 pts/cluster, increasing clusters")
+    print(f"  Scenario 4: Single cluster scaling (Numba prange within one cluster)")
     print(f"{'='*70}")
-    print(f"  {'Clusters':>8} {'Total':>8} {'Seq':>8} {'Par':>8} {'Speedup':>8}")
-    print(f"  {'-'*8} {'-'*8} {'-'*8} {'-'*8} {'-'*8}")
-
-    for n_clusters in [2, 4, 8, 16]:
-        data, labels = generate_clustered_data(1000, n_clusters)
-
-        seq_times = []
-        par_times = []
-        for _ in range(N_REPEATS):
-            seq_times.append(time_fit(data, labels, N_NEIGHBORS, False, 1))
-            par_times.append(time_fit(data, labels, N_NEIGHBORS, False, -1))
-
-        seq_best = min(seq_times)
-        par_best = min(par_times)
-        speedup = seq_best / par_best if par_best > 0 else float('inf')
-        total = n_clusters * 1000
-
-        print(f"  {n_clusters:>8} {total:>7,} {seq_best:>7.3f}s {par_best:>7.3f}s {speedup:>7.2f}x")
-
-    # Scenario 5: Single large cluster — no parallelism benefit expected
-    data, labels = generate_clustered_data(10000, 1)
-    run_benchmark(data, labels, N_NEIGHBORS,
-        "Scenario 5: Single cluster (10,000 pts) — no multi-cluster parallelism",
-        [
-            ("Vectorized (baseline)",       False, 1),
-            ("Multiprocessing n_jobs=-1",    False, -1),
-            ("Numba sequential",             True,  1),
-            ("Numba parallel (prange)",      True,  -1),
-        ]
-    )
+    print(f"  {'Points':>8} {'Vectorized':>10} {'Numba Seq':>10} {'Numba Par':>10} {'Speedup':>8}")
+    print(f"  {'-'*8} {'-'*10} {'-'*10} {'-'*10} {'-'*8}")
+
+    for n in [5000, 10000, 20000]:
+        data, labels = generate_clustered_data(n, 1)
+        vec_t = min(time_fit(data, labels, N_NEIGHBORS, False, 1) for _ in range(N_REPEATS))
+
+        if HAS_NUMBA:
+            small_data, small_labels = generate_clustered_data(50, 1, seed=99)
+            time_fit(small_data, small_labels, 10, True, 1)
+            time_fit(small_data, small_labels, 10, True, -1)
+
+            nb_s = min(time_fit(data, labels, N_NEIGHBORS, True, 1) for _ in range(N_REPEATS))
+            nb_p = min(time_fit(data, labels, N_NEIGHBORS, True, -1) for _ in range(N_REPEATS))
+            speedup = vec_t / nb_p
+            print(f"  {n:>7,} {vec_t:>9.2f}s {nb_s:>9.2f}s {nb_p:>9.2f}s {speedup:>7.2f}x")
+        else:
+            print(f"  {n:>7,} {vec_t:>9.2f}s {'N/A':>10} {'N/A':>10} {'N/A':>8}")
 
     print(f"\n{'='*70}")
     print("  Benchmark complete.")
diff --git a/tests/test_loop.py b/tests/test_loop.py
index 5a7106a..939a686 100644
--- a/tests/test_loop.py
+++ b/tests/test_loop.py
@@ -842,22 +842,19 @@ def test_vectorized_1d_data() -> None:
     assert scores[-1] > 0
 
 
-def test_parallel_with_progress_bar(X_n140_outliers) -> None:
+def test_n_jobs_without_numba_warns(X_n120) -> None:
     """
-    Tests that parallel processing works with the progress bar enabled,
-    covering the progress_bar branch inside _distances_parallel.
+    Tests that n_jobs > 1 without use_numba=True warns the user and
+    falls back to sequential vectorized processing.
     """
-    a = [0] * 120
-    b = [1] * 20
-    cluster_labels = a + b
+    with pytest.warns(UserWarning, match="n_jobs > 1 requires use_numba=True"):
+        clf = loop.LocalOutlierProbability(
+            X_n120, n_neighbors=10, n_jobs=2
+        )
+        scores = clf.fit().local_outlier_probabilities
 
-    clf = loop.LocalOutlierProbability(
-        X_n140_outliers, n_neighbors=10, cluster_labels=cluster_labels,
-        n_jobs=2, progress_bar=True
-    )
-    scores = clf.fit().local_outlier_probabilities
     assert scores is not None
-    assert len(scores) == len(X_n140_outliers)
+    assert len(scores) == len(X_n120)
 
 
 def test_n_jobs_negative_two() -> None:
@@ -887,42 +884,6 @@ def test_vectorized_progress_bar_single_cluster(X_n120) -> None:
     assert len(scores) == len(X_n120)
 
 
-def test_n_jobs_equivalence(X_n140_outliers) -> None:
-    """
-    Tests that n_jobs > 1 produces equivalent results to n_jobs=1 when
-    using cluster labels (multiple clusters processed in parallel).
-    """
-    a = [0] * 120
-    b = [1] * 20
-    cluster_labels = a + b
-
-    clf_seq = loop.LocalOutlierProbability(
-        X_n140_outliers, n_neighbors=10, cluster_labels=cluster_labels, n_jobs=1
-    )
-    scores_seq = clf_seq.fit().local_outlier_probabilities
-
-    clf_par = loop.LocalOutlierProbability(
-        X_n140_outliers, n_neighbors=10, cluster_labels=cluster_labels, n_jobs=2
-    )
-    scores_par = clf_par.fit().local_outlier_probabilities
-
-    assert_array_almost_equal(scores_seq, scores_par, decimal=10)
-
-
-def test_n_jobs_single_cluster(X_n120) -> None:
-    """
-    Tests that n_jobs=-1 works correctly with a single cluster (falls back
-    to sequential since there is only one cluster to process).
-    """
-    clf1 = loop.LocalOutlierProbability(X_n120, n_neighbors=10, n_jobs=1)
-    scores1 = clf1.fit().local_outlier_probabilities
-
-    clf2 = loop.LocalOutlierProbability(X_n120, n_neighbors=10, n_jobs=-1)
-    scores2 = clf2.fit().local_outlier_probabilities
-
-    assert_array_almost_equal(scores1, scores2, decimal=10)
-
-
 def test_n_jobs_invalid() -> None:
     """
     Tests that invalid n_jobs values produce a warning and default to 1.