"""Functions operating on arrays."""
# Any commits made to this module between 2015-05-01 and 2017-03-01
# by Gerrit Holl are developed for the EC project “Fidelity and
# Uncertainty in Climate Data Records from Earth Observations (FIDUCEO)”.
# Grant agreement: 638822
#
# All those contributions are dual-licensed under the MIT license for use
# in typhon, and the GNU General Public License version 3.
import numpy as np
import scipy.stats
[docs]def localmin(arr):
"""Find local minima for 1-D array
Given a 1-dimensional numpy.ndarray, return the locations of any local
minimum as a boolean array. The first and last item are always
considered False.
Arguments:
localmin (numpy.ndarray): 1-D ndarray for which to find local
minima. Should have a numeric dtype.
Returns:
numpy.ndarray with dtype `bool`. True for any element that is
strictly smaller than both neighbouring elements. First and last
element are always False.
"""
localmin = np.hstack(
(False, (arr[1:-1] < arr[0:-2]) & (arr[1:-1] < arr[2:]), False)
)
return localmin
[docs]def limit_ndarray(M, limits):
"""Select elements from structured ndarray based on value ranges
This function filters a structured ndarray based on ranges defined for
zero or more fields. For each field f with limits (lo, hi), it will
select only those elements where lo<=X[f]<hi.
>>> X = array([(2, 3), (4, 5), (8, 2), (5, 1)],
dtype=[("A", "i4"), ("B", "i4")])
>>> print(limit_ndarray(X, {"A": (2, 5)}))
[(2, 3) (4, 5)]
>>> X = array([([2, 3], 3), ([4, 6], 5), ([8, 3], 2), ([5, 3], 1)],
dtype=[("A", "i4", 2), ("B", "i4")])
>>> print(limit_ndarray(X, {"A": (2, 5, "all")}))
[([2, 3], 3)]
Arguments:
M (numpy.ndarray): 1-D structured ndarray
limits (dict): Dictionary with limits. Keys must correspond to
fields in M. If this is a scalar field
(`M.dtype[field].shape==()`), values are tuples (lo, hi).
If this is a multidimensional field, values are tuples (lo,
hi, mode), where mode must be either `all` or `any`.
Values in the range [lo, hi) are retained, applying all or any
when needed.
Returns:
ndarray subset of M. This is a view, not a copy.
"""
selection = np.ones(shape=M.shape, dtype="?")
for (field, val) in limits.items():
ndim = len(M.dtype[field].shape)
if ndim == 0:
(lo, hi) = val
selection = selection & (M[field] >= lo) & (M[field] < hi)
else:
(lo, hi, mode) = val
lelo = M[field] >= lo
sthi = M[field] < hi
while lelo.ndim > 1:
lelo = getattr(lelo, mode)(-1)
sthi = getattr(sthi, mode)(-1)
selection = selection & lelo & sthi
return M[selection]
[docs]def parity(v):
"""Vectorised parity-checking.
For any ndarray with an nd.integer dtype, return an equally shaped
array with the bit parity for each element.
Arguments:
v (numpy.ndarray): Array of integer dtype
Returns:
ndarray with uint8 dtype with the parity for each value in v
"""
v = v.copy() # don't ruin original
parity = np.zeros(dtype=">u1", shape=v.shape)
while v.any():
parity[v != 0] += 1
v &= v - 1
return parity
[docs]def mad_outliers(arr, cutoff=10, mad0="raise"):
"""Mask out mad outliers
Mask out any values that are more than N times the median absolute
devitation from the median.
Although I (Gerrit Holl) came up with this myself, it's also
documented at:
http://eurekastatistics.com/using-the-median-absolute-deviation-to-find-outliers/
except that I rolled by own approach for "what if mad==0".
Note: If all values except one are constant, it is not possible to
determine whether the remaining one is an outlier or “reasonably
close” to the rest, without additional hints. In this case, some
outliers may go unnoticed.
Arguments:
arr (numpy.ndarray): n-D array with numeric dtype
cutoff (int): Maximum tolerable normalised fractional distance
mad0 (str): What to do if mad=0. Can be 'raise', 'ignore', or
'perc'. In case of 'perc', will search for the lowest
percentile at which the percentile absolute deviation is
nonzero, increase the cutoff by the fractional approach toward
percentile 100, and use that percentile instead. So if the
first non-zero is at percentile 75%, it will use the
75th-percntile-absolute-deviation and increase the cutoff by
a factor (100 - 50)/(100 - 75).
Returns:
ndarray with bool dtype, True for outliers
"""
if arr.ptp() == 0:
return np.zeros(shape=arr.shape, dtype="?")
ad = abs(arr - np.ma.median(arr))
mad = np.ma.median(ad)
if mad == 0:
if mad0 == "raise":
raise ValueError("Cannot filter outliers, MAD=0")
elif mad0 == "perc":
# try other percentiles
perc = np.r_[np.arange(50, 99, 1), np.linspace(99, 100, 100)]
pad = scipy.stats.scoreatpercentile(ad, perc)
if (pad == 0).all(): # all constant…?
raise ValueError("These data are weird!")
p_i = pad.nonzero()[0][0]
cutoff *= (100 - 50) / (100 - perc[p_i])
return (ad / pad[p_i]) > cutoff
elif mad is np.ma.masked:
# all are masked already…
return np.ones(shape=ad.shape, dtype="?")
else:
return (ad / mad) > cutoff
[docs]def argclosest(array, value, retvalue=False):
"""Returns the index of the closest value in array.
Parameters:
array (ndarray): Input array.
value (float): Value to compare to.
retvalue (bool): If True, return the index and the closest value.
Returns:
int, float:
Index of closest value, Closest value (if ``retvalue`` is True)
"""
idx = np.abs(np.asarray(array) - value).argmin()
return (idx, array[idx]) if retvalue else idx