fwhm and height calculation of BreitWignerModel mistakenly removed?

[hmoldenhauer]

I'm using the BreitWignerModel of the lmfit package and i was wondering, why the fwhm and height calculation have been removed in the recent lmfit-release by commit d086e62

For me there seems to be no reason to do so, as there are analytical solutions to the parameters (see here for example). This was already discussed, but unfortunately the parameters have been removed.

I even did the math to verify the solutions suggested by Austin Fox and came to the very same analytical solutions.

So here's the question:

Could the fwhm and height calculation be included again in future lmfit versions?

If so, could i help anyhow to get this fixed in the near future?

I'm curious to read your responses

[Matt Newville]

On Tue, Jul 14, 2020 at 7:06 AM hmoldenhauer <henning.m...@gmail.com> wrote:

I'm using the BreitWignerModel of the lmfit package and i was wondering, why the fwhm and height calculation have been removed in the recent lmfit-release by commit d086e62

Because that calculation is wrong.

For me there seems to be no reason to do so, as there are analytical solutions to the parameters (see here for example). This was already discussed, but unfortunately the parameters have been removed.

 

I even did the math to verify the solutions suggested by Austin Fox and came to the very same analytical solutions.

The attached plot shows the result of

   x = np.linspace(-100, 100, 2001)

   y = breit_wigner(x, center=1, sigma=3, q=1.5)

   plot(x, y)

The FWHM calculation fo this gives  ~37.  

So here's the question:

Could the fwhm and height calculation be included again in future lmfit versions?

If so, could i help anyhow to get this fixed in the near future?

I don't think that calculation is right. 

--Matt

[hmoldenhauer]

Am Mittwoch, 15. Juli 2020 00:44:18 UTC+2 schrieb Matt Newville:

On Tue, Jul 14, 2020 at 7:06 AM hmoldenhauer <henning.m...@gmail.com> wrote:

I'm using the BreitWignerModel of the lmfit package and i was wondering, why the fwhm and height calculation have been removed in the recent lmfit-release by commit d086e62

Because that calculation is wrong.

i plotted several bwf distributions (see attached image), as visible, the calculations are correct, except for the region around [-2, 2]. This is based in the root of the fwhm calculation (1/(q**2 - 2)). So for small q a different solution has to be developed, but for q beyond this region the calculations are correct. I try to find a solution for the exceptional region.

The plots were plotted using:

import numpy as np
import matplotlib.pyplot as plt
from lmfit.lineshapes import breit_wigner

qs = [-5, -2, -np.sqrt(2), -0.5, 0, 0.5, np.sqrt(2), 2, 5]

sigma = 3.0
mu = 0

fig, ax = plt.subplots(3, 3)

for i, q in enumerate(qs):
    fwhm = np.sqrt(q**2 * sigma**2 * (2+q**2)) / (q**2-2)

    x_min = (q * sigma - 2 * mu + q**2 * mu) / (q**2 -2) - fwhm/2
    x_max = (q * sigma - 2 * mu + q**2 * mu) / (q**2 -2) + fwhm/2

    x = np.linspace(-100, 100, 2001)

    y = breit_wigner(x, center=mu, sigma=sigma, q=q)

    print(fwhm)

    ax[int((i-i%3)/3)][i%3].plot(x,y)
    if abs(q) == np.sqrt(2):
        pass
    else:
        ax[int((i-i%3)/3)][i%3].plot(x_min, breit_wigner(x_min, center=mu, sigma=sigma, q=q), 'r.')
        ax[int((i-i%3)/3)][i%3].plot(x_max, breit_wigner(x_max, center=mu, sigma=sigma, q=q), 'r.')
    ax[int((i-i%3)/3)][i%3].set_title(f'breit_wigner(x, c = {mu:3.1f}, s = {sigma:3.1f}, q = {q:3.1f}) fwhm = {fwhm:1.1e}')

plt.show()

[Matt Newville]

On Wed, Jul 15, 2020 at 4:20 AM hmoldenhauer <henning.m...@gmail.com> wrote:

Am Mittwoch, 15. Juli 2020 00:44:18 UTC+2 schrieb Matt Newville:

On Tue, Jul 14, 2020 at 7:06 AM hmoldenhauer <henning.m...@gmail.com> wrote:

I'm using the BreitWignerModel of the lmfit package and i was wondering, why the fwhm and height calculation have been removed in the recent lmfit-release by commit d086e62

Because that calculation is wrong.

i plotted several bwf distributions (see attached image), as visible, the calculations are correct, except for the region around [-2, 2].

Well, I think that means the calculations are wrong.   I would also question whether the calculations are "correct" when |q| > 2.  Try q=5 and sigma=50.  The value for FWHM is a decent estimate, but off by several percents.   You are certainly free to use that calculation in your own code, but I'm pretty sure we do not want to report those values as FWHM.

If there was a better calculation of FWHM, we would be willing to use it.  But the Breit-Wigner distribution isn't really a function with a single peak that goes to zero at "very many sigmas from the center".  The meanings of "width" and "max" in "full width at half max" are confusing (and the meaning of "height" is also confusing), even in the plots you post.    So, I kind of doubt that the even could be a calculation for FWHM that works.

