ENH: Improve plateaus identification in `_fcutoff`

[physicien]

Description of the problem

At present, plateau detection in _fcutoff is carried out by imposing a narrow tolerance band around 0 on both the first (see the yellow [loose tolerance] and dashed orange [tight tolerance] regions in the autocorrelation plot) and second derivatives (see the blue region in the autocorrelation plot), and then designating as plateaus those points where these two criteria coincide.

I tested several alternative strategies that rely on features of the autocorrelation curve (such as inflection points or extrema of its first derivative), but these turned out to be overly sensitive to the instabilities and discontinuities that may appear in the autocorrelation plot. Consequently, a more robust approach is needed to reliably detect all plateaus.

Description of a possible solution or alternative

As @derb12 pointed out in an email exchange, instead of examining the first two derivatives, it may be more effective to consider something like the rolling standard deviation of the autocorrelation values. For this purpose, a rolling standard deviation is expected to yield a unimodal distribution of values, from which an appropriate threshold for detecting flat regions could be determined automatically using scikit-image’s threshold_triangle. That said, this approach still cannot guarantee selection of the correct plateau when multiple plausible plateaus exist. Even so, being able to robustly detect plateau regions at all would represent a significant improvement.

[derb12]

I spent a little time looking at this, and unfortunately the rolling standard deviation idea doesn't work very well in some cases, such as when a plateau has some amount of slope like the one you show above (thanks to the details you listed in fcut.md, I was able to understand how to make synthetic data that make varying amounts of plateaus). Plus, it loses all of the shape information that the derivative signs give. The auto-thresholding using scikit-image may still help you pick thresholds for the derivatives, though. All my testing code is below.

code

import matplotlib.pyplot as plt
import numpy as np
from pybaselines import Baseline
from pybaselines.utils import gaussian
from scipy import ndimage
from skimage.filters import threshold_yen, try_all_threshold


def make_data(x, baseline_type=0):
    """Creates the data and baseline where the endpoints are either zero or nonzero."""
    signal = (
        gaussian(x, 9, 100, 12)
        + gaussian(x, 6, 180, 5)
        + gaussian(x, 8, 350, 11)
        #+ gaussian(x, 15, 400, 18)
        + gaussian(x, 6, 550, 6)
        + gaussian(x, 13, 700, 8)
        + gaussian(x, 9, 800, 9)
        + gaussian(x, 9, 880, 7)
    )
    if baseline_type == 0:  # simple parabola that ends at 0 on both ends
        baseline = 2e-5 * (x - 500)**2 - 5
    elif baseline_type == 1:  # exponentially decaying baseline
        baseline = 10 - 10 * np.exp(-x / 600)
    elif baseline_type == 2:  # very complicated baseline
        # approximate a logistic baseline by integrating a Gaussian
        baseline = -np.cumsum(gaussian(x, 0.05, 400, 100)) + gaussian(x, 3, 800, 100) - 5
    else:  # a sinusoidal baseline
        baseline = 10 + 1 * np.sin(x / 50)
    noise = np.random.default_rng(0).normal(0, 0.05, len(x))
    y = signal + baseline + noise

    return y, baseline


#x = np.linspace(0, 1000, 1000)  # produces one plateau
x = np.linspace(0, 4000, 1000)  # produces two plateaus
#x = np.linspace(-4000, 1000, 1000)  # produces three plateaus

baseline_fitter = Baseline(x_data=x)
y, baseline = make_data(x, baseline_type=2)

half_window = 3  # for the rolling standard deviation; need to increase as x.size or cutoff_freqs increases
cutoff_freqs = np.geomspace(1e-4, 0.45, 30)

min_y = y.min()
y_transf = np.log(y - min_y + 1)
baseline_transf = np.log(baseline - min_y + 1)

lam_0 = 1 / abs(y_transf).sum()
lam_1 = 1 / abs(np.diff(y_transf, 1)).sum()
lam_2 = 1 / abs(np.diff(y_transf, 2)).sum()
alpha = 1

