diff --git a/pybaselines/polynomial.py b/pybaselines/polynomial.py index c6b7bc7..912d961 100644 --- a/pybaselines/polynomial.py +++ b/pybaselines/polynomial.py @@ -364,11 +364,11 @@ def penalized_poly(self, data, poly_order=2, tol=1e-3, max_iter=250, weights=Non return_coef=return_coef ) - @_Algorithm._handle_io(sort_keys=('weights', 'coef'), require_unique=True) + @_Algorithm._handle_io(sort_keys=('weights', 'coef'), require_unique=True, mask_support=2) def loess(self, data, fraction=0.2, total_points=None, poly_order=1, scale=3.0, tol=1e-3, max_iter=10, symmetric_weights=False, use_threshold=False, num_std=1, use_original=False, weights=None, return_coef=False, - conserve_memory=True, delta=None): + conserve_memory='deprecated', delta=None, sigma_func=None): """ Locally estimated scatterplot smoothing (LOESS). @@ -421,15 +421,13 @@ def loess(self, data, fraction=0.2, total_points=None, poly_order=1, scale=3.0, If True, will convert the polynomial coefficients for the fit baseline to a form that fits the input x_data and return them in the params dictionary. Default is False, since the conversion takes time. - conserve_memory : bool, optional - If False, will cache the distance-weighted kernels for each value - in `x_data` on the first iteration and reuse them on subsequent iterations to - save time. The shape of the array of kernels is (len(`x_data`), `total_points`). - If True (default), will recalculate the kernels each iteration, which uses very - little memory, but is slower. Can usually set to False unless `x_data` and`total_points` - are quite large and the function causes memory issues when caching the kernels. If - numba is installed, there is no significant time difference since the calculations are - sped up. + conserve_memory + + .. deprecated:: 1.3 + `conserve_memory` is deprecated and no longer has any effect. If + calculation time is too slow, increase `delta` instead. Will be removed + in version 1.5. + delta : float, optional If `delta` is > 0, will skip all but the last x-value in the range `x_last + delta`, where `x_last` is the last x-value to be fit using weighted least squares, and instead @@ -438,6 +436,18 @@ def loess(self, data, fraction=0.2, total_points=None, poly_order=1, scale=3.0, statsmodels [6]_ and Cleveland's original Fortran lowess implementation [7]_). Fits all x-values if `delta` is <= 0. Default is None, which sets `delta` to `0.01 * (max(x_data) - min(x_data))`. + sigma_func : Callable, optional + The function that calculates the estimated experimental noise by which the residuals + are scaled before weighting each iteration. Must have the call signature:: + + sigma_func(residual: numpy.ndarray) -> float + + where the input `residual` array is ``data - baseline`` for that iteration, and the + output is a float. If None (default), the function is the standardized median + absolute value, ``median(abs(residual)) / 0.6745`` [2]_. Other options, as suggested + by the REBS (robust extraction of baseline signal) method [8]_, are to use the + standard deviation of residual values less than 0 or less than the mode if values + are known to follow some distribution. Returns ------- @@ -464,8 +474,19 @@ def loess(self, data, fraction=0.2, total_points=None, poly_order=1, scale=3.0, ------ ValueError Raised if the number of points per window for the fitting is less than - `poly_order` + 1 or greater than the total number of points, or if the - values in `self.x` are not strictly increasing. + ``poly_order + 2`` or greater than the total number of points, or if the + x-values are not strictly increasing. Also raised if a window had too few + non-zero weights during fitting when ``Baseline.mask`` is not None and + the object was initialized with `strict_mask=True`. + TypeError + Raised if the input `sigma_func` does not return a float. + + Warns + ----- + ParameterWarning + Emitted if `poly_order` is greater than 2, which can introduce numerical + issues with small window sizes. Also emitted if a window had too few + non-zero weights during fitting when ``Baseline.mask`` is None. Notes ----- @@ -474,7 +495,7 @@ def loess(self, data, fraction=0.2, total_points=None, poly_order=1, scale=3.0, fit data based on the residuals, as proposed by [3]_, similar to the modpoly and imodpoly techniques. - In baseline literature, this procedure is sometimes called "rbe", meaning + In baseline literature, this procedure is often called "rbe", meaning "robust baseline estimate". References @@ -496,23 +517,35 @@ def loess(self, data, fraction=0.2, total_points=None, poly_order=1, scale=3.0, 1363-1367. .. [6] https://github.com/statsmodels/statsmodels. .. [7] https://www.netlib.org/go (lowess.f is the file). + .. [8] Ruckstuhl, A.F., et al. Robust extraction of baseline signal of atmospheric trace + species using local regression. Atmospheric Measurement Techniques, 2012, + 5(11), 2613-2624. """ if total_points is None: total_points = ceil(fraction * self._size) - if total_points < poly_order + 1: - raise ValueError('total points must be greater than polynomial order + 1') + # cutoff is poly_order + 2 rather than poly_order + 1 for other polynomials since + # furthest point in each window has a weight of 0 + if total_points < poly_order + 2: + raise ValueError('total points must be greater than polynomial order + 2') elif total_points > self._size: raise ValueError(( 'points per window is higher than total number of points; lower either ' '"fraction" or "total_points"' )) - elif poly_order > 2: + if poly_order > 2: warnings.warn( ('polynomial orders greater than 2 can have numerical issues;' ' consider using a polynomial order of 1 or 2 instead'), ParameterWarning, stacklevel=2 ) + if conserve_memory != 'deprecated': + warnings.warn( + ('conserve_memory is deprecated and no longer has any effect'), + DeprecationWarning, stacklevel=2 + ) + if sigma_func is None: + sigma_func = _median_absolute_value y, weight_array = self._setup_polynomial(data, weights, poly_order, calc_vander=True) if use_original: @@ -546,19 +579,19 @@ def loess(self, data, fraction=0.2, total_points=None, poly_order=1, scale=3.0, # do max_iter + 1 since a max_iter of 0 would return y as baseline otherwise for i in range(max_iter + 1): baseline_old = baseline - if conserve_memory: - baseline = _loess_low_memory( - x, y, sqrt_w, coefs, vandermonde, self._size, windows, fits - ) - elif i == 0: - kernels, baseline = _loess_first_loop( - x, y, sqrt_w, coefs, vandermonde, total_points, self._size, windows, fits - ) - else: - baseline = _loess_nonfirst_loops( - y, sqrt_w, coefs, vandermonde, kernels, windows, self._size, fits - ) - + baseline, zero_wt_window = _loess_low_memory( + x, y, sqrt_w, coefs, vandermonde, self._size, windows, fits + ) + if zero_wt_window: + if self.mask is not None and self._strict_mask: + # hard to know if zero-weighted window was from masking or too small a window, + # so err on the side of caution + raise ValueError('A window had too few non-zero weights to fit') + else: + warnings.warn( + ('A window had too few non-zero weights to fit; likely need to ' + 'increase "fraction" or "total_points"'), ParameterWarning, stacklevel=2 + ) _fill_skips(x, baseline, skips) calc_difference = relative_difference(baseline_old, baseline) @@ -572,11 +605,18 @@ def loess(self, data, fraction=0.2, total_points=None, poly_order=1, scale=3.0, y0 if use_original else y, baseline + num_std * np.std(residual[pos_wt_mask]) ) else: - # TODO median_absolute_value can be 0 if more than half of residuals are - # 0 (perfect fit); can that ever really happen? if so, should prevent dividing by 0 - sqrt_w = _tukey_square( - residual / _median_absolute_value(residual), scale, symmetric_weights - ) + noise_sigma = sigma_func(residual) + if not isinstance(noise_sigma, float): + raise TypeError('"sigma_func" must return a float') + elif noise_sigma < np.finfo(float).eps: + # break rather than setting noise_sigma to eps since any further reweighting + # will still be unstable; see statsmodels issue#2108 + warnings.warn( + 'calculated noise scale is near 0; terminating early', ParameterWarning, + stacklevel=3 + ) + break + sqrt_w = _tukey_square(residual, scale * noise_sigma, symmetric_weights) params = {'weights': sqrt_w**2, 'tol_history': tol_history[:i + 1]} if return_coef: @@ -1181,12 +1221,12 @@ def _tukey_square(residual, scale=3, symmetric=False): """ if symmetric: inner = residual / scale - weights = np.maximum(0, 1 - inner * inner) + weights = np.maximum(0., 1. - inner**2) else: weights = np.ones_like(residual) mask = residual > 0 inner = residual[mask] / scale - weights[mask] = np.maximum(0, 1 - inner * inner) + weights[mask] = np.maximum(0., 1. - inner**2) return weights @@ -1328,6 +1368,8 @@ def _loess_low_memory(x, y, weights, coefs, vander, num_x, windows, fits): ------- baseline : numpy.ndarray, shape (N,) The calculated baseline. + zero_wt_window : bool + Whether any window had too few non-zero values during fitting. Notes ----- @@ -1337,6 +1379,8 @@ def _loess_low_memory(x, y, weights, coefs, vander, num_x, windows, fits): baseline = np.empty(num_x) y_fit = y * weights vander_fit = vander.T * weights + min_nonzero = vander.shape[1] + zero_wt_window = False for idx in range(fits.shape[0]): i = fits[idx] window = windows[idx] @@ -1349,137 +1393,20 @@ def _loess_low_memory(x, y, weights, coefs, vander, num_x, windows, fits): difference = 1 - difference kernel = np.sqrt(difference * difference * difference) - coef = _loess_solver( - kernel * vander_fit[:, left:right], kernel * y_fit[left:right] - ) - baseline[i] = vander[i].dot(coef) - coefs[i] = coef - - return baseline - - -# adapted from (https://gist.github.com/agramfort/850437); see license above -@jit(nopython=True, cache=True) -def _loess_first_loop(x, y, weights, coefs, vander, total_points, num_x, windows, fits): - """ - The initial fit for loess that also caches the window values for each x-value. - - Parameters - ---------- - x : numpy.ndarray, shape (N,) - The x-values of the measured data, with N data points. - y : numpy.ndarray, shape (N,) - The y-values of the measured data, with N points. - weights : numpy.ndarray, shape (N,) - The array of weights. - coefs : numpy.ndarray, shape (N, ``poly_order + 1``) - The array of polynomial coefficients (with polynomial order poly_order), - for each value in `x`. - vander : numpy.ndarray, shape (N, ``poly_order + 1``) - The Vandermonde matrix for the `x` array. - total_points : int - The number of points to include when fitting each x-value. - num_x : int - The number of data points in `x`, also known as N. - windows : numpy.ndarray, shape (F, 2) - An array of left and right indices that define the fitting window for each fit - x-value. The length is F, which is the total number of fit points. If `fit_dx` - is <= 0, F is equal to N, the total number of x-values. - fits : numpy.ndarray, shape (F,) - The array of indices indicating which x-values to fit. - - Returns - ------- - kernels : numpy.ndarray, shape (N, total_points) - The array containing the distance-weighted kernel for each x-value. - baseline : numpy.ndarray, shape (N,) - The calculated baseline. - - Notes - ----- - The coefficient array, `coefs`, is modified inplace. - - """ - kernels = np.empty((num_x, total_points)) - baseline = np.empty(num_x) - y_fit = y * weights - vander_fit = vander.T * weights - for idx in range(fits.shape[0]): - i = fits[idx] - window = windows[idx] - left = window[0] - right = window[1] - - difference = np.abs(x[left:right] - x[i]) - difference = difference / max(difference[0], difference[-1]) - difference = difference * difference * difference - difference = 1 - difference - kernel = np.sqrt(difference * difference * difference) - - kernels[i] = kernel - coef = _loess_solver( - kernel * vander_fit[:, left:right], kernel * y_fit[left:right] - ) - baseline[i] = vander[i].dot(coef) - coefs[i] = coef - - return kernels, baseline - - -@jit(nopython=True, cache=True) -def _loess_nonfirst_loops(y, weights, coefs, vander, kernels, windows, num_x, fits): - """ - The loess fit to use after the first loop that uses the cached window values. - - Parameters - ---------- - y : numpy.ndarray, shape (N,) - The y-values of the measured data, with N points. - weights : numpy.ndarray, shape (N,) - The array of weights. - coefs : numpy.ndarray, shape (N, ``poly_order + 1``) - The array of polynomial coefficients (with polynomial order poly_order), - for each value in `x`. - vander : numpy.ndarray, shape (N, ``poly_order + 1``) - The Vandermonde matrix for the `x` array. - kernels : numpy.ndarray, shape (N, total_points) - The array containing the distance-weighted kernel for each x-value. Each - kernel has a length of total_points. - windows : numpy.ndarray, shape (F, 2) - An array of left and right indices that define the fitting window for each fit - x-value. The length is F, which is the total number of fit points. If `fit_dx` - is <= 0, F is equal to N, the total number of x-values. - num_x : int - The total number of values, N. - fits : numpy.ndarray, shape (F,) - The array of indices indicating which x-values to fit. - - Returns - ------- - baseline : numpy.ndarray, shape (N,) - The calculated baseline. - - Notes - ----- - The coefficient array, `coefs`, is modified inplace. - - """ - baseline = np.empty(num_x) - y_fit = y * weights - vander_fit = vander.T * weights - for idx in range(fits.shape[0]): - i = fits[idx] - window = windows[idx] - left = window[0] - right = window[1] - kernel = kernels[i] - coef = _loess_solver( - kernel * vander_fit[:, left:right], kernel * y_fit[left:right] - ) - baseline[i] = vander[i].dot(coef) - coefs[i] = coef + # statsmodels uses 1e-12 guard against 0 weights (statsmodels issue #7700), + # but pybaselines uses sqrt(weights), so guard against sqrt(1e-12) + if (kernel * weights[left:right] > 1e-6).sum() < min_nonzero: + zero_wt_window = True + baseline[i] = y[i] + coefs[i] = np.nan + else: + coef = _loess_solver( + kernel * vander_fit[:, left:right], kernel * y_fit[left:right] + ) + baseline[i] = vander[i].dot(coef) + coefs[i] = coef - return baseline + return baseline, zero_wt_window @jit(nopython=True, cache=True) @@ -1602,7 +1529,8 @@ def _determine_fits(x, num_x, total_points, delta): @_polynomial_wrapper def loess(data, x_data=None, fraction=0.2, total_points=None, poly_order=1, scale=3.0, tol=1e-3, max_iter=10, symmetric_weights=False, use_threshold=False, num_std=1, - use_original=False, weights=None, return_coef=False, conserve_memory=True, delta=None): + use_original=False, weights=None, return_coef=False, conserve_memory='deprecated', + delta=None, sigma_func=None): """ Locally estimated scatterplot smoothing (LOESS). @@ -1675,6 +1603,18 @@ def loess(data, x_data=None, fraction=0.2, total_points=None, poly_order=1, scal statsmodels [14]_ and Cleveland's original Fortran lowess implementation [15]_). Fits all x-values if `delta` is <= 0. Default is None, which sets `delta` to `0.01 * (max(x_data) - min(x_data))`. + sigma_func : callable, optional + The function that calculates the estimated experimental noise by which the residuals + are scaled before weighting each iteration. Must have the call signature:: + + sigma_func(residual: numpy.ndarray) -> float + + where the input `residual` array is ``data - baseline`` for that iteration, and the + output is a float. If None (default), the function is the standardized median + absolute value, ``median(abs(residual)) / 0.6745`` [10]_. Other options, as suggested + by the REBS (robust extraction of baseline signal) method [20]_, are to use the + standard deviation of residual values less than 0 or less than the mode if values + are known to follow some distribution. Returns ------- @@ -1732,6 +1672,9 @@ def loess(data, x_data=None, fraction=0.2, total_points=None, poly_order=1, scal 1363-1367. .. [14] https://github.com/statsmodels/statsmodels. .. [15] https://www.netlib.org/go (lowess.f is the file). + .. [20] Ruckstuhl, A.F., et al. Robust extraction of baseline signal of atmospheric trace + species using local regression. Atmospheric Measurement Techniques, 2012, + 5(11), 2613-2624. """ diff --git a/tests/data/lowess_zero_weights_iter1.csv b/tests/data/lowess_zero_weights_iter1.csv new file mode 100644 index 0000000..aec7276 --- /dev/null +++ b/tests/data/lowess_zero_weights_iter1.csv @@ -0,0 +1,240 @@ +2.953578502381795090e+01 +2.977249488319329629e+01 +3.000832015933116281e+01 +3.024228863322001004e+01 +3.047427316477404702e+01 +3.070789998897616968e+01 +3.090987519675835316e+01 +3.108066701034287860e+01 +3.120768611426809969e+01 +3.126626931507050600e+01 +3.127134083831706945e+01 +3.125264276656585238e+01 +3.120309709351434435e+01 +3.112053142017954244e+01 +3.101408040917840836e+01 +3.089204569903130349e+01 +3.075147353254864413e+01 +3.059130452522651211e+01 +3.037901188956348264e+01 +3.009714168364473252e+01 +2.972868946430979520e+01 +2.922052138743913474e+01 +2.873767976984627381e+01 +2.823758890625844487e+01 +2.763918369507365469e+01 +2.691735833950625079e+01 +2.612400169777359693e+01 +2.524842025851127048e+01 +2.433710441681530767e+01 +2.340402944259044205e+01 +2.244650223366981834e+01 +2.147155656440631688e+01 +2.048942146412870358e+01 +1.953132259832426598e+01 +1.859336881684354736e+01 +1.764532944138871429e+01 +1.671425000005206130e+01 +1.715412999999999855e+01 +1.702539000000000158e+01 +1.723645000000000138e+01 +1.769518000000000058e+01 +1.847265000000000157e+01 +1.949915999999999983e+01 +2.087391999999999825e+01 +2.247628999999999877e+01 +2.381081000005909232e+01 +2.629326490731614285e+01 +2.870220960640705243e+01 +3.112292949874907677e+01 +3.356936456910040079e+01 +3.602488357238080852e+01 +3.848459738422875631e+01 +4.096045475950738535e+01 +4.348111340749738929e+01 +4.602679741264559254e+01 +4.838567999997947311e+01 +4.950215000000000032e+01 +5.128018000000000143e+01 +5.267683000000000249e+01 +5.387601000000000084e+01 +5.498996000000000350e+01 +5.589578999999999809e+01 +5.645094999999999885e+01 +5.688656000000000290e+01 +5.715155000000000030e+01 +5.716919000000000040e+01 +5.704115000000000180e+01 +5.760801999982504640e+01 +5.682428000013054259e+01 +5.606303300595880756e+01 +5.532078025365849783e+01 +5.447862945390608047e+01 +5.361051539546360090e+01 +5.272475199384587796e+01 +5.179868661718538192e+01 +5.083917009175981150e+01 +4.982694806680797939e+01 +4.878002474467118788e+01 +4.771582169963846098e+01 +4.663326090424532566e+01 +4.554475258641577540e+01 +4.445931169190671994e+01 +4.337231786354050911e+01 +4.229865951191304418e+01 +4.124887440950502082e+01 +4.022181352592309622e+01 +3.921158199806066591e+01 +3.822033561833716675e+01 +3.725205718834188673e+01 +3.631742859851436833e+01 +3.540548355696436289e+01 +3.450428803578780190e+01 +3.361714876814359343e+01 +3.276811575403154109e+01 +3.194614974383823380e+01 +3.114494959849039191e+01 +3.037215090864775391e+01 +2.964575695261322608e+01 +2.895530646558572840e+01 +2.828993215765453684e+01 +2.764745674543856424e+01 +2.702498954521917796e+01 +2.643215047422939179e+01 +2.587044403118083835e+01 +2.532427434921262943e+01 +2.480056591527436893e+01 +2.431172020661229283e+01 +2.384654111802724330e+01 +2.339846243623420818e+01 +2.296998829611773729e+01 +2.256705668683257926e+01 +2.218968903409366789e+01 +2.182979360955506110e+01 +2.148440336667389872e+01 +2.115530849910749467e+01 +2.084409887327100819e+01 +2.054953248706542013e+01 +2.027090555752983647e+01 +2.000819757976653790e+01 +1.976007985947177659e+01 +1.952546223725597940e+01 +1.930407964641751306e+01 +1.909678298185639633e+01 +1.890191751911923390e+01 +1.871934246268384072e+01 +1.854909561127855255e+01 +1.839208268340921038e+01 +1.824711159921893255e+01 +1.811353929137003860e+01 +1.799210053805350640e+01 +1.787985609939927656e+01 +1.777580151777698347e+01 +1.767820444984456785e+01 +1.758725843490891805e+01 +1.750367752772641694e+01 +1.742370253641675149e+01 +1.734966843724850705e+01 +1.728858277984226532e+01 +1.723776232116927432e+01 +1.719451485349423336e+01 +1.715815302042136992e+01 +1.712942097176374290e+01 +1.710828820072620360e+01 +1.708827868641234460e+01 +1.706882142878700037e+01 +1.705610247796730761e+01 +1.705000853823282014e+01 +1.704836296815279439e+01 +1.705074724486661708e+01 +1.705936892679287098e+01 +1.707381547628461149e+01 +1.709202406843397881e+01 +1.711356702319685041e+01 +1.713896150216336522e+01 +1.716743238897197799e+01 +1.719892119998947067e+01 +1.723445988339812374e+01 +1.727435230705458835e+01 +1.731623579938147017e+01 +1.736169533481457705e+01 +1.741258501218886678e+01 +1.746617085249124912e+01 +1.752005780025160675e+01 +1.757325105267359433e+01 +1.762781834071987319e+01 +1.768285756562687894e+01 +1.773721383612251046e+01 +1.779372287198829383e+01 +1.785461533595937311e+01 +1.792030483617186221e+01 +1.798877758719162756e+01 +1.805753095330695857e+01 +1.812777920009659027e+01 +1.819826251659685923e+01 +1.826761635230401026e+01 +1.833840016915023341e+01 +1.841172777922988857e+01 +1.848807472921430062e+01 +1.856627855365931268e+01 +1.864552593462992292e+01 +1.872575212154766078e+01 +1.880760678114398488e+01 +1.889231629788845268e+01 +1.898030419731082219e+01 +1.907055998143158959e+01 +1.916309835434974218e+01 +1.925813641161091638e+01 +1.935362204398268915e+01 +1.944764070993049287e+01 +1.954152508744462224e+01 +1.963579563532072925e+01 +1.972749270083013684e+01 +1.981545623459203753e+01 +1.990181929115959392e+01 +1.998878889069638021e+01 +2.007828341304221453e+01 +2.017096362230518380e+01 +2.026545526792562057e+01 +2.035965107923766482e+01 +2.045269495683356453e+01 +2.054471722800417055e+01 +2.063398681090252751e+01 +2.072091210442592413e+01 +2.080871381821784993e+01 +2.090019597076333113e+01 +2.099331472677294741e+01 +2.108635623630926403e+01 +2.117887841835197449e+01 +2.127208782007945231e+01 +2.136601250846063493e+01 +2.146013808813087564e+01 +2.155325379786478734e+01 +2.164879087596473184e+01 +2.174740200572712112e+01 +2.184489028561688428e+01 +2.193939342196937048e+01 +2.203383237238606185e+01 +2.212997818179000475e+01 +2.222603722328866382e+01 +2.232062057891663542e+01 +2.241589381052817487e+01 +2.251315763642333323e+01 +2.261085063517487725e+01 +2.270947099929177782e+01 +2.280910545823071089e+01 +2.290750334509312225e+01 +2.300448605766455756e+01 +2.310084200668690357e+01 +2.319400637907992646e+01 +2.328342707107736942e+01 +2.336990132689829380e+01 +2.345754963367423684e+01 +2.354608038667946346e+01 +2.363520033128583364e+01 +2.372359756253728236e+01 +2.380753055201298451e+01 +2.389079108330385992e+01 +2.397383923039965836e+01 +2.405622078098507899e+01 +2.413776219608735119e+01 diff --git a/tests/data/lowess_zero_weights_iter2.csv b/tests/data/lowess_zero_weights_iter2.csv new file mode 100644 index 0000000..60bbc94 --- /dev/null +++ b/tests/data/lowess_zero_weights_iter2.csv @@ -0,0 +1,240 @@ +2.954333540733417962e+01 +2.977345908398207897e+01 +3.000329607571281088e+01 +3.023224414555364703e+01 +3.046083671103644974e+01 +3.068833167069062640e+01 +3.087823602352483121e+01 +3.101923330184969174e+01 +3.116668062745565848e+01 +3.126632007289023463e+01 +3.126695278079656504e+01 +3.124842547819765670e+01 +3.119838275355169444e+01 +3.111514889320455168e+01 +3.100900758701208204e+01 +3.089054251540774843e+01 +3.075621153078220971e+01 +3.060923170981066121e+01 +3.044480627542310103e+01 +3.027347033883485139e+01 +2.985513105896750474e+01 +2.896372303987824637e+01 +2.859337894400377422e+01 +2.831571701010567921e+01 +2.782218443758982573e+01 +2.699476516311882790e+01 +2.614097923652488475e+01 +2.523716742694081816e+01 +2.431619574438242637e+01 +2.338546170973097560e+01 +2.243561170552031570e+01 +2.147408533643909578e+01 +2.050891104068051973e+01 +1.974078176421846820e+01 +1.916077129646758692e+01 +1.857944440100555639e+01 +1.804196596340431213e+01 +1.738148701530593598e+01 +1.717608438996403564e+01 +1.757854748707794457e+01 +1.818165132122755878e+01 +1.894482797014354247e+01 +1.986766849759607112e+01 +2.114496286050675877e+01 +2.296487627129967990e+01 +2.490287016951809562e+01 +2.690810203469436246e+01 +2.900619076277325448e+01 +3.119487244967429262e+01 +3.358042669058319518e+01 +3.604354682523317877e+01 +3.850837893124736411e+01 +4.089739976086291051e+01 +4.308275306062339638e+01 +4.519029824122080186e+01 +4.720663827849227800e+01 +4.914670686971172842e+01 +5.098796804009593586e+01 +5.242594682652680405e+01 +5.359482680966491586e+01 +5.461642693035567930e+01 +5.546035201285100413e+01 +5.610458461964563526e+01 +5.657524502276548617e+01 +5.690167684889248534e+01 +5.709234788917573411e+01 +5.696909408833425914e+01 +5.641856893479644697e+01 +5.597967962547853915e+01 +5.554633229471774314e+01 +5.507004783884892163e+01 +5.444256682438926020e+01 +5.358044374919304431e+01 +5.267724629498445665e+01 +5.171721498493612756e+01 +5.078084731061822765e+01 +4.979581800458463903e+01 +4.874854080866000317e+01 +4.769235927465463476e+01 +4.661790511726701425e+01 +4.553535918949314976e+01 +4.445471460152489840e+01 +4.337135171749193319e+01 +4.230155396006081503e+01 +4.125387941770051725e+01 +4.023004628132930804e+01 +3.922026193615140954e+01 +3.822388520529001710e+01 +3.726044961884582563e+01 +3.634436487185843845e+01 +3.544577370593157895e+01 +3.455037819659293064e+01 +3.365417603143591663e+01 +3.278472566548025924e+01 +3.195426594860956016e+01 +3.114738703852019697e+01 +3.036896582723511173e+01 +2.965075631526649147e+01 +2.897087346019050003e+01 +2.830756564461610125e+01 +2.766452148953775136e+01 +2.703476880467837873e+01 +2.643582349158529965e+01 +2.587321969905113761e+01 +2.532329651122887526e+01 +2.480412363237737594e+01 +2.432904668770021317e+01 +2.386852745877022386e+01 +2.341843870194417931e+01 +2.298069079666377590e+01 +2.256820402765366751e+01 +2.219106355860583335e+01 +2.183293315872107243e+01 +2.148771072425676820e+01 +2.115671976623101003e+01 +2.084456154022127805e+01 +2.055010032634399053e+01 +2.027239747257317859e+01 +2.001133723180966228e+01 +1.976397415794552259e+01 +1.952848253032126280e+01 +1.930543232267530129e+01 +1.909725241058439948e+01 +1.890197165397735901e+01 +1.871945420179826769e+01 +1.854970055462140621e+01 +1.839271580806445527e+01 +1.824742560526746971e+01 +1.811375537241946176e+01 +1.799252718360282444e+01 +1.788067416882206118e+01 +1.777716312306801072e+01 +1.767970310223441643e+01 +1.758837708314481318e+01 +1.750457974457459542e+01 +1.742469235594921884e+01 +1.735142344913385060e+01 +1.729239710989837420e+01 +1.724290698777066666e+01 +1.719933992063044670e+01 +1.716167366018315832e+01 +1.713121775069875952e+01 +1.710896233177097159e+01 +1.708840192881550379e+01 +1.706889148773266385e+01 +1.705632394984988309e+01 +1.705030694005287373e+01 +1.704863828282210747e+01 +1.705096025811450033e+01 +1.705945124928798506e+01 +1.707382369477891260e+01 +1.709203970624782087e+01 +1.711363478556197393e+01 +1.713911573194200244e+01 +1.716766452935855725e+01 +1.719926134764100212e+01 +1.723485573700850537e+01 +1.727468075830168814e+01 +1.731645662879468972e+01 +1.736185313273895048e+01 +1.741274678380737129e+01 +1.746633413233232091e+01 +1.752022394119313731e+01 +1.757348961718096092e+01 +1.762813554936333915e+01 +1.768317898016754697e+01 +1.773759075474497493e+01 +1.779435530440911961e+01 +1.785528097160956662e+01 +1.792072632166382107e+01 +1.798899488731089846e+01 +1.805760253041958308e+01 +1.812774366792217506e+01 +1.819810902722508317e+01 +1.826745854795061774e+01 +1.833830283331045052e+01 +1.841172258105241966e+01 +1.848810590686263922e+01 +1.856629794132614109e+01 +1.864554577489320408e+01 +1.872577325533481840e+01 +1.880762276582374426e+01 +1.889232793384308096e+01 +1.898031531432458152e+01 +1.907058833804358855e+01 +1.916314052050764261e+01 +1.925816561153528994e+01 +1.935362393839323047e+01 +1.944762308049480737e+01 +1.954148823154421777e+01 +1.963573243649991085e+01 +1.972739943705089161e+01 +1.981536243734687730e+01 +1.990175922874347592e+01 +1.998876747275941668e+01 +2.007821937874720675e+01 +2.017084586770025467e+01 +2.026522540750643486e+01 +2.035919573528644833e+01 +2.045218027771252522e+01 +2.054431002261370764e+01 +2.063366408809014274e+01 +2.072065607408875820e+01 +2.080866152555616111e+01 +2.090022468147584789e+01 +2.099330491416094446e+01 +2.108632127832482084e+01 +2.117883467374159423e+01 +2.127203585574512701e+01 +2.136596061816178960e+01 +2.146018103599320526e+01 +2.155355498862900276e+01 +2.164930102010314172e+01 +2.174794412642263808e+01 +2.184531575645304002e+01 +2.193970814881419429e+01 +2.203402679012394572e+01 +2.213000412056726418e+01 +2.222589488162934401e+01 +2.232045590936911239e+01 +2.241584861368352577e+01 +2.251315759770086089e+01 +2.261083323081822982e+01 +2.270944011050998768e+01 +2.280906287664970478e+01 +2.290743744776830937e+01 +2.300440675362142784e+01 +2.310087636881958773e+01 +2.319428444199270345e+01 +2.328387088514272563e+01 +2.337043887932165021e+01 +2.345806543072868777e+01 +2.354642086755476171e+01 +2.363513426474200330e+01 +2.372298214174934472e+01 +2.380670011036356470e+01 +2.388985365371507541e+01 +2.397279028478686769e+01 +2.405510036841111088e+01 +2.413662983738944945e+01 diff --git a/tests/test_polynomial.py b/tests/test_polynomial.py index bae0f64..321ccb2 100644 --- a/tests/test_polynomial.py +++ b/tests/test_polynomial.py @@ -7,10 +7,12 @@ """ from math import ceil +from pathlib import Path import numpy as np from numpy.testing import assert_allclose, assert_array_equal import pytest +from scipy import stats from pybaselines import polynomial from pybaselines.utils import ParameterWarning @@ -345,18 +347,13 @@ class TestLoess(IterativePolynomialTester, RecreationMixin, WeightMaskingMixin): func_name = 'loess' allows_zero_iteration = False requires_unique_x = True - supports_mask = False @pytest.mark.parametrize('use_class', (True, False)) @pytest.mark.parametrize('delta', (0, 0.01)) - @pytest.mark.parametrize('conserve_memory', (True, False)) @pytest.mark.parametrize('use_threshold', (True, False)) - def test_unchanged_data(self, use_class, use_threshold, conserve_memory, delta): + def test_unchanged_data(self, use_class, use_threshold, delta): """Ensures that input data is unchanged by the function.""" - super().test_unchanged_data( - use_class, use_threshold=use_threshold, - conserve_memory=conserve_memory, delta=delta - ) + super().test_unchanged_data(use_class, use_threshold=use_threshold, delta=delta) @pytest.mark.parametrize('use_threshold', (True, False)) @pytest.mark.parametrize('use_original', (True, False)) @@ -372,8 +369,8 @@ def test_wrong_fraction_fails(self, fraction): @pytest.mark.parametrize('poly_order', (0, 1, 2, 3)) def test_too_small_window_fails(self, poly_order): - """Ensures a window smaller than poly_order + 1 raises an exception.""" - for num_points in range(poly_order + 1): + """Ensures a window smaller than poly_order + 2 raises an exception.""" + for num_points in range(poly_order + 2): with pytest.raises(ValueError): self.class_func(self.y, total_points=num_points, poly_order=poly_order) @@ -406,8 +403,7 @@ def test_output_coefs(self, poly_order, delta): assert_allclose(baseline, recreated_poly) - @pytest.mark.parametrize('conserve_memory', (True, False)) - def test_compare_to_statsmodels(self, conserve_memory): + def test_compare_to_statsmodels(self): """ Compares the output of loess to the output of statsmodels.lowess. @@ -451,9 +447,8 @@ def test_compare_to_statsmodels(self, conserve_memory): # test several iterations to ensure weighting is correct for iterations in range(4): output = self.algorithm_base(x, check_finite=False, assume_sorted=True).loess( - y, conserve_memory=conserve_memory, total_points=total_points, - max_iter=iterations, tol=-1, scale=4.0469385011764905, symmetric_weights=True, - delta=0.0 + y, total_points=total_points, max_iter=iterations, tol=-1, + scale=4.0469385011764905, symmetric_weights=True, delta=0.0 ) assert_allclose( @@ -517,11 +512,149 @@ def test_input_weights(self, use_threshold): super().test_input_weights(use_threshold=use_threshold) @pytest.mark.threaded_test - @pytest.mark.parametrize('conserve_memory', (True, False)) - def test_threading(self, conserve_memory): + def test_threading(self): """Tests the different possible computation routes under threading.""" - delta = 0.05 * (self.x.max() - self.x.min()) # use a larger delta to speed up method - super().test_threading(conserve_memory=conserve_memory, delta=delta) + # use a larger delta to speed up method + super().test_threading(delta=0.05 * (self.x.max() - self.x.min())) + + def test_custom_sigma_func(self): + """Ensures input sigma_func modifies the reweighting.""" + baseline, params = self.class_func(self.y) + baseline2, params2 = self.class_func(self.y, sigma_func=lambda vals: np.std(vals[vals < 0])) + + # simple check that different sigma calcs produced different baselines and weights + with pytest.raises(AssertionError): + assert_allclose(baseline2, baseline, rtol=1e-4, atol=1e-3) + with pytest.raises(AssertionError): + assert_allclose(params2['weights'], params['weights'], rtol=1e-1, atol=1e-1) + + def test_incorrect_sigma_func_fails(self): + """Ensures an exception is raised if input sigma_func does not return a float.""" + with pytest.raises(TypeError, match='"sigma_func" must return a float'): + self.class_func(self.y, sigma_func=lambda vals: 'a') + + def test_zero_sigma_exits(self): + """Ensures the method exits early when the calculated noise sigma is ~0. + + Replicates statsmodels issue #2108. + + """ + x = np.arange(20) + y = np.array([0] * 10 + [1] * 10, dtype=float) + fraction = 2 / 3 + total_points = int(len(x) * fraction) + with pytest.warns(ParameterWarning, match='calculated noise scale is near 0'): + output = self.algorithm_base(x).loess( + y, total_points=total_points, max_iter=3, delta=0, + scale=4.0469385011764905, symmetric_weights=True, tol=-1 + )[0] + expected_output = np.array([ + 0, 0, 0, 0, 0, 0, 0, 0.03796574, 0.29511209, 0.44982749, 0.55017251, 0.70488791, + 0.96203426, 1, 1, 1, 1, 1, 1, 1 + ]) + assert_allclose(output, expected_output, rtol=1e-6, atol=1e-9) + + def test_zero_sigma_exits_2(self): + """Ensures the method exits early when the calculated noise sigma is ~0. + + Replicates the first part of statsmodels issue #1798. + + """ + x = np.arange(20) + y = np.arange(20, dtype=float) + fraction = 0.4 + total_points = int(len(x) * fraction) + with pytest.warns(ParameterWarning, match='calculated noise scale is near 0'): + output = self.algorithm_base(x).loess( + y, total_points=total_points, max_iter=3, delta=0, + scale=4.0469385011764905, symmetric_weights=True, tol=-1 + )[0] + expected_output = y # should be a perfect fit + assert_allclose(output, expected_output, rtol=1e-14, atol=1e-13) + + @pytest.mark.parametrize('max_iter', (1, 2)) + def test_zero_weights_fill(self, max_iter): + """Ensures a window with zero weights with fill with y instead of causing numerical issues. + + Dataset is adapted from statsmodels issue #7700. The data files for statsmodels's output + were created using:: + + from statsmodels.nonparametric.smoothers_lowess import lowess + output = lowess(y, x, frac=11 / len(x), it=max_iter, delta=0).T[1] + + with statsmodels version 0.14.6. + + """ + y = np.array([ + 29.60046, 29.70066, 29.99869, 30.18495, + 30.52497, 30.88539, 31.06073, 31.16298, 31.3087, 31.34476, 31.4047, 31.27913, + 31.29533, 31.14104, 31.033, 30.95522, 30.7452, 30.6161, 30.48558, 30.20304, + 29.94876, 29.49816, 28.99673, 28.47641, 27.75036, 26.98692, 26.22662, 25.29733, + 24.45699, 23.47883, 22.421, 21.46149, 20.50521, 19.55747, 18.71905, 17.97059, + 17.4616, 17.15413, 17.02539, 17.23645, 17.69518, 18.47265, 19.49916, 20.87392, + 22.47629, 24.34076, 26.46264, 28.66842, 31.13522, 33.57669, 35.95129, 38.50984, + 40.9788, 43.45954, 45.54811, 47.72132, 49.50215, 51.28018, 52.67683, 53.87601, + 54.98996, 55.89579, 56.45095, 56.88656, 57.15155, 57.16919, 57.04115, 56.87761, + 56.42096, 55.93649, 55.2568, 54.47306, 53.79956, 52.8701, 51.84985, 50.93586, + 49.95632, 48.73087, 47.77627, 46.75819, 45.54977, 44.36957, 43.32188, 42.29313, + 41.24385, 40.14291, 39.15614, 38.17805, 37.27126, 36.13561, 35.32942, 34.35569, + 33.69126, 32.67565, 31.91131, 31.0636, 30.32011, 29.60982, 28.88217, 28.10989, + 27.56996, 27.03619, 26.36284, 25.82758, 25.27555, 24.80477, 24.25029, 23.74979, + 23.31028, 22.95834, 22.56406, 22.13128, 21.81209, 21.42739, 21.12386, 20.8205, + 20.52693, 20.26264, 19.94682, 19.74871, 19.47004, 19.28826, 19.09282, 18.8813, + 18.69543, 18.51512, 18.37025, 18.21213, 18.09597, 18.00692, 17.84771, 17.7365, + 17.70439, 17.54311, 17.50521, 17.42641, 17.32607, 17.29374, 17.17156, 17.14076, + 17.18559, 17.12909, 17.11519, 17.06809, 17.05098, 17.06691, 17.02511, 17.01555, + 17.07787, 17.05032, 17.05407, 17.06751, 17.12841, 17.12312, 17.16593, 17.21924, + 17.19979, 17.25681, 17.31144, 17.36246, 17.43259, 17.43767, 17.5086, 17.58345, + 17.62989, 17.70608, 17.70383, 17.81441, 17.82661, 17.8836, 18.00816, 18.05311, + 18.16044, 18.19468, 18.24426, 18.32978, 18.41256, 18.47817, 18.57559, 18.6523, + 18.71417, 18.79602, 18.89392, 18.96791, 19.0598, 19.17692, 19.25897, 19.33334, + 19.45276, 19.56273, 19.63092, 19.71592, 19.83377, 19.91831, 19.97547, 20.07111, + 20.15791, 20.23325, 20.38081, 20.49393, 20.54687, 20.62749, 20.70332, 20.81285, + 20.87916, 21.01356, 21.07556, 21.19642, 21.26882, 21.35373, 21.45083, 21.55625, + 21.66463, 21.75115, 21.8033, 21.9497, 22.06961, 22.1253, 22.20523, 22.32333, + 22.41526, 22.50364, 22.62715, 22.70702, 22.80392, 22.89037, 23.02072, 23.12152, + 23.18633, 23.29179, 23.39558, 23.4171, 23.56042, 23.59962, 23.76348, 23.7985, + 23.93591, 23.97028, 24.04745, 24.12475 + ]) + x = np.linspace(2160, 2559, len(y)) / 60 + with pytest.warns(ParameterWarning, match='A window had too few non-zero weights'): + output = self.algorithm_base(x).loess( + y, poly_order=1, total_points=11, max_iter=max_iter, delta=0, + scale=4.0469385011764905, symmetric_weights=True, tol=-1 + )[0] + expected_output = np.loadtxt( + Path(__file__).parent.joinpath(f'data/lowess_zero_weights_iter{max_iter}.csv') + ) + assert_allclose(output, expected_output, rtol=1e-11, atol=1e-11) + + def test_zero_weights(self): + """Simpler version of test_zero_weights_fill, using input all-zeros weights. + + Allows testing for which indices should fail. + + """ + with pytest.warns(ParameterWarning, match='A window had too few non-zero weights'): + output, params = self.class_func( + self.y, delta=0, weights=np.zeros_like(self.y), max_iter=0, return_coef=True + ) + + assert_allclose(output, self.y, rtol=1e-14, atol=1e-14) + assert np.isnan(params['coef']).all() + + @pytest.mark.parametrize('strict_mask', (True, False)) + def test_zero_weights_masking(self, strict_mask): + """Ensures behavior when a zero-weight window occurs when using masking.""" + mask = np.zeros_like(self.y, dtype=bool) + fitter = self.algorithm_base(self.x, mask=mask, strict_mask=strict_mask) + if strict_mask: + context = pytest.raises(ValueError, match='A window had too few non-zero weights') + else: + context = pytest.warns(ParameterWarning, match='A window had too few non-zero weights') + + with context: + fitter.loess(self.y, delta=0, weights=np.zeros_like(self.y), max_iter=0) @pytest.mark.parametrize('use_threshold', (True, False)) def test_weight_masking(self, use_threshold): @@ -532,6 +665,13 @@ def test_weight_masking(self, use_threshold): with pytest.raises(AssertionError): super().test_weight_masking(use_threshold=use_threshold) + @ensure_deprecation(1, 5) + @pytest.mark.parametrize('conserve_memory', (True, False)) + def test_conserve_memory_deprecation(self, conserve_memory): + """Ensures a warning if emitted if conserve_memory is input.""" + with pytest.warns(DeprecationWarning, match='conserve_memory is deprecated'): + self.class_func(self.y, conserve_memory=conserve_memory) + class TestQuantReg(IterativePolynomialTester, RecreationMixin): """Class for testing quant_reg baseline.""" @@ -787,7 +927,13 @@ def test_median_absolute_value(values): mav_calc = polynomial._median_absolute_value(values) mav_actual = np.median(np.abs(values)) / 0.6744897501960817 - assert_allclose(mav_calc, mav_actual) + assert_allclose(mav_calc, mav_actual, rtol=1e-14, atol=1e-14) + + # also compare against scipy, set center to 0 to make the MAD into MAV + mav_scipy = stats.median_abs_deviation( + values, scale='normal', center=lambda *args, **kwargs: 0 + ) + assert_allclose(mav_calc, mav_scipy, rtol=1e-14, atol=1e-14) def test_loess_solver(): @@ -800,7 +946,7 @@ def test_loess_solver(): solved_coefs = polynomial._loess_solver(vander.T, y) - assert_allclose(solved_coefs, coefs) + assert_allclose(solved_coefs, coefs, rtol=1e-12, atol=1e-14) def test_determine_fits_simple():