Skip to content

adding peak_voigt model #743

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
marimperorclerc wants to merge 2 commits into
base: master
Choose a base branch
from
Open

adding peak_voigt model #743

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d295689
adding peak_voigt model
marimperorclerc Jul 7, 2026
44f9c50
ruff fix
krzywon Jul 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
149 changes: 149 additions & 0 deletions sasmodels/models/peak_voigt.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,149 @@
r"""
This model describes a pseudo-Voigt shaped peak on a flat background.

Definition
----------

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.

@pkienzle pkienzle Jul 7, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 $(L_γ \star G_σ)(q)$. pseudo-Voigt approximates this with a linear combination of a lorentzian and a gaussian.

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 $f(q) = (L_γ \star G_σ \star G_{Δq})(q) = (L_γ \star G_σ')(q)$ with σ'=√(σ² + Δq²) where σ is the width returned from the fit. The fitted σ is intrinsic to the sample. You will need to subtract the modified resolution width Δq/√(2 \ln 2) from the HWHM when guessing the peak_hwhm parameter value.


The scattering intensity $I(q)$ is calculated as

.. math::

I(q) = scale \cdot \left[ w_f \cdot I(q)_L + (1 - w_f) \cdot I(q)_G \right] + background

where $w_f$ is a weighting factor and

.. math::

I(q)_L = \frac{1}{1 + \left( \frac{q - q_0}{HWHM} \right)^2}

I(q)_G = \exp\left[ -\frac{1}{2} (q - q_0)^2 / \sigma^2 \right]

The peak is taken to be centered at $q_0$ with a HWHM (half-width
half-maximum) of $1.17741\,\sigma$, 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:

.. math::

\sigma = HWHM / \sqrt{2 \ln 2} = HWHM / 1.17741

When $w_f = 1$ a Lorentzian peak is returned, and when $w_f = 0$ a
Gaussian peak is returned.

For 2D data the scattering intensity is calculated in the same way as 1D,
where the $q$ vector is defined as

.. math::

q = \sqrt{q_x^2 + q_y^2}


Validation
----------

The pseudo-Voigt peak reduces exactly to a pure Lorentzian for $w_f = 1$
and to a pure Gaussian for $w_f = 0$; both limits were checked against their
analytic values (see tests section at the end).
The full pseudo-Voigt shape has also been compared, for identical
parameters, against a slightly different SasView implementation (https://marketplace.sasview.org/models/127/)
of the same function and gives the same result.


References
----------

1. L A Feigin, D I Svergun, G W Taylor
Structure Analysis by Small-Angle X-ray and Neutron Scattering
Springer (1987)

2. Aaron L. Stancik, Eric B. Brauns
A simple asymmetric lineshape for fitting infrared absorption spectra
Vibrational Spectroscopy 47 (2008) 66-69


Authorship and Verification
----------------------------

* **Author:** Steve King **Date:** 24 June 2020

* **Authors:** Marianne Imperor-Clerc (marianne.imperor@cnrs.fr)
Anirban Mandal (mandalanirban2023@gmail.com)

* **Last Modified by:** Anirban Mandal **Date:** 06 July 2026

* **Last Reviewed by:** Steve King **Date:**

"""

import numpy as np
from numpy import errstate, inf

name = "peak_voigt"
title = "Single pseudo-Voigt peak"
description = """\
I(q) = scale*peak + background
"""

category = "shape-independent"

parameters = [["w_f", "", 0.8, [0, 1], "", "lorentzian/gaussian weighting factor"],
["peak_pos", "1/Ang", 0.05, [0, inf], "", "Position of the peak"],
["peak_hwhm", "1/Ang", 0.01, [0, 1], "", "HWHM of the peak"]]


def Ipeak(q, wf, q0, hwhm):
"""
When $w_f$ = 1 a Lorentzian peak is returned, and when $w_f$ = 0 a
Gaussian peak is returned.

The peak is taken to be centered at $q_0$ with a HWHM (half-width
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.

@pkienzle pkienzle Jul 7, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

"""
cste = np.sqrt(2 * np.log(2))
# 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)))

@pkienzle pkienzle Jul 7, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

return intensity

@pkienzle pkienzle Jul 7, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.



def Iq(q, w_f, peak_pos, peak_hwhm):
"""
w_f: weighting coefficient in the pseudo-Voigt peak function;
w_f = 1 for a Lorentzian and w_f = 0 for a Gaussian peak.
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

Iq.vectorized = True # Iq accepts an array of q values

tests = [
# pure Lorentzian (w_f = 1): peak centre, half-width, and 2 x HWHM
[{"scale": 1.0, "background": 0.0, "w_f": 1.0,
"peak_pos": 0.05, "peak_hwhm": 0.01}, 0.05, 1.0],
[{"scale": 1.0, "background": 0.0, "w_f": 1.0,
"peak_pos": 0.05, "peak_hwhm": 0.01}, 0.06, 0.5],
[{"scale": 1.0, "background": 0.0, "w_f": 1.0,
"peak_pos": 0.05, "peak_hwhm": 0.01}, 0.07, 0.2],
# pure Gaussian (w_f = 0): half-width is 0.5 by definition, 2 x HWHM = 1/16
[{"scale": 1.0, "background": 0.0, "w_f": 0.0,
"peak_pos": 0.05, "peak_hwhm": 0.01}, 0.06, 0.5],
[{"scale": 1.0, "background": 0.0, "w_f": 0.0,
"peak_pos": 0.05, "peak_hwhm": 0.01}, 0.07, 0.0625],
# mixed pseudo-Voigt (w_f = 0.8) away from the centre
[{"scale": 1.0, "background": 0.0, "w_f": 0.8,
"peak_pos": 0.05, "peak_hwhm": 0.01}, 0.07, 0.1725],
]
Loading