fig1, ax = plt.subplots()
fig2, ax2 = plt.subplots()
ax.plot(y_transf)
ax2.plot(y)
r_squares = []
min_mse = np.inf
for i, cutoff_freq in enumerate(cutoff_freqs):
    log_baseline, params = baseline_fitter.beads(
        y_transf, freq_cutoff=cutoff_freq, lam_0=alpha * lam_0, lam_1=alpha * lam_1,
        lam_2=alpha * lam_2, asymmetry=6
    )
    ax.plot(log_baseline, '--', label=cutoff_freq)
    fit = np.exp(log_baseline) + min_y - 1
    noise = y - fit  # omit subtracting signal since using log transform
    diff_noise = np.diff(noise, 1)
    durbin_watson = (diff_noise @ diff_noise) / (noise @ noise)
    r_sq = 0.25 * (2 - durbin_watson)**2
    r_squares.append(r_sq)
    ax2.plot(fit, '--', label=cutoff_freq)
    if (mse := ((log_baseline - baseline_transf)**2).sum()) < min_mse:
        best_idx = i
        min_mse = mse

ax.plot(baseline_transf, ':', label='true baseline')
ax2.plot(baseline, ':', label='true baseline')

ax.legend()
ax2.legend()
ax.set_title('log transformed')
ax2.set_title('transformed back')

r_squares = np.array(r_squares)
r_sq_prime = np.gradient(r_squares)
r_sq_double_prime = np.gradient(r_sq_prime)
if hasattr(ndimage, 'vectorized_filter'):  # introduced in scipy v1.16.0; better edge handling
    rolling_std = ndimage.vectorized_filter(
        r_squares,
        lambda vals, axis: np.std(vals, ddof=1, axis=axis),
        size=[2 * half_window + 1]
    )
else:
    # really goofy solution, not equal to vectorized_filter on the edges
    rolling_std = np.std(
        np.pad(
            np.lib.stride_tricks.sliding_window_view(r_squares, 2 * half_window + 1),
            [[half_window, half_window], [0, 0]], mode='reflect'
        ),
        axis=1, ddof=1
    )

plt.figure()
plt.plot(cutoff_freqs, r_squares, 'o-', label=r"$r^2$")
plt.plot(cutoff_freqs[best_idx], r_squares[best_idx], 'ro', label='best fit')
plt.plot(cutoff_freqs, r_sq_prime, 's--', label=r"$r^2$'")
plt.plot(cutoff_freqs, r_sq_double_prime, '.:', label=r"$r^2$''")
plt.plot(cutoff_freqs, rolling_std, '.:', label='rolling std')
plt.semilogx()
plt.legend()


### Thresholding based on rolling standard deviation ###
threshold = threshold_yen(rolling_std)
mask = rolling_std < threshold

plt.figure()
plt.plot(cutoff_freqs, r_squares, 'o-')
plt.plot(cutoff_freqs[mask], r_squares[mask], '.', label='std mask')
plt.legend()
plt.semilogx()

### Thresholding based on derivatives ###
deriv_threshold = threshold_yen(abs(r_sq_prime))
mask_deriv = abs(r_sq_prime) < deriv_threshold
deriv2_threshold = threshold_yen(abs(r_sq_double_prime))
mask_deriv2 = abs(r_sq_double_prime) < deriv2_threshold
total_mask = mask_deriv & mask_deriv2

plt.figure()
plt.plot(cutoff_freqs, r_squares, 's-', ms=10, label=r"$r^2$")
plt.plot(cutoff_freqs[mask_deriv], r_squares[mask_deriv], 's', label='first deriv mask')
plt.plot(cutoff_freqs[mask_deriv2], r_squares[mask_deriv2], 'o', label='second deriv mask')
plt.plot(cutoff_freqs[total_mask], r_squares[total_mask], '.', label='both masks')
plt.legend()
plt.semilogx()

