diff --git a/PyNomaly/loop.py b/PyNomaly/loop.py index 48f499b..31a32b1 100644 --- a/PyNomaly/loop.py +++ b/PyNomaly/loop.py @@ -1,17 +1,83 @@ 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 __author__ = "Valentino Constantinou" -__version__ = "0.3.5" +__version__ = "0.4.0" __license__ = "Apache License, Version 2.0" @@ -74,6 +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: 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: """ """ @@ -315,7 +384,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 +418,7 @@ def new_f(*args, **kwds): (int, np.integer), list, bool, + (int, np.integer), bool, ) def __init__( @@ -359,6 +430,7 @@ def __init__( n_neighbors=10, cluster_labels=None, use_numba=False, + n_jobs=1, progress_bar=False, ) -> None: self.data = data @@ -368,6 +440,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 +456,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 +524,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 +539,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 +645,77 @@ 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(self, progress_bar: bool = False) -> None: """ Provides the distances between each observation and it's closest @@ -576,27 +728,42 @@ 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 + + if self.use_numba: + self._distances_numba( + clusters, distances, indexes, 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 + ) self.distance_matrix = distances self.neighbor_matrix = indexes @@ -626,18 +793,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 +811,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..8b46a45 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) @@ -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,26 +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. -### Utilizing Numba and Progress Bars +### Speeding Things Up with Numba -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`: +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, progress_bar=True).fit() +m = loop.LocalOutlierProbability(data, extent=2, n_neighbors=20, use_numba=True).fit() 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 avaialble in `examples/numba_speed_diff.py`. +To go further, set `n_jobs=-1` to enable Numba's thread-level parallelism +(`prange`), which distributes work across all available CPU cores: -You may also choose to print progress bars _with our without_ the use of numba -by passing `progress_bar=True` to the `LocalOutlierProbability()` method as above. +```python +from PyNomaly import loop +m = loop.LocalOutlierProbability( + data, extent=2, n_neighbors=20, + use_numba=True, n_jobs=-1 +).fit() +scores = m.local_outlier_probabilities +print(scores) +``` + +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. + +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 `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. + +### Progress Bars + +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, use_numba=True, n_jobs=-1, progress_bar=True).fit() +``` + +Progress bars are supported in both sequential and Numba execution modes. ### Choosing Parameters @@ -431,7 +466,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 ee92cdd..9b41565 100644 --- a/changelog.md +++ b/changelog.md @@ -4,6 +4,34 @@ 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 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)). +- 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, 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 +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. +- 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 - Refactored the `Validate` class by dissolving it and moving validation methods diff --git a/examples/numba_speed_diff.py b/examples/numba_speed_diff.py index 2d4f622..25260ee 100644 --- a/examples/numba_speed_diff.py +++ b/examples/numba_speed_diff.py @@ -2,30 +2,48 @@ 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 +def main(): + # generate a large set of data + data = np.ones(shape=(10000, 4)) -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.") + # 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.") -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.") + # 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 (sequential).") + + # Numba JIT (parallel, n_jobs=-1) + t5 = time.time() + scores_numba_par = loop.LocalOutlierProbability( + data, + n_neighbors=3, + use_numba=True, + n_jobs=-1, + progress_bar=True + ).fit().local_outlier_probabilities + 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__': + main() diff --git a/examples/parallel_benchmark.py b/examples/parallel_benchmark.py new file mode 100644 index 0000000..432a165 --- /dev/null +++ b/examples/parallel_benchmark.py @@ -0,0 +1,153 @@ +""" +Benchmark of PyNomaly's Numba parallel processing (prange). + +Tests across multiple dimensions: + - Data sizes (total observations) + - 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 +""" + +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): + 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() + return time.perf_counter() - t0 + + +def run_benchmark(data, labels, n_neighbors, label, configs): + 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 + + for mode_name, use_numba, n_jobs in configs: + if use_numba and not HAS_NUMBA: + continue + + 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) + + 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') + print(f" {mode_name:<35} {best:>7.3f}s {mean:>7.3f}s {speedup:>7.2f}x") + + +def main(): + cpu_count = os.cpu_count() or 1 + 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}") + + 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)", + configs) + + 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)", + configs) + + data, labels = generate_clustered_data(2000, 8) + run_benchmark(data, labels, N_NEIGHBORS, + "Scenario 3: Large (2,000 pts x 8 clusters = 16,000 total)", + configs) + + # Single cluster scaling — shows prange benefit within a cluster + print(f"\n{'='*70}") + print(f" Scenario 4: Single cluster scaling (Numba prange within one cluster)") + print(f"{'='*70}") + 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.") + print(f"{'='*70}") + + +if __name__ == '__main__': + main() 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=[], diff --git a/tests/test_loop.py b/tests/test_loop.py index b9f4aaf..939a686 100644 --- a/tests/test_loop.py +++ b/tests/test_loop.py @@ -827,3 +827,155 @@ 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_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_n_jobs_without_numba_warns(X_n120) -> None: + """ + Tests that n_jobs > 1 without use_numba=True warns the user and + falls back to sequential vectorized processing. + """ + 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 + + assert scores is not None + assert len(scores) == len(X_n120) + + +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_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 + + +# --- 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)