Basically, "sometimes wrong" means we should probably not be reporting it as FWHM. 

--Matt

This is based in the root of the fwhm calculation (1/(q**2 - 2)). So for small q a different solution has to be developed, but for q beyond this region the calculations are correct. I try to find a solution for the exceptional region.

The plots were plotted using:

import numpy as np
import matplotlib.pyplot as plt
from lmfit.lineshapes import breit_wigner

qs = [-5, -2, -np.sqrt(2), -0.5, 0, 0.5, np.sqrt(2), 2, 5]

sigma = 3.0
mu = 0

fig, ax = plt.subplots(3, 3)

for i, q in enumerate(qs):
    fwhm = np.sqrt(q**2 * sigma**2 * (2+q**2)) / (q**2-2)

    x_min = (q * sigma - 2 * mu + q**2 * mu) / (q**2 -2) - fwhm/2
    x_max = (q * sigma - 2 * mu + q**2 * mu) / (q**2 -2) + fwhm/2

    x = np.linspace(-100, 100, 2001)
    y = breit_wigner(x, center=mu, sigma=sigma, q=q)

    print(fwhm)

    ax[int((i-i%3)/3)][i%3].plot(x,y)
    if abs(q) == np.sqrt(2):
        pass
    else:
        ax[int((i-i%3)/3)][i%3].plot(x_min, breit_wigner(x_min, center=mu, sigma=sigma, q=q), 'r.')
        ax[int((i-i%3)/3)][i%3].plot(x_max, breit_wigner(x_max, center=mu, sigma=sigma, q=q), 'r.')
    ax[int((i-i%3)/3)][i%3].set_title(f'breit_wigner(x, c = {mu:3.1f}, s = {sigma:3.1f}, q = {q:3.1f}) fwhm = {fwhm:1.1e}')

plt.show()

 

For me there seems to be no reason to do so, as there are analytical solutions to the parameters (see here for example). This was already discussed, but unfortunately the parameters have been removed.

 

I even did the math to verify the solutions suggested by Austin Fox and came to the very same analytical solutions.

The attached plot shows the result of

   x = np.linspace(-100, 100, 2001)

   y = breit_wigner(x, center=1, sigma=3, q=1.5)

   plot(x, y)

The FWHM calculation fo this gives  ~37.  

So here's the question:

Could the fwhm and height calculation be included again in future lmfit versions?

If so, could i help anyhow to get this fixed in the near future?

I don't think that calculation is right. 

--Matt

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/e752869c-6d3d-468f-b713-b21b415c37b3o%40googlegroups.com.

[Eugene Grey]

This is an old thread, but I thought I could help clear some things up. Austin Fox's derivation is incorrect because it assumes that the maximum occurs at x = center, which is only true if q is infinite. The actual extrema for the Breit-Wigner-Fano (BWF) function occur at:
x_1 = center + sigma / (2 * q), giving BWF(x=x_1) = amplitude * (1 + q**2)
x_2 = center – q * sigma / 2, giving BWF(x=x_2) = 0.

Using amplitude * (1 + q**2) as the maximum, then: fwhm = abs(sigma * (1 + q**2) / (-1 + q**2)),
where the absolute value is needed because x_fwhm_1 > x_fwhm_2 for q > 1 and x_fwhm_1 < x_fwhm_2 for q < 1. For q = 1, there is only a single x value at which BWF = max/2, because the second value is asymptotically approached as abs(x) goes to infinity, making fwhm undefined for q = 1.

Rather than having lmfit add back in the height and fwhm calculations, I would recommend looking at how the BWF function is defined/used within your field, so that your fit results can correctly be communicated to others. For example, in Raman spectroscopy, the BWF function is often defined as:

BWF_alternate = H * (1 + (x - center) / (q * width))**2 / (1 + ((x - center) /  (width))**2)

where H = amplitude * q**2 and width = sigma / 2 compared to the lmfit definition. In this definition of the BWF function, the H variable is the height of the upward-facing peak compared to the BWF’s baseline (so the max of BWF_alternate = H * (1 + 1 / q**2), and the baseline = H / q**2 = amplitude), or in other words, as the max height of the function when q is infinite.

The Origin software also uses the BWF_alternate definition, and further defines the fwhm using the maximum of the BWF function, so that ultimately their calculated fwhm = abs(2 * width * (1 + q**2) / (-1 + q**2)), which is the same result as the fwhm given above with width = sigma / 2. However, for abs(q) <= 1, Origin’s output differs from that calculation by an order of magnitude, so I’m not sure what they’re doing for that case.

So, in summary, height could be defined as the maximum of the function or as the height of the upward facing peak, while the fwhm is usually consistently calculated using the x-values where y = max/2, but causes issues for abs(q) <= 1. Due to these ambiguities and issues, I think it was correct for them to not be included for the BWF model in lmfit since users can do any further calculations when appropriate.