# Try out different thresholding algorithms from scikit-image.
# Makes an "image" of the mask, where filled/black squares denote True on the mask.
# Using this on a few different combinations of autocorrelation plots, threshold_yen seems
# a bit better than threshold_triangle
fig, _ = try_all_threshold(rolling_std.reshape(1, -1), verbose=False)
fig2, _ = try_all_threshold(abs(r_sq_prime).reshape(1, -1), verbose=False)
fig.suptitle('rolling std')
fig2.suptitle('first derivative')

plt.show()

The only other thought I have is something I skimmed when I was working with L-curve optimization, which produces curves similar in shape to the autocorrelation plot. The paper "An Adaptive Pruning Algorithm for the Discrete L-Curve Criterion" by P.C. Hansen describes how to adaptively find the elbow in the "L" of the curve; maybe following similar logic to their idea would help to identify only plateaus. Given that derivative values, and hence any thresholding, for the autocorrelation plot are significantly dependant on data size, sampling of cutoff frequences, etc., basing the selection on known geometry of the curve might be an additional/alternative method to using the derivative?

[derb12]

Also to maybe step back a little, if you're not fixed on using beads, you could try mixture_model or loess as other baseline correction methods that work with both positive and negative peaks. mixture_model would have some dependence on data size for the lam value to use, but you could combine with optimize_extended_range to maybe work around that; loess should be fairly insensitive to data size since it's a rolling window method, so depending on how the peak widths vary with data size, you'd need to use either fraction or total_points (or combine with custom_bc like you did with beads). They don't provide the smoothed signal with the output like beads does, but I would think smoothing is a bit easier to standardize compared to baseline correction for these.

[physicien]

@derb12 Thanks for running these tests and sharing them. Interestingly, while you concluded that thresholding the rolling standard deviation (r-std) doesn’t work, I actually arrived at the opposite conclusion yesterday. Let me walk through my progress, and then I’ll return to your findings.

The following concerns the commit I just pushed. In short, I added the find_plateau2 function (in utils.py), which detects plateaus by applying a threshold to the r-std. Overall, it appears to perform well at picking up sufficiently large drops between distinct flat regions. As an example, here are a few representative plots (bottom-left panel; the purple hatched areas correspond to the detected plateaus):

As you can see, this approach appears to detect the plateaus quite effectively, including their endpoints, in the regions where the derivative begins to change gradually. In the 3rd figure, you can also see that when the vertical gap between two plateaus is too small, the method merges them and returns a single, continuous plateau spanning both. Nevertheless, my data seems to always show a noticeable difference in the “stability” of the r-std for different plateaus, which suggests that we might be able to exploit this feature afterward to further refine and shorten the list of candidate plateaus.

The only substantial failure mode I’ve observed occurs when the assumption that the distribution is unimodal is significantly violated. As an illustration, consider the following figure where the two plateaus situated between the two main peaks are not detected by threshold_triangle:

