kneeliverse.kneedle

  1# coding: utf-8
  2
  3'''
  4The following module provides knee detection method
  5based on Kneedle algorithm.
  6'''
  7
  8__author__ = 'Mário Antunes'
  9__version__ = '1.0'
 10__email__ = 'mario.antunes@ua.pt'
 11__status__ = 'Development'
 12__license__ = 'MIT'
 13__copyright__ = '''
 14Copyright (c) 2021-2023 Stony Brook University
 15Copyright (c) 2021-2023 The Research Foundation of SUNY
 16'''
 17
 18import math
 19import enum
 20import logging
 21import numpy as np
 22import uts.ema as ema
 23import uts.peak_detection as pd
 24import kneeliverse.multi_knee as mk
 25import kneeliverse.linear_fit as lf
 26
 27
 28logger = logging.getLogger(__name__)
 29
 30
 31class Direction(enum.Enum):
 32    """
 33    Enum data type that represents the direction of a concavity.
 34    """
 35    Increasing = 'increasing'
 36    Decreasing = 'decreasing'
 37
 38    def __str__(self):
 39        return self.value
 40
 41
 42class Concavity(enum.Enum):
 43    """
 44    Enum data type that represents the rotation of a concavity.
 45    """
 46    Counterclockwise = 'counter-clockwise'
 47    Clockwise = 'clockwise'
 48
 49    def __str__(self):
 50        return self.value
 51
 52
 53class PeakDetection(enum.Enum):
 54    """
 55    Enum data type that identifies the peak selection algorithm.
 56    """
 57    Kneedle = 'Kneedle'
 58    ZScore = 'ZScore'
 59    Significant = 'Significant'
 60    All = 'All'
 61
 62    def __str__(self):
 63        return self.value
 64
 65
 66def differences(points: np.ndarray, cd: Direction, cc: Concavity) -> np.ndarray:
 67    """
 68    Computes the differences from the y axis.
 69
 70    These differences represent a rotation within the original algorithm.
 71
 72    Args:
 73        points (np.ndarray): numpy array with the points (x, y)
 74        cd (Direction): direction of the concavity
 75        cc (Concavity): rotation of the concavity
 76
 77    Returns:
 78        np.ndarray: the points array with the differences
 79    """
 80
 81    rv = np.empty(points.shape)
 82
 83    if cd is Direction.Decreasing and cc is Concavity.Clockwise:
 84        for i in range(0, len(points)):
 85            rv[i][0] = points[i][0]
 86            rv[i][1] = points[i][0] + points[i][1]  # x + y
 87    elif cd is Direction.Decreasing and cc is Concavity.Counterclockwise:
 88        for i in range(0, len(points)):
 89            rv[i][0] = points[i][0]
 90            rv[i][1] = 1.0 - (points[i][0] + points[i][1])  # 1.0 - (x + y)
 91    elif cd is Direction.Increasing and cc is Concavity.Clockwise:
 92        for i in range(0, len(points)):
 93            rv[i][0] = points[i][0]
 94            rv[i][1] = points[i][1] - points[i][0]  # y - x
 95    else:
 96        for i in range(0, len(points)):
 97            rv[i][0] = points[i][0]
 98            rv[i][1] = math.fabs(points[i][1] - points[i][0])  # abs(y - x)
 99
