adding peak_voigt model - #743
Conversation
|
Purpose of this PR is to add pseudo Voigt peak shape model. If someone could help with the RUFF issue, I would be grateful ! |
|
@marimperorclerc - I've applied the ruff fixes. If you would like to not have this type of issue in the future, please follow the instructions at https://github.com/orgs/SasView/discussions/3171#discussioncomment-17148772 to get ruff and pre-commit set up. |
| sigma = hwhm / cste | ||
| intensity = (wf * (1 / (1 + ((q - q0)**2.0 / hwhm**2.0)))) + \ | ||
| ((1.0 - wf) * np.exp((-0.5 * (q - q0)**2.0) / (sigma**2.0))) | ||
| return intensity |
There was a problem hiding this comment.
You can use Voigt rather than pseudo-Voigt:
def voigt(x, sigma, gamma):
"""
Return the voigt function, which is the convolution of a Lorentz
function with a Gaussian.
:Parameters:
gamma : real
The half-width half-maximum of the Lorentzian
sigma : real
The 1-sigma width of the Gaussian, which is one standard deviation.
Ref: W.I.F. David, J. Appl. Cryst. (1986). 19, 63-64
Note: adjusted to use stddev and HWHM rather than FWHM parameters
"""
# wofz function = w(z) = Fad[d][e][y]eva function = exp(-z**2)erfc(-iz)
from scipy.special import wofz
# TODO: if sigma == 0: return Lorentzian
z = (x + 1j * gamma) / (sigma * np.sqrt(2))
V = wofz(z) / (np.sqrt(2 * pi) * sigma)
return V.real
For a 1000 point function on my mac this takes 47.5 μs vs 12.5 μs for the pseudo-Voigt function.
There was a problem hiding this comment.
if we are going full Voigt why not just
from scipy.special import voigt_profile
V = voigt_profile(x, sigma, gamma)
Return VThat said, Pseudo Voigt seems to be quite popular as @marimperorclerc points out, probably due to the factor of 4 you measured. Going where the users are does not seem like a terrible plan to me? Perhaps we could add the full Voigt as an option to minimize the number of new models? Though I think that may get too complicated.
There was a problem hiding this comment.
I wrote voigt before voigt_profile was available. I didn't know they added it to scipy (Dec. 2019).
The pseudo-Voigt was popular because it was readily available. The linear combination of Gauss and Lorentz distributions was easier to calculate than the w(z) function. Other than the rare situation where you have pure Gaussian and Lorentzian peak broadening from different parts of the sample adding incoherently, I don't know that it has any physical meaning.
The parameterization of pseudo-Voigt using HWHM was a further convenience, being something that can be read directly off the graph.
Now that we can fit the Gaussian and Lorentzian parameters of the Voigt profile directly, that seems more useful than continuing with pseudo-Voigt based on peak width.
There was a problem hiding this comment.
I would argue that that decision is "above my pay grade." It really needs input from the community that uses it. That said I would argue a Pseudo Voigt is better than none at all and currently we only have a PR for a Pseudo Voigt - though as we just showed above writing a true Voigt would be trivial should someone be so inclined? Then the question would become a different one: do we need both?
| # cste = 1.17741 | ||
| sigma = hwhm / cste | ||
| intensity = (wf * (1 / (1 + ((q - q0)**2.0 / hwhm**2.0)))) + \ | ||
| ((1.0 - wf) * np.exp((-0.5 * (q - q0)**2.0) / (sigma**2.0))) |
There was a problem hiding this comment.
Rather than the "\" line extender use parentheses around the entire expression so it can break across multiple lines without the backslash.
Move the "+" operator to the beginning of the second line.
| half-maximum) for the Lorentzian and sigma = HWHM / 1.17741 for the | ||
| Gaussian, where sigma is the standard deviation of the Gaussian. In | ||
| other words, the widths of the Lorentzian and the Gaussian have been | ||
| coupled for convenience of parameterisation. |
There was a problem hiding this comment.
Convenient for whom? The width of the Gaussian and Lorentzian are easier to interpret physically. I don't think the optimizer will care which pair it is fitting. The only advantage seems to be that it is easier to guess the initial value of HWHM from the graph.
There was a problem hiding this comment.
I note that the reference whence this was taken (according to the reference cited here) does exactly this (create a single HFHM term for both line shapes) as does the marketplace model. It also reduces the number of fit parameters. Not sure "convenience of parameterization" is the correct wording. I also don't know if this coupling is general enough for most scattering use cases or if there is a real call for decoupling them.
There was a problem hiding this comment.
The number of fit parameters is the same: center, width, gauss:lorentz ratio vs. center, σ (gauss), γ (lorentz)
There was a problem hiding this comment.
Except that there is only the center and width no ratio right? Unless you are talking about weighting ratio and are comparing the Voigt to Pseudo Voigt? I'm confused.
There was a problem hiding this comment.
The w_f parameter gives the relative portion of gaussian to lorentzian. The fit can adjust center, width and proportion.
The gaussian and the lorentzian have the same peak_hwhm value, so adjusting w_f effectively changes the tail weight without changing the width.
I haven't checked in detail that every choice of (σ, γ) for voigt can be approximated by some choice of peak_hwhm and w_f and vice versa, but the parameter count is the same.
| This pseudo-Voigt peak function is a weighted linear summation of | ||
| Lorentzian (L) and Gaussian (G) peak shapes. | ||
| It is a popular function for modelling peak shape. | ||
| It can be tailored to any specific peak shape and it can also produce a peak shape with asymmetry. |
There was a problem hiding this comment.
Both the gaussian and the lorentzian are symmetric about peak_pos, and so is the linear combination. It make look asymmetric due to log q scale, or because the resolution function Δq is increasing with q.
Mention that the voigt function is a convolution
Add a note that the instrument resolution function contributes to the convolution. That is, if the resolution is a gaussian of approximately constant width Δq then peak_hwhm parameter value.
There was a problem hiding this comment.
Agreed -- asymmetric Pseudo Voigt is a thing but requires some modification such as a q dependent FWHM. see for example https://analyticalsciencejournals.onlinelibrary.wiley.com/doi/10.1002/sia.5521 or chrome-extension://efaidnbmnnnibpcajpcglclefindmkaj/https://arxiv.org/pdf/1804.06083. This function does not seem to be that.
On the other hand, I see that the reference from which @smk78 took the equation he placed in the marketplace and from which this seems to be derived, is specifically about an asymmetric Pseudo Voigt which does exactly what I state above (https://www.sciencedirect.com/science/article/pii/S0924203108000453). I it is probably also why that paper also couples the width of the two peaks?
However I suspect that is more complexity than we may want here. At least for now I'd recommend keeping it simple as done in the marketplace did, and if there is a clamoring for allowing asymmetry we can add that later (the asymmetry parameter should be adjustable so that it has no effect)?
There was a problem hiding this comment.
Regarding convolution with resolution... isn't that generically true for any peak function? Is there a standard blurb we use for all those?
There was a problem hiding this comment.
Specifically, I'm asking that the clause "it can also produce a peak shape with asymmetry" be removed from the documentation. No combination of parameters will lead to an asymmetric peak.
There was a problem hiding this comment.
Regarding peak broadening through resolution it is not mentioned in other models. It seemed relevant here because the measured peak HWHM will be larger than the fitted peak_hwhm value.
For the other models:
- broad peak uses
correlation_lengthbut the arbitrary power changes the width - gaussian peak uses
sigmawith measured FWHM = 2.35σ + resolution broadening - lorentz uses
cor_lengthwith measured HWHM = γ + resolution broadening - two lorentzian uses
lorentz_length_[12]but the arbitrary power changes the width
| Authorship and Verification | ||
| ---------------------------- | ||
|
|
||
| * **Author:** Steve King **Date:** 24 June 2020 |
There was a problem hiding this comment.
There should only be one Author line. Otherwise it gets confusing. In principle I would say that if this was just a minor edits to the code in the marketplace, the author normally would remain the original author with the new authors as "last modified by."
On the other hand, if the code was inspired by the marketplace but significantly rewritten here then I would just replace the marketplace authors with both your names.
I do not think there is precedent in sasmodels yet, but one could add a line at the top called "original author." This could then be used for the marketplace author and authors would be the two of you (normally one line with "and") and a date which could just be "June 2026" if a specific date doesn't make sense.
| 1. L A Feigin, D I Svergun, G W Taylor | ||
| Structure Analysis by Small-Angle X-ray and Neutron Scattering | ||
| Springer (1987) |
There was a problem hiding this comment.
Is this reference truly relevant to this model? I don't find any Voigt or Pseudo Voigt in that book. Could you had a page number if it is? and if not maybe remove. It is an excellent textbook on scattering but Voigt peaks are fundamentally just math functions with no intrinsic scattering knowledge needed to write this?
There was a problem hiding this comment.
No quality gates enabled for this code.
See analysis details in CodeScene
Quality Gate Profile: Custom Configuration
Install CodeScene MCP: safeguard and uplift AI-generated code. Catch issues early with our IDE extension and CLI tool.
|
Oops... I accidentally clicked on merge with master. @marimperorclerc, you will need to pull before making more changes. Sorry! |
| peak_pos: position of the peak | ||
| peak_hwhm: HWHM of the peak | ||
| """ | ||
|
|
||
| with errstate(divide='ignore'): | ||
| L = Ipeak(q, w_f, peak_pos, peak_hwhm) | ||
|
|
||
| return L | ||
|
|
There was a problem hiding this comment.
Helper functions are certainly allowed, but in this case it seems completely redundant and just makes it slightly more convoluted to read. Just rename Ipeak to I(q) and delete this function which does nothing at all from what I can tell? That function already returns a sensible variable (intensity instead of L) and doesn't try to "ignore divide by zero which cannot happen from what I can tell s the denominator is always >1?
The only other thing that will need to be done when changing the Ipeak function name to Iq is to also change the parameter list names. For some reason Ipeak renames the variables currently.
butlerpd
left a comment
There was a problem hiding this comment.
The Voigt function is a commonly used function which was never implemented in sasmodels and requires relying on a marketplace plugin. As such this is a great addition and it would be nice to see it included in the next point release. I've tested functionally and it performs as advertised.
Despite the voluminous comments, all of which are useful, The important issues that should be addressed before a merge IMO are fairly simple and I list them here.
- Correct the documentation about asymmetry (unless of course the choice is to rewrite the function to include asymmetry).
- Correct the references and authorship as appropriate (see comments/questions)
- Remove the Ipeak (actually it becomes Iq and the old Iq code is just deleted) -- or explain why it is really necessary to keep separate.. in which case other fixes will be needed in that function IMO.
|
Thank you for reviewing this model ! |
|
I would hazard a guess that the pseudo-Voigt gained traction amongst crystallographers in those heady days when processors were slow and memory was expensive and limited... I recall it was implemented in the CCP13 XFIX program at the start of the 90's. |
No description provided.