"""
Bivariate statistical tools.
* elbow detection in an (x,y) plot
* peaks: finding peaks, quantifying their height, width, center, area, left & right boundaries
* area under curve
"""
import numpy as np
from ..regression.methods import fit_robust_lm
[docs]
def find_elbow_point(x: np.ndarray, y: np.ndarray, max_iter: int = 41) -> int | float:
"""
Find the elbow point when plotting numeric entries in `x` vs numeric values in list `y`.
Using a robust linear fit, sorts the samples in X (independent variable)
and takes the first 5 samples from the left, and the last 5 from the right,
then fits two linear regressions and computes the intersection of the two
fitted lines. The window size is then grown over `max_iter` (default 41)
evenly spaced steps, via `numpy.linspace`, up to roughly half the data,
accumulating one intersection point per step.
The elbow is taken as the data point whose (x, y) location is closest to
the median of the accumulated intersection points; the median location is
where the intersections should stabilise.
Both axes are first scaled to span 0 to 1, so that "closest" does not add up
distances in two different units. The elbow is therefore the same whatever
units `x` and `y` are in: time in seconds or in hours gives the same index.
Will probably not work well on few data points. If so, try fitting a spline
to the raw data and then repeat with the interpolated data.
Parameters
----------
x : np.ndarray
Independent-variable values. Must have the same length as `y`. NaN
entries are dropped pairwise with `y` before fitting; more than 10
finite pairs must remain.
y : np.ndarray
Dependent-variable values. Must have the same length as `x`.
max_iter : int, optional
Number of evenly spaced window sizes to sweep when accumulating
intersection points. Default is 41.
Returns
-------
int or float
The 0-based index of the elbow point in the (sorted) vectors.
Returns ``-1`` if every value in ``x`` or ``y`` is missing.
Returns ``np.nan`` when the intersection sweep produced only
NaNs (for example, when every candidate line pair was
near-parallel), so no consensus intersection could be formed.
Raises
------
ValueError
If `x` and `y` have different lengths, or if 10 or fewer finite
pairs remain after NaN entries are dropped.
"""
start = 5
# assert divmod(max_iter, 2)[1] # must be odd number; to ensure we calculate the median later
def calculate_line_length(x1: float, y1: float, x2: np.ndarray, y2: np.ndarray) -> float | np.ndarray:
r"""Return the length of the line between 2 points (x1, y1) and (x2, y2).
Defined as :math:`\\sqrt{(x2 - x1)^2 + (y2 - y1)^2}`.
"""
return np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
# Ensure it is a Numpy array: pandas objects and lists are correctly handled.
x = np.array(x.copy())
y = np.array(y.copy())
if len(x) != len(y):
raise ValueError(f"x and y must have the same length; got len(x)={len(x)}, len(y)={len(y)}.")
# Stop if everything is missing:
if np.isnan(np.nanmedian(x)) or np.isnan(np.nanmedian(y)):
return -1
# Eliminate missing values in x and y simultaneously.
x, y = x[~(np.isnan(x) | np.isnan(y))], y[~(np.isnan(x) | np.isnan(y))]
if len(x) <= 10:
raise ValueError("Requires more than 10 values in the vectors (not including missing data).")
idx_sort = x.argsort()
x = _scale_to_unit_range(x[idx_sort])
y = _scale_to_unit_range(y[idx_sort])
N = x.size
int_x_list: list[float] = []
int_y_list: list[float] = []
# On fewer than about 2 * max_iter points the floored window sizes repeat, so each size is
# fitted once. A repeat still adds its intersection again, which keeps the median unchanged.
intersections: dict[int, tuple] = {}
for i in np.floor(np.linspace(0, int(N / 2 - start) + 1, max_iter)):
idx = int(start + i)
if idx not in intersections:
lft = fit_robust_lm(x[0:idx], y[0:idx])
rgt = fit_robust_lm(x[-idx:], y[-idx:])
intersections[idx] = find_line_intersection(lft[1], lft[0], rgt[1], rgt[0])
int_x, int_y = intersections[idx]
int_x_list.append(int_x)
int_y_list.append(int_y)
# The elbow is the raw data point closest to the consensus intersection
# point. Taking the median of every accumulated intersection - in both x
# and y - keeps the estimate robust to the occasional spurious window fit
# (e.g. near-parallel lines yielding a NaN or far-off intersection),
# instead of relying on a single, arbitrarily selected intersection.
mid_x = np.nanmedian(int_x_list)
mid_y = np.nanmedian(int_y_list)
if np.isnan(mid_x) or np.isnan(mid_y):
return np.nan
return int(np.argmin(calculate_line_length(mid_x, mid_y, x, y)))
def _scale_to_unit_range(values: np.ndarray) -> np.ndarray:
"""Shift and scale ``values`` to span 0 to 1; a constant vector is only shifted, to 0."""
span = float(np.ptp(values))
return (values - values.min()) / (span if span > 0 else 1.0)
[docs]
def find_line_intersection(m1: float, b1: float, m2: float, b2: float) -> tuple:
"""
Find the intersection point of two lines.
From Stackoverflow:
stackoverflow.com/questions/20677795/how-do-i-compute-the-intersection-point-of-two-lines
Returns a tuple: (x, y) where the two lines intersect, given slopes `m1` and
`m2`, and intercepts `b1` and `b2`.
"""
# These lines are essentially parallel!
if np.abs(m1 - m2) < np.sqrt(np.finfo(float).eps):
return np.nan, np.nan
# y = mx + b
# Derivation: Set both lines equal to find the intersection point in the x direction
# m1 * x + b1 = m2 * x + b2
# m1 * x - m2 * x = b2 - b1
# x * (m1 - m2) = b2 - b1
# Solving the above equation for x:
x = (b2 - b1) / (m1 - m2)
# Now solve it for y: use either line, because they are equal here:
y = m1 * x + b1
return x, y
# ENG-23 (#305): explicit ``__all__`` so the thin re-exporter ``methods.py``
# can do ``from ._elbow_peak import *`` without triggering CodeQL's
# py/polluting-import warning.
__all__ = [
"find_elbow_point",
"find_line_intersection",
"fit_robust_lm",
]
