pyod/models/rod.py
# -*- coding: utf-8 -*-
"""Rotation-based Outlier Detector (ROD)
"""
# Author: Yahya Almardeny <almardeny@gmail.com>
# License: BSD 2 clause
from __future__ import division
from __future__ import print_function
import multiprocessing
from itertools import combinations as com
from multiprocessing import Pool
import numba
import numpy as np
from sklearn.preprocessing import MinMaxScaler, RobustScaler
from sklearn.utils import check_array
from .base import BaseDetector
@numba.njit
def mad(costs, median=None):
"""Apply the robust median absolute deviation (MAD)
to measure the inconsistency/variability of the
rotation costs.
Parameters
----------
costs : list of rotation costs
median: float (default=None), MAD median
Returns
-------
z : float
the modified z scores
"""
costs_ = np.reshape(costs, (-1, 1))
median = np.nanmedian(costs_) if median is None else median
diff = np.abs(costs_ - median)
return np.ravel(0.6745 * diff / np.median(diff)), median
def angle(v1, v2):
"""find the angle between two 3D vectors.
Parameters
----------
v1 : list, first vector
v2 : list, second vector
Returns
-------
angle : float, the angle
"""
return np.arccos(np.dot(v1, v2) /
(np.linalg.norm(v1) * np.linalg.norm(v2)))
def geometric_median(x, eps=1e-5):
"""
Find the multivariate geometric L1-median by applying
Vardi and Zhang algorithm.
Parameters
----------
x : array-like, the data points
eps: float (default=1e-5), a threshold to indicate when to stop
Returns
-------
gm : array, Geometric L1-median
"""
points = np.unique(x, axis=0)
gm_ = np.mean(points, 0) # initialize geometric median
while True:
D = euclidean(points, gm_, c=True)
non_zeros = (D != 0)[:, 0]
Dinv = 1 / D[non_zeros]
Dinvs = np.sum(Dinv)
W = Dinv / Dinvs
T = np.sum(W * points[non_zeros], 0)
num_zeros = len(points) - np.sum(non_zeros)
if num_zeros == 0:
gm1 = T
elif num_zeros == len(points):
return gm_
else:
R = (T - gm_) * Dinvs
r = np.linalg.norm(R)
r_inv = 0 if r == 0 else num_zeros / r
gm1 = max(0, 1 - r_inv) * T + min(1, r_inv) * gm_
if euclidean(gm_, gm1) < eps:
return gm1
gm_ = gm1
def scale_angles(gammas, scaler1=None, scaler2=None):
"""
Scale all angles in which angles <= 90
degree will be scaled within [0 - 54.7] and
angles > 90 will be scaled within [90 - 126]
Parameters
----------
gammas : list, angles
scaler1: obj (default=None), MinMaxScaler of Angles group 1
scaler2: obj (default=None), MinMaxScaler of Angles group 2
Returns
-------
scaled angles, scaler1, scaler2
"""
first, second = [], []
first_ind, second_ind = [], []
q1 = np.pi / 2.
for i, g in enumerate(gammas):
if g <= q1:
first.append(g)
first_ind.append(i)
else:
second.append(g)
second_ind.append(i)
if scaler1 is None: # this indicates the `fit()`
min_f, max_f = 0.001, 0.955
scaler1 = MinMaxScaler(feature_range=(min_f, max_f))
# min_f and max_f are required to be fit by scaler for consistency between train and test sets
scaler1.fit(np.array(first + [min_f, max_f]).reshape(-1, 1))
first = scaler1.transform(np.array(first).reshape(-1, 1)).reshape(-1) if first else []
else:
first = scaler1.transform(np.array(first).reshape(-1, 1)).reshape(-1) if first else []
if scaler2 is None: # this indicates the `fit()`
min_s, max_s = q1 + 0.001, 2.186
scaler2 = MinMaxScaler(feature_range=(min_s, max_s))
# min_s and max_s are required to be fit by scaler for consistency between train and test sets
scaler2.fit(np.array(second + [min_s, max_s]).reshape(-1, 1))
second = scaler2.transform(np.array(second).reshape(-1, 1)).reshape(-1) if second else []
else:
second = scaler2.transform(np.array(second).reshape(-1, 1)).reshape(-1) if second else []
# restore original order
return np.concatenate([first, second])[np.argsort(first_ind + second_ind)], scaler1, scaler2
def euclidean(v1, v2, c=False):
"""
Find the euclidean distance between two vectors
or between a vector and a collection of vectors.
Parameters
----------
v1 : list, first 3D vector or collection of vectors
v2 : list, second 3D vector
c : bool (default=False), if True, it means the v1 is a list of vectors.
Returns
-------
list of list of euclidean distances if c==True.
Otherwise float: the euclidean distance
"""
if c:
res = []
for _v in v1:
res.append([np.sqrt((_v[0] - v2[0]) ** 2 +
(_v[1] - v2[1]) ** 2 +
(_v[2] - v2[2]) ** 2)])
return np.array(res, copy=False)
return np.sqrt((v1[0] - v2[0]) ** 2 +
(v1[1] - v2[1]) ** 2 +
(v1[2] - v2[2]) ** 2)
def rod_3D(x, gm=None, median=None, scaler1=None, scaler2=None):
"""
Find ROD scores for 3D Data.
note that gm, scaler1 and scaler2 will be returned "as they are"
and without being changed if the model has been fit already
Parameters
----------
x : array-like, 3D data points.
gm: list (default=None), the geometric median
median: float (default=None), MAD median
scaler1: obj (default=None), MinMaxScaler of Angles group 1
scaler2: obj (default=None), MinMaxScaler of Angles group 2
Returns
-------
decision_scores, gm, scaler1, scaler2
"""
# find the geometric median if it is not already fit
gm = geometric_median(x) if gm is None else gm
# find its norm and center data around it
norm_ = np.linalg.norm(gm)
_x = x - gm
# calculate the scaled angles between the geometric median and each data point vector
v_norm = np.linalg.norm(_x, axis=1)
gammas, scaler1, scaler2 = scale_angles(np.arccos(np.clip(np.dot(_x, gm) / (v_norm * norm_), -1, 1)),
scaler1=scaler1, scaler2=scaler2)
# apply the ROD main equation to find the rotation costs
costs = np.power(v_norm, 3) * np.cos(gammas) * np.square(np.sin(gammas))
# apply MAD to calculate the decision scores
decision_scores, median = mad(costs, median=median)
return decision_scores, list(gm), median, scaler1, scaler2
@numba.njit
def sigmoid(x):
"""
Implementation of Sigmoid function
Parameters
----------
x : array-like, decision scores
Returns
-------
array-like, x after applying sigmoid
"""
return 1 / (1 + np.exp(-x))
def process_sub(subspace, gm, median, scaler1, scaler2):
"""
Apply ROD on a 3D subSpace then process it with sigmoid
to compare apples to apples
Parameters
----------
subspace : array-like, 3D subspace of the data
gm: list, the geometric median
median: float, MAD median
scaler1: obj, MinMaxScaler of Angles group 1
scaler2: obj, MinMaxScaler of Angles group 2
Returns
-------
ROD decision scores with sigmoid applied, gm, scaler1, scaler2
"""
mad_subspace, gm, median, scaler1, scaler2 = rod_3D(subspace, gm=gm,
median=median,
scaler1=scaler1,
scaler2=scaler2)
return sigmoid(np.nan_to_num(np.array(mad_subspace))), gm, median, scaler1, scaler2
def rod_nD(X, parallel, gm=None, median=None, data_scaler=None, angles_scalers1=None, angles_scalers2=None):
"""
Find ROD overall scores when Data is higher than 3D:
# scale dataset using Robust Scaler
# decompose the full space into a combinations of 3D subspaces,
# Apply ROD on each combination,
# squish scores per subspace, so we compare apples to apples,
# calculate average of ROD scores of all subspaces per observation.
Note that if gm, data_scaler, angles_scalers1, angles_scalers2 are None,
that means it is a `fit()` process and they will be calculated and returned
to the class to be saved for future prediction. Otherwise, if they are not None,
then it is a prediction process.
Parameters
----------
X : array-like, data points
parallel: bool, True runs the algorithm in parallel
gm: list (default=None), the geometric median
median: list (default=None), MAD medians
data_scaler: obj (default=None), RobustScaler of data
angles_scalers1: list (default=None), MinMaxScalers of Angles group 1
angles_scalers2: list (default=None), MinMaxScalers of Angles group 2
Returns
-------
ROD decision scores, gm, median, data_scaler, angles_scalers1, angles_scalers2
"""
if data_scaler is None: # for fitting
data_scaler = RobustScaler()
X = data_scaler.fit_transform(X)
else: # for prediction
X = data_scaler.transform(X)
dim = X.shape[1]
all_subspaces = [X[:, _com] for _com in com(range(dim), 3)]
all_gms = [None] * len(all_subspaces) if gm is None else gm
all_meds = [None] * len(all_subspaces) if median is None else median
all_angles_scalers1 = [None] * len(all_subspaces) if angles_scalers1 is None else angles_scalers1
all_angles_scalers2 = [None] * len(all_subspaces) if angles_scalers2 is None else angles_scalers2
if parallel:
p = Pool(multiprocessing.cpu_count())
args = [[a, b, c, d, e] for a, b, c, d, e in zip(all_subspaces, all_gms, all_meds,
all_angles_scalers1, all_angles_scalers2)]
results = p.starmap(process_sub, args)
subspaces_scores, gm, median, angles_scalers1, angles_scalers2 = [], [], [], [], []
for res in results:
subspaces_scores.append(list(res[0]))
gm.append(res[1])
median.append(res[2])
angles_scalers1.append(res[3])
angles_scalers2.append(res[4])
scores = np.average(np.array(subspaces_scores).T, axis=1).reshape(-1)
p.close()
p.join()
return scores, gm, median, data_scaler, angles_scalers1, angles_scalers2
subspaces_scores, gm, median, angles_scalers1, angles_scalers2 = [], [], [], [], []
for subspace, _gm, med, ang_s1, ang_s2 in zip(all_subspaces, all_gms, all_meds, all_angles_scalers1,
all_angles_scalers2):
scores_, gm_, med_, ang_s1_, ang_s2_ = process_sub(subspace=subspace, gm=_gm, median=med,
scaler1=ang_s1, scaler2=ang_s2)
subspaces_scores.append(scores_)
gm.append(gm_)
median.append(med_)
angles_scalers1.append(ang_s1_)
angles_scalers2.append(ang_s2_)
scores = np.average(np.array(subspaces_scores).T, axis=1).reshape(-1)
return scores, gm, median, data_scaler, angles_scalers1, angles_scalers2
class ROD(BaseDetector):
"""Rotation-based Outlier Detection (ROD), is a robust and parameter-free
algorithm that requires no statistical distribution assumptions and
works intuitively in three-dimensional space, where the 3D-vectors,
representing the data points, are rotated about the geometric median
two times counterclockwise using Rodrigues rotation formula.
The results of the rotation are parallelepipeds where their volumes are
mathematically analyzed as cost functions and used to calculate the
Median Absolute Deviations to obtain the outlying score.
For high dimensions > 3, the overall score is calculated by taking the
average of the overall 3D-subspaces scores, that were resulted from
decomposing the original data space.
See :cite:`almardeny2020novel` for details.
Parameters
----------
contamination : float in (0., 0.5), optional (default=0.1)
The amount of contamination of the data set, i.e.
the proportion of outliers in the data set. Used when fitting to
define the threshold on the decision function.
parallel_execution: bool, optional (default=False).
If set to True, the algorithm will run in parallel,
for a better execution time. It is recommended to set
this parameter to True ONLY for high dimensional data > 10,
and if a proper hardware is available.
Attributes
----------
decision_scores_ : numpy array of shape (n_samples,)
The outlier scores of the training data.
The higher, the more abnormal. Outliers tend to have higher
scores. This value is available once the detector is
fitted.
threshold_ : float
The threshold is based on ``contamination``. It is the
``n_samples * contamination`` most abnormal samples in
``decision_scores_``. The threshold is calculated for generating
binary outlier labels.
labels_ : int, either 0 or 1
The binary labels of the training data. 0 stands for inliers
and 1 for outliers/anomalies. It is generated by applying
``threshold_`` on ``decision_scores_``.
"""
def __init__(self, contamination=0.1, parallel_execution=False):
super(ROD, self).__init__(contamination=contamination)
if not isinstance(parallel_execution, bool):
raise TypeError("parallel_execution should be bool. "
"Got {}".format(type(parallel_execution)))
self.parallel_execution = parallel_execution
def fit(self, X, y=None):
"""Fit detector. y is ignored in unsupervised methods.
Parameters
----------
X : numpy array of shape (n_samples, n_features)
The input samples.
y : Ignored
Not used, present for API consistency by convention.
Returns
-------
self : object
Fitted estimator.
"""
X = check_array(X)
self._set_n_classes(y)
# reset learning parameters after each fit
self.gm_ = None # geometric median(s)
self.median_ = None # MAD median(s)
self.data_scaler_ = None # data scaler (in case of d>3)
self.angles_scaler1_ = None # scaler(s) of Angles Group 1
self.angles_scaler2_ = None # scaler(s) of Angles Group 2
self.decision_scores_ = self.decision_function(X)
self._process_decision_scores()
return self
def decision_function(self, X):
"""Predict raw anomaly score of X using the fitted detector.
The anomaly score of an input sample is computed based on different
detector algorithms. For consistency, outliers are assigned with
larger anomaly scores.
Parameters
----------
X : numpy array of shape (n_samples, n_features)
The training input samples. Sparse matrices are accepted only
if they are supported by the base estimator.
Returns
-------
anomaly_scores : numpy array of shape (n_samples,)
The anomaly score of the input samples.
"""
X = check_array(X)
if X.shape[1] < 3:
X = np.hstack((X, np.zeros(shape=(X.shape[0], 3 - X.shape[1]))))
if X.shape[1] == 3:
scores, self.gm_, self.median_, self.angles_scaler1_, \
self.angles_scaler2_ = rod_3D(x=X, gm=self.gm_,
median=self.median_,
scaler1=self.angles_scaler1_,
scaler2=self.angles_scaler2_)
return scores
scores, self.gm_, self.median_, self.data_scaler_, \
self.angles_scaler1_, self.angles_scaler2_ = rod_nD(X=X,
parallel=self.parallel_execution,
gm=self.gm_,
median=self.median_,
data_scaler=self.data_scaler_,
angles_scalers1=self.angles_scaler1_,
angles_scalers2=self.angles_scaler2_)
return scores