Split Pseudo-Voigt Modelling with Drop-In Modules (WIP)

[JRLarkin]

Hi authors, I'm attempting to accurately model experimental XRD line-profiles and although pseudo-Voigt and Pearson7 functions have typically been used in the past, these are less able to approximate profiles arising from heavily radiation-damaged samples.
The asymmetry in experimental diffraction peaks can be best approximated by a split Pseudo-Voigt model (although I am also working with more sophisticated models, i.e. CMWP).
I have written drop-in functions for the lineshapes.py and models.py modules that attempt a split-PV fit, based heavily on the logic in your own SplitLorentzian model, but I am struggling with ensuring continuity at the peak centre (I(0) left = I(0) right). As you can see, the final fits have a 'shelf' at the top of the peak which implies that the height parameter is not conserved across the two pseudo-Voigt functions.

Here is my drop-in function for the lineshapes.py module:

def split_pvoigt(x, amplitude=1.0, center=0.0, sigma=1.0, sigma_r=1.0, fraction=0.5, fraction_r=0.5):

    """
    split_pvoigt(x, amplitude, center, sigma, sigma_r, fraction, fraction_r) =
    pv(i_1) + pv(i_2) where i_1 = f(x, amp, center, sigma, frac) and i_2 = f(x, amp, center, sigma_r, frac_r)

    Definitions for split pseudo-Voigt functions and I_hkl(I_1, I_2) from:
    T. Delgado, C. Enachescu, A. Tissot, A. Hauser, L. Guénée and C. Besnard, Evidencing size-dependent cooperative effects on spin 
    crossover nanoparticles following their HS - LS relaxation, J. Mater. Chem. C, 2018, 6, 12698–12706.

    """

    amp = amplitude
    frac = fraction
    frac_r = fraction_r

    g = sqrt((log2/pi))
    sigma_g = sigma / sqrt(2 * log2)
    sigma_g_r = sigma_r / sqrt(2 * log2)
    h = max(0, 1)

    i_denom = ((frac*2)/(pi*sigma) + (1-frac)*(2*g/sigma)) + ((frac_r*2)/(pi*sigma_r) + (1-frac_r)*(2*g/sigma_r))

    i_1 = amp * ((frac_r*(2/(pi*sigma_r))) + ((1-frac_r)*(2*g/sigma_r))) / i_denom

    i_2 = amp * ((frac*(2/pi*sigma)) + ((1-frac)*(2*g/sigma))) / i_denom

    # final pv fit is either pv(sigma, frac) OR pv(sigma_r, frac_r).

    pv_1 = (h * (x < center)) * ((1 - frac) * gaussian(x, i_1, center, sigma_g) + frac * lorentzian(x, i_1, center, sigma))

    pv_2 = (h * (x >= center)) * ((1 - frac_r) * gaussian(x, i_2, center, sigma_g_r) + frac_r * lorentzian(x, i_2, center, sigma_r))

    return pv_1 + pv_2

And here is my drop-in class for the models.py module:

class SplitPseudoVoigtModel(Model):
    r"""A model based on a pseudo-Voigt distribution function.

    The left-hand PV function has the parameters (x, A, center, sigma, fraction), and
    the right-hand PV function replaces sigma & fraction with sigma_r & fraction_r.

    The intensity of the left-hand PV function is multiplied by zero if the condition (x < center) is not met.
    Likewise for the right-hand PV function [(x >= center) is False].

    fwhm = sigma + sigma_r
    """

    fwhm_factor = 2.0

    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
                 **kwargs):
        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                       'independent_vars': independent_vars})
        super().__init__(split_pvoigt, **kwargs)
        self._set_paramhints_prefix()

    def _set_paramhints_prefix(self):
        fwhm_expr = '{pre:s}sigma+{pre:s}sigma_r'
        height_expr = '2*{pre:s}amplitude/{0:.7f}/max({1:.7f}, ({pre:s}sigma+{pre:s}sigma_r))'
        self.set_param_hint('sigma', min=0)
        self.set_param_hint('sigma_r', min=0)
        self.set_param_hint('fraction', value=0.5, min=0.0, max=1.0)
        self.set_param_hint('fraction_r', value=0.5, min=0.0, max=1.0)
        self.set_param_hint('fwhm', expr=fwhm_expr.format(pre=self.prefix))
        fmt = ("(((1-{prefix:s}fraction)*{prefix:s}amplitude)/"
               "max({0}, ({prefix:s}sigma*sqrt(pi/log(2))))+"
               "({prefix:s}fraction*{prefix:s}amplitude)/"
               "max({0}, (pi*{prefix:s}sigma)))")

        # self.set_param_hint('height', expr=fmt.format(tiny, prefix=self.prefix))
        self.set_param_hint('height', expr=height_expr.format(np.pi, tiny, pre=self.prefix))

    def guess(self, data, x, negative=False, **kwargs):
        """Estimate initial model parameter values from data."""
        pars = guess_from_peak(self, data, x, negative, ampscale=1.25)

        sigma = pars[f'{self.prefix}sigma']
        pars[f'{self.prefix}sigma_r'].set(value=sigma.value, min=sigma.min, max=sigma.max)

        fraction = pars[f'{self.prefix}fraction']
        pars[f'{self.prefix}fraction_r'].set(value=fraction.value, min=fraction.min, max=fraction.max)

        return update_param_vals(pars, self.prefix, **kwargs)

    __init__.__doc__ = COMMON_INIT_DOC
    guess.__doc__ = COMMON_GUESS_DOC

Here's a plot of the model versus an experimental diffraction peak:

Close-up of the continuity issue at the top of the peak:

Admittedly, I have reached the limit of my LMFIT knowledge! Would you be able to help me solve this continuity problem?

Kind regards,
-Jake