However, I handled this by first assessing whether the distribution is unimodal using Hartigan’s dip test (https://doi.org/10.1214/aos/1176346577), as implemented in the diptest package. To do so, I compute the p-value for the distribution, and if pval < 0.002, I compute the threshold on the "corrected" r-std (by subtracting local threshold values from it), instead of applying the procedure directly to the r-std. The cutoff 0.002 was chosen for convenience based on my own dataset, although, more formally, Freeman and Dale (https://doi.org/10.3758/s13428-012-0225-x) recommend using pval < 0.05. I also tested this latter criterion, and it appears to work well in practice.

To sum up my tests from the last two days, they suggest that the plateaus are correctly detected, and that the remaining task is to cut out regions where the r-std either fluctuates strongly or, conversely, becomes several orders of magnitude closer to zero (on the final plateau). For the strong fluctuations, my idea was to do this by examining the difference between the r-std and a “smoothed” version of it (for example, the orange curve, which corresponds to the local thresholds and is fairly smooth). Regions where this difference exceeds a certain cutoff would be discarded, which would cause some plateaus to become hatched. The resulting collection of small, leftover plateaus could then be eliminated using the _long_plateaus function.

As a side note, I also try using the r-mad, but it doesn’t seem effective for identifying the plateau of a clearly bimodal distribution, and it reduces the noise that I actually want to use for filtering out multiple regions.

[physicien]

Now, returning to your findings, I tested the synthetic data you suggested using my implementation, and here are the BEADS results:

BEADS

1 plateau

2 plateaus

3 plateaus

At first glance, it looks to me like the plateaus are correctly detected and that the methodology outlined at the end of my previous message could be viable. However, I see that your examples with 2 and 3 plateaus push the limits of my implementation because the features are clustered at either the very beginning or the very end of the signal. This issue arises mainly because I currently only trim the excess signal on the right-hand side.

Custom BEADS

1 plateau

2 plateaus

3 plateaus

When I instead apply the method using beads through custom_bc, your examples expose a major failure mode of this approach (or at least of my implementation of it): most of the regions containing features appear to have merged, so the algorithm recognizes them only as a single broad peak.

Now, regarding your conclusions, if the strategy based on the r-std ultimately proves ineffective, a hybrid method that combines the shape information from the derivative’s sign with the plateau regions of the r-std might be viable. Failing that, your suggestion of using an L-curve criterion could indeed offer a solution, though I’m not in a position to comment on it in more detail at this stage.

Regarding the possibility of moving to a different baseline-correction algorithm, I’m still open to it, but for now I’d prefer to keep refining how we use BEADS. It still seems very well suited to chromatography, as long as we can settle on a suitable value for fcut. Also, in its current state, my implementation is at least usable as an “assisted” method (and in many cases even effectively “automatic”), because the correct fcut can still be identified fairly easily from the autocorrelation plot. Looking further ahead, I’d be more inclined to add alternative methods alongside beads and custom_beads as complementary options within Weaselytics, rather than as direct replacements.

[derb12]

That's great, I gave up on the r-std too early it seems. I was a bit afraid I was wasting your time leading down dead-ends. And the check for unimodality coupled with the local thresholding is really smart. I wouldn't worry too much about my synthetic data, it was mostly just to see how the shape of the autocorrelation plot changes; the true benchmark is your real chromatography data, which it looks like this process works very well on now.

[physicien]

No worries, it was a very good idea you had. I'm used to disappearing down rabbit holes anyway, and hitting dead ends is part of the fun when you're hunting for solutions. About the idea of checking for unimodality and using local thresholding, I actually got it after reading this really informative page about thresholding.

Regarding your synthetic signals, I agree that we don’t need to worry about them too much, and that the real benchmark, ultimately, is my chromatographic data. While I’d like my method to be fully general, each type of measurement inevitably comes with its own reasonably acceptable assumptions. Let’s just hope that the possibility of a sinusoidal baseline in a chromatogram is not one of them hahaha.

By the way, for future reference, I’ve uploaded to Drive the images that show both the signal and its autocorrelation plot for all the signals in my dataset, and this for the present implementation. In case you are interested in checking out a set of real, highly varied signals, I will email you the link so you can look through them.

[physicien]

Progress update

Since my last update on this issue, I have experimented with a variety of methods to better pinpoint the plateau that contains the optimal value of fcut. To recap, the algorithm has so far been able, via thresholding on the rolling standard deviation (r-std), to reliably detect flat regions where r2 varies only slightly. We refer to these regions as “global plateaus,” as they are determined by the global parameter threshold. The current objective is to eliminate some of these plateaus, either entirely or in segments, to narrow the search to the specific plateau that includes the optimal fcut. In particular, the plateaus located at both ends of the autocorrelation plot are typically unstable or very close to 0, and thus must be filtered out as they correspond to a baseline that is either overly rigid or overly flexible.

Of the methods tested to date, most are either adaptations of the r-std approach that instead use a different statistic, such as the rolling median absolute deviation (r-mad) or the difference between r-std and r-mad, or they rely on manually defined thresholds applied to one or more of these quantities.

# Example
not_too_flat = rolling_std > 1.0E-06
plateaus = np.logical_and(plateaus, not_too_flat)

For completeness, it is worth mentioning a somewhat unconventional method that nearly worked. It stems from the observation that, when we compare the local thresholds of the r-std curve with the curve itself, the most stable regions correspond to segments where the two curves do not intersect and run almost in parallel. We can then efficiently locate all intersection points in much the same way as identifying inflection points of a function, by monitoring sign changes in their difference. A key benefit of this technique is that it does not rely on any explicit threshold, so it remains applicable regardless of how large or small the instabilities are. Afterward, we can discard all plateaus that are too short (for instance, those spanning fewer than 10 data points in an array of 1000), as these are likely spurious.

# Unconventional method
local_threshold = threshold_sauvola(rolling_std)
corrected = rolling_std - local_threshold

crossings = np.where(np.diff(np.sign(corrected)))[0]
plateaus[crossings] = False
plateaus = _long_segments(plateaus)

However, in a few specific situations, this method completely missed the region containing the optimal fcut value, which ultimately makes it unreliable.

plotted example of the unconventional method

• In the bottom-left plot, the blue line corresponds to rolling_std, while the orange line corresponds to local_threshold.

Breakthrough

Over the past few days, I did make some progress when I noticed that, in the earlier derivative-based method, the first and last plateaus (whose r-std is either unstable or nearly zero) were usually identified quite accurately. This gave me the idea to combine the derivative approach with rolling statistics, which resulted in the following function:

def find_plateaus(x, window=3, nbins=256, pval_cutoff=0.002):
    """
    NOTE: CHANGE pval_cutoff to 0.05 or 0.10 later
    """
    # Rolling statistics
    rolling_std = _rolling_std(x, window=window)
    rolling_mad = _rolling_mad(x, window=window)
    diff_std_mad = rolling_std - rolling_mad
    
    # Test if the distribution is unimodal (p=1)  
    _, pval = diptest(rolling_std)
    print(f"{'pval:':<20}{pval:0.4f}")

    # Find the threshold value
    if pval < pval_cutoff:
        # In case of significant multimodality
        local_threshold = threshold_sauvola(rolling_std)
        corrected = rolling_std - local_threshold
        threshold = threshold_triangle(corrected, nbins=nbins)
        plateaus = corrected < threshold
    else:
        threshold = threshold_triangle(rolling_std, nbins=nbins)
        plateaus = rolling_std < threshold
    print(f"{'Threshold:':<20}{threshold:0.4E}")
    
    # Discard plateaus at both ends
    ends = _flat_ends(x, diff_std_mad)
    plateaus = np.logical_and(plateaus, ~ends)

    # EXPERIMENTAL (mainly intended to prevent discarding plateaus beyond the
    # maximum of rolling_std when only a small noisy segment remains directly
    # before it)
    test = diff_std_mad < 5.0E-05
    plateaus = np.logical_and(plateaus, test)

    # Discard shorter plateaus
    plateaus = _long_segments(plateaus)
    
    # Discard regions beyond the max of rolling_std if nothing before it
    rstd_argmax = np.argmax(rolling_std)
    if plateaus[:rstd_argmax].any():
        plateaus[rstd_argmax:] = False

    return plateaus, ends, rolling_std, diff_std_mad

For reference, here is the implementation of _flat_ends, the function that detects the plateaus at both ends:

code

def _flat_ends(x, rdiff, smoothing_window=15, tol0=1.0E-03, tol1=1.0E-05,
               tol_rdiff=1.0E-04):
    """
    Identify the plateau regions at both ends of `x`, where `x` is the array of
    autocorrelation coefficients ,`r2`, as a function of a given parameter. In
    situations where the values of `x` increase again toward 1 after hitting
    their minimum (instead of continuing to decrease toward 0), this final
    region is also identified. Here, we assume that the left (right) side
    corresponds to parameter values for which the baseline is most rigid
    (flexible).

    Parameters
    ----------
    x : array-like, shape (N,)
        The y-values of the data.
    rdiff : array-like, shape (N,)
        The difference between the rolling standard deviation and the rolling
        median absolute deviation of `x`.
    smoothing_window : int, optional
        Standard deviation for the Gaussian kernel used to smooth the data.
        Default is 15.
    tol0 : float, optional
        Threshold applied to detect the first plateau of `x`, i.e., the region
        where the baseline exhibits the lowest flexibility. This parameter
        specifies the maximum allowed decrease from the value at which this
        plateau begins. Default is 1.0E-03.
    tol1 : float, optional
        Threshold applied to detect plateaus in the first derivative of the
        smoothed data. Default is 1.0E-05.
    tol_rdiff : float, optional
        Threshold applied to `rdiff` to locate the final plateau of `x`, i.e.,
        the segment where the baseline has the greatest flexibility. This
        parameter represents the maximum deviation of a point on that last
        plateau. Default is 1.0E-04.

    Returns
    -------
    ends : array-like, shape (N,)
        The boolean array whose `True` entries mark the regions to be excluded
        at both ends of `x`.

    """
    # Smoothed data and derivatives
    smooth_d0 = gaussian_filter1d(x, smoothing_window)
    smooth_d1 = np.gradient(smooth_d0)
    argmin_d1 = np.argmin(smooth_d1)

    # Initialization of the array to return
    ends = np.full(len(x), False, dtype=bool)

    # Proto-plateaus from d1
    tight_d1_flats = find_flat(smooth_d1, tol1)
    tight_continuous = continuous_ranges(tight_d1_flats)

    # Initial plateau
    starting_r2 = np.median(smooth_d0[tight_continuous[0]])   # mean or median?
    starting_end = np.where(
            np.absolute(starting_r2 - x[:argmin_d1]) < tol0)[0][-1]
    starting_plateau = starting_end + 1 
    ends[:starting_plateau] = True

    # Final plateau (or climbing final region)
    last_r2 = len(x) - 1
    if rdiff[last_r2] <= tol_rdiff:
        ending_r2 = np.where(rdiff > tol_rdiff)[0][-1]
        last_r2 = ending_r2 - 1
    else:
        last_r2 = np.argmin(x)
    ends[last_r2:] = True
    return ends

Results

To demonstrate how this new approach performs, the following examples are provided, highlighting the main scenarios identified in my dataset:

Example 1: one large drop

• In the bottom-left plot, the blue line corresponds to rolling_std, while the orange line corresponds to diff_std_mad.

Example 2: two large drops

• In the bottom-left plot, the blue line corresponds to rolling_std, while the orange line corresponds to diff_std_mad.

Example 3: short signal with very sharp low-amplitude peaks

• In the bottom-left plot, the blue line corresponds to rolling_std, while the orange line corresponds to diff_std_mad.

Example 4: significant bimodality of the r-std

• In the bottom-left plot, the blue line corresponds to rolling_std, while the orange line corresponds to diff_std_mad.

Example 5: very broad first plateau

• In the bottom-left plot, the blue line corresponds to rolling_std, while the orange line corresponds to diff_std_mad.

Example 6: case that could be challenging to improve

• In the bottom-left plot, the blue line corresponds to rolling_std, while the orange line corresponds to diff_std_mad.

In summary, a hybrid strategy that combines derivatives to detect flat regions at the boundaries with rolling statistics to remove segments where the variability is too high appears promising. That said, I still find it difficult to fully grasp how the various criteria within find_plateaus interact with one another. I am also skeptical about the reliability of parts of the code that depend on manually chosen thresholds. Lastly, the plateaus produced by the current method often appear overly fragmented, and I am unsure how to mitigate this without beginning to accept regions that are clearly unreasonable.