100    return rv
101
102
103def _knee(points: np.ndarray, t: float, cd: Direction, cc: Concavity) -> int:
104    """
105    Returns the index of the knee point based on the Kneedle method.
106
107    Args:
108        points (np.ndarray): numpy array with the points (x, y)
109        t (float): tau of the side window used to smooth the curve
110        cd (Direction): direction of the concavity
111        cc (Concavity): rotation of the concavity
112
113    Returns:
114        int: the index of the knee point
115    """
116
117    Ds = ema.linear(points, t)
118    pmin = Ds.min(axis=0)
119    pmax = Ds.max(axis=0)
120    diff = pmax - pmin
121    diff[diff == 0] = 1.0
122
123    Dn = (Ds - pmin)/diff 
124
125    Dd = differences(Dn, cd, cc)
126
127    peaks = pd.all_peaks(Dd)
128    
129    idx = pd.highest_peak(Dd, peaks)
130    if idx == -1:
131        return None
132    else:
133        return idx
134
135
136def _knees(points: np.ndarray, t: float, cd: Direction, cc: Concavity, sensitivity:float=1.0, p:PeakDetection=PeakDetection.Kneedle, debug:bool=False) -> np.ndarray:
137    """
138    Returns the index of the knees point based on the Kneedle method.
139
140    This implementation uses an heuristic to automatically define
141    the direction and rotation of the concavity.
142
143    Furthermore, it support three different methods to select the 
144    relevant knees:
145    1. Kneedle    : classical algorithm
146    2. Significant: significant knee peak detection
147    3. ZScore     : significant knee peak detection based on zscore
148
149    Args:
150        points (np.ndarray): numpy array with the points (x, y)
151        t (float): tau of the sliding window used to smooth the curve
152        cd (Direction): direction of the concavity
153        cc (Concavity): rotation of the concavity
154        sensitivity (float): controls the sensitivity of the peak detection (default 1.0)
155        p (PeakDetection): selects the peak detection method (default PeakDetection.Kneedle)
156        debug (bool): debug flag; when True the algorithm returns more information
157
158    Returns:
159        np.ndarray: the indexes of the knee points
160    """
161    Ds = ema.linear(points, t)
162
163    pmin = Ds.min(axis=0)
164    pmax = Ds.max(axis=0)
165    Dn = (Ds - pmin)/(pmax - pmin)
166
167    Dd = differences(Dn, cd, cc)
168
169    knees = []
170
171    peaks_idx = pd.all_peaks(Dd)
172
173    if p is PeakDetection.Kneedle:
174        knees = pd.kneedle_peak_detection(Dd, peaks_idx, sensitivity)
175    elif p is PeakDetection.Significant:
176        knees = pd.significant_peaks(Dd, peaks_idx, sensitivity)
177    elif p is PeakDetection.ZScore:
178        knees = pd.significant_zscore_peaks(Dd, peaks_idx, sensitivity)
179    else:
180        knees = peaks_idx
181
182    if debug is True:
183        return {'knees': knees, 'dd': Dd, 'peaks': pd.all_peaks(Dd)}
184    else:
185        return knees
186
187
188def knees(points: np.ndarray,  t: float = 1.0, sensitivity: float = 1.0, p: PeakDetection = PeakDetection.Kneedle) -> np.ndarray:
189    """
190    Returns the index of the knees point based on the Kneedle method.
191
192    This implementation uses an heuristic to automatically define
193    the direction and rotation of the concavity.
194
195    Furthermore, it support three different methods to select the 
196    relevant knees:
197    1. Kneedle    : classical algorithm
198    2. Significant: significant knee peak detection
199    3. ZScore     : significant knee peak detection based on zscore
200
201    Args:
202        points (np.ndarray): numpy array with the points (x, y)
203        t (float): tau of the side window used to smooth the curve
204        sensitivity (float): controls the sensitivity of the peak detection
205        p (PeakDetection): selects the peak detection method
206
207    Returns:
208        np.ndarray: the indexes of the knee points
209    """
210    _, m = lf.linear_fit_points(points)
211
212    if m > 0.0:
213        cd = Direction.Increasing
214    else:
215        cd = Direction.Decreasing
216
217    knees_1= _knees(points, t, cd, Concavity.Counterclockwise, sensitivity, p)
218    knees_2 = _knees(points, t, cd, Concavity.Clockwise, sensitivity, p)
219
220    knees_idx = np.concatenate((knees_1, knees_2))
221    # np.concatenate generates float array when one is empty (see https://github.com/numpy/numpy/issues/8878)
222    knees_idx = knees_idx.astype(int)
223    knees_idx = np.unique(knees_idx)
224    knees_idx.sort()
225
226    return knees_idx
227
228
229def knee(points: np.ndarray, t: float = 1.0) -> int:
230    """
231    Returns the index of the knee point based on the Kneedle method.
232
233    This implementation uses an heuristic to automatically define
234    the direction and rotation of the concavity.
235
236    Args:
237        points (np.ndarray): numpy array with the points (x, y)
238        t (float): tau of the side window used to smooth the curve
239
240    Returns:
241        int: the index of the knee point
242    """
243    b, m = lf.linear_fit_points(points)
244
245    if m > 0.0:
246        cd = Direction.Increasing
247    else:
248        cd = Direction.Decreasing
249
250    y = points[:, 1]
251    yhat = np.empty(len(points))
252    for i in range(0, len(points)):
253        yhat[i] = points[i][0]*m+b
254
255    vote = np.sum(y - yhat)
256
257    if cd is Direction.Increasing and vote > 0:
258        cc = Concavity.Clockwise
259    elif cd is Direction.Increasing and vote <= 0:
260        cc = Concavity.Counterclockwise
261    elif cd is Direction.Decreasing and vote > 0:
262        cc = Concavity.Clockwise
263    else:
264        cc = Concavity.Counterclockwise
265
266    return _knee(points, t, cd, cc)
267
268
269def multi_knee(points, t1=0.01, t2=3):
270    """
271    Recursive knee point detection based on Kneedle.
272
273    It returns the knee points on the curve.
274
275    Args:
276        points (np.ndarray): numpy array with the points (x, y)
277        t1 (float): coefficient of determination threshold (default 0.01)
278        t2 (int): number of points threshold (default 3)
279
280    Returns:
281        np.ndarray: knee points on the curve
282
283    """
284    knees = mk.multi_knee(knee, points, t1, t2)
285    return knees