Add possibility to modify initial guesses / limits of parameters of analyze_peak() function (DRT)

[nathan-kuehr]

Hey, the minimising algorithms that are used for the peak analysis usually work well, but actually are highly dependent on the data you put in. If the sigma and / or alpha values are too small / close to ±1, then the slope of the skew normals may attain numerical zero at the next data point, i.e. we have vanishing gradients in the optimiser. This is problematic as non-zero DRT datapoints which are close by won't have any effects now. The problem is illustrated below. I plotted the estimated polarisation as well.

For the second image, I modified the minimal possible sigma value to 1e-3, which disallows for too small sigmas. We can see that now, it correctly fits the skew normal to the data points. In the first one, the closest local minimum seems to be a value of sigma which is very very small. The notion of "too small", however, is dependent on the spacing / sampling of the test data.

That's why it would be nice to be able to adapt these parameters, without having to manually modify the local copy of the library. Maybe via a kwargs that is passed to the _generate_parameters function.
Or even cooler, adaptive sigmas that depend on the spacing of the data points.

[nathan-kuehr]

And thank you very much for the package, it is very well done :)

[vyrjana]

Thanks for the feedback. If you have an impedance or DRT spectrum (simulated or real) that I can test against and you would be willing to share, then that would be helpful. I'm working on an overhaul of parts of both pyimpspec and DearEIS so this feedback arrived at a good time since I haven't yet started working on the DRT part of the codebase.

Just as a side note and for future reference, you don't necessarily need to modify a local copy of the pyimpspec package. You can also use monkey patching in your script (or your own package) to add a wrapper around the _generate_parameters function at runtime in order to adjust the parameters. This monkey patching approach makes it easier to re-use your fix across multiple virtual environments and across versions of pyimpspec. For example:

import pyimpspec

function_to_wrap = pyimpspec.analysis.drt.peak_analysis._generate_parameters

def wrapper(*args, **kwargs):
    retval = function_to_wrap(*args, **kwargs)
    parameters = retval[0]

    # Do what you need here
    for key in parameters:
        if "sigma_" in key:
            parameters[key].set(min=1e-3)

    return retval

# Default function
peaks = [(0.3, 0.8), (0.8, 1.3)]
before = pyimpspec.analysis.drt.peak_analysis._generate_parameters(peaks=peaks, disallow_skew=False)
print(f"{before=}")

# Wrapped function
pyimpspec.analysis.drt.peak_analysis._generate_parameters = wrapper
after = pyimpspec.analysis.drt.peak_analysis._generate_parameters(peaks=peaks, disallow_skew=False)
print(f"{after=}")

[nathan-kuehr]

Thanks for the monkey patching advice, this is actually a way more elegant way of doing so!

Sure thing, I attached the data of the three samples that I measured. For context, these are EIS measurements from PEDOT:PSS-coated ITO electrodes in DPBS, presenting important reflective Warburg like behaviour for low frequencies.

• 20260306_PEDOT_PSS_on_ITO_MG2_NoAct_S42: A (control) sample where the skew normal fitting went well !
• 20260306_PEDOT_PSS_on_ITO_MG2_Act_S92 A sample where σ vanishes for the peak at ~ 1.2e-3 s
• 20260306_PEDOT_PSS_on_ITO_MG2_NoAct_S41 A sample where α gets really close to 1 for the peak at ~ 8.69e-4 s, presenting the vanishing gradient for the left neighbouring data points.

I attached both the raw data as well as the calculated DRT values. In case you want to work with the raw data, I excluded all frequencies slower than 15 Hz in the DRT calculation, and used tr-nnls.

processed_drt_data.py

raw_data.zip

[nathan-kuehr]

Dear Ville,
thanks again for the monkey-patching suggestion ! In case it's useful to you or anyone running into the same issue, here's the workaround that ended up working for me.

In the wrapper around _generate_parameters instead of fitting σ and α directly, I introduce two reparametrized parameters widthl and widthr, which are the left and right flank distances from the peak center (by default: 3σ distance). These become the free fit parameters, while σ and α become derived parameters.
Each width is then bound from below by the distance to the nearest neighbouring data point, which prevents a flank from collapsing onto the peak and thus avoids the vanishing gradients. Works like a charm :)

Here's my code:

def _make_generate_parameters_wrapper(taus: np.ndarray, sigma_multiplier: float = 3.0):
    # Transform to logged x axis as in analyze_peaks()
    x = np.log(taus)
    x -= x.min()
    x /= x.max()

    def patch(peaks, disallow_skew):
        parameters, num_variables = _original_generate_parameters(
            peaks, disallow_skew
        )

        for i, (px, _) in enumerate(peaks):
            neighbour_dist_left = (px - x[x < px][-1]) if np.any(x < px) else x[1] - x[0]
            neighbour_dist_right = (x[x > px][0] - px) if np.any(x > px) else x[-1] - x[-2]

            sigma_init = parameters[f"sigma_{i}"].value
            # Alpha is 0 initialized, so ignore

            widthl_init = max(sigma_multiplier * sigma_init, neighbour_dist_left)
            widthr_init = max(sigma_multiplier * sigma_init, neighbour_dist_right)

            if disallow_skew:
                parameters.add(f"widthl_{i}", value=widthl_init, min=max(neighbour_dist_left, neighbour_dist_right))
                parameters.add(f"widthr_{i}", expr=f"widthl_{i}") # symmetric, so follow left width
            else:
                parameters.add(f"widthl_{i}", value=widthl_init, min=neighbour_dist_left)
                parameters.add(f"widthr_{i}", value=widthr_init, min=neighbour_dist_right)

            # Deactivate alpha and sigma
            parameters[f"sigma_{i}"].set(
                vary=False,
                expr=f"(2.0/{sigma_multiplier}) * (widthl_{i} * widthr_{i}) / (widthl_{i} + widthr_{i})",
            )
            parameters[f"alpha_{i}"].set(
                vary=False,
                expr=f"(widthl_{i} - widthr_{i}) / (widthl_{i} + widthr_{i})",
            )

        return parameters, num_variables

    return patch