Peak profile convolution

[Jonas Beyer]

Hello all,

I am fairly new to diffpy and have a question regarding convolutions to the peak profile. What I am looking for is the possibility to convolve inherent gaussian peaks (either with peak widths from DebyeWallerPeakWidth or JeongPeakWidth) with a custom peak profile function that I define myself. So far, I have been able to find a work-around by defining a custom PeakProfile() object. My code looks roughly likes this:

--------------------------------------------------------------------------

class CustomPeakProfile(PeakProfile):

def __call__(self, x, fwhm):
        xrange = np.arange(-3*fwhm, 3*fwhm, rstep)
        gauss = 2/fwhm*np.sqrt(np.log(2)/np.pi)*np.exp(-4*np.log(2)*xrange**2/fwhm**2)
        

        custom function = <custom function with same length as gauss>

                   
        conv_prod_ = np.convolve(gauss, custom, 'same')
        conv_prod = conv_prod_/abs(np.trapz(conv_prod_, x=xrange))
       
        xidx = np.abs(xrange - x).argmin()
     
        return conv_prod[xidx]

def clone(self):
        import copy
        return copy.copy(self)

def create(self):
        return XiaProfile()

def type(self):
        return "custom"

def xboundhi(self, fwhm):
        return 3*fwhm

def xboundlo(self, fwhm):
        return -3*fwhm

--------------------------------------------------------------------------

The problem with this is that it is rather stupid: It calculates the convolution product for every "x" for every peak. Even for small structures, it is very slow.

Is there any way to convolute the peaks in a more efficient manner? Maybe even some kind of functionality which allows for this?

All the best,

Jonas Beyer

Postdoc at Aarhus University

[Simon Billinge]

Thanks for the question Jonas.  I am not sure how to do this, does anyone know?  But I would think a friendly code design would be if we could sequentially apply convolutions within the calculator.  I am not sure if that is possible...

--
You received this message because you are subscribed to the Google Groups "diffpy-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to diffpy-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/diffpy-users/80c95e12-b21f-4319-9173-3f467fbcb186n%40googlegroups.com.

--

Simon Billinge

Professor, Columbia University

Physicist, Brookhaven National Laboratory

[Peter Metz]

A clarifying question: by "It calculates the convolution product for every "x" for every peak," do you mean that you are convoluting your custom distribution in direct space?

If so, I'm not sure I understand the rationale. On the other hand if you're applying the convolution theorem to introduce aberrations in Q-space, you should only need one pass.

In either case, can you speed up your calculation by masking your distribution below a certain threshold? i.e. many of your points are probably within a floating point tolerance of zero.

Happy computing,
Peter

[Jonas Beyer]

To clarify, I want to convolute every PDF peak with the custom function in direct space, and the custom function actually changes as a function of r according to some parameters.
Because of the latter, the final PDF I want to compute (after the conv. prod. has been calculated for every peak) is strictly speaking not a convolution, but rather a 'broadening integral'. As such, the convolution theorem does not apply.

I want to access every peak in the PDF as a list of r-coordinates contributing to the peak and the corresponding intensity.

However, The "__call__" method receives the "x" (i.e. r-coordinate) of a specific peak one-by-one, which is then repeated for all peaks.

Best,

Jonas

[Peter Metz]

Hi Jonas, 

I'm guessing you found your way to the example here? https://github.com/diffpy/cmi_exchange/blob/master/cmi_scripts/pdfrectprofile/rectangleprofile.py 

I had forgotten that in a previous life I had considered a similar question-- it looks like it didn't get much traction: https://groups.google.com/g/diffpy-users/c/HGPlHXtJjJQ/m/YRwr3JixHAAJ 

I think there are two things to consider:

(1) DiffPy is boosted (i.e. a python wrapping to C++ code), so the built in peak generation is going to be faster. It seems the default profile generated is a truncated gaussian (https://github.com/diffpy/libdiffpy/blob/master/src/diffpy/srreal/PDFCalculator.cpp) I notice you already specified cutoffs in your custom class, so my earlier advice was probably not useful.

(2)  Possibly you don't actually need convolution for the way the PDF is synthesized. If your peak width model is length dependent, returning an x-dependent peak width in from the `__call__` method of your custom class might be sufficient? This might be a fundamental roadblock to implementation in DiffPy, on the other hand-- you would essentially be writing functionals, which may be incompatible with the PDFCalculator design.

That said, I abandoned this project in DiffPy (sorry devs!) and deserted to Topas.

Here are two examples attempting a range dependent peak width using DiffPy, either modifying (subclassing) or replacing the standard model:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

"""
Example use:

from diffpy_helpers import CustomPeakWidth, CustomConstantPeakWidth

calc = rec.zif8.G._calc
calc.peakwidthmodel = CustomConstantPeakWidth()
calc.peakwidthmodel._deregisterType('jeong')
calc.peakwidthmodel.width = 0.5
calc.rcut = 6.5
calc.c1 = 0.03
print(calc.peakwidthmodel)
evaluate(rec)
ax.plot(rec.zif8.G.profile.x, rec.zif8.G.profile.ycalc, 'k--')

Created on Sun Mar  1 10:45:26 2020

@author: peter
"""
from copy import copy

from diffpy.srreal.peakprofile import PeakProfile
from diffpy.srreal.peakwidthmodel import PeakWidthModel, JeongPeakWidth, DebyeWallerPeakWidth, ConstantPeakWidth

class CustomConstantPeakWidth(ConstantPeakWidth):
    """ """
    def model(self, pt):
        """
        model(pt) multiplicatively scales the Jeong peak width model
       
        must be overriden in custom class instance.
       
        Args:
            - pt (float)
       
        Returns:
            model(pt) (float) | 1.
        """
        if not hasattr(self, 'rcut'):
            self.rcut = 0.
        if not hasattr(self, 'c1'):
            self.c1 = 0.
        if pt <= self.rcut:
            return 1.
        elif pt > self.rcut:
            return 1 + (pt - self.rcut) * self.c1
       
    def calculate(self, bnds) -> float:
        """" """
        return self.width * self.model(bnds.distance())

# End CustomPeakWidth

class CustomPeakWidth(JeongPeakWidth):
    """ """
       
    def model(self, pt):
        """
        model(pt) multiplicatively scales the Jeong peak width model
       
        must be overriden in custom class instance.
       
        Args:
            - pt (float)
       
        Returns:
            model(pt) (float) | 1.
        """
        if not hasattr(self, 'rcut'):
            self.rcut = 0.
        if not hasattr(self, 'c1'):
            self.c1 = 0.
        if pt <= self.rcut:
            return 1.
        elif pt > self.rcut:
            return 1 + (pt - self.rcut) * self.c1

    def type(self):
        """ """
        return 'custompeakwidth'
       
    def calculate(self, bnds) -> float:
        """" """
        return super().calculate(bnds) *  self.model(bnds.distance())

# End CustomPeakWidth

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

And here's a similar strategy implemented in Topas

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

/*
Developing a range-dependent / correlated motion model for molecular crystals
 
Author: Peter Metz
Date: August 18 2020 11:18 EST
 
*/
 
' --- func defs ---
fn jeong(sigma, del1, del2, Qb)
   {
   ' Jeong peak width model a.k.a. beqpdfgui
   def default = 0;
   ' handle div by zero
   def arg = If(X <= 0, default, (1 - (del1)/X - (del2)/X^2 + (Qb)^2 X^2)) ;
   ' handle imaginary
   def eva = If(arg <= 0, default, sigma Sqrt(arg)) ;
   ' return
   return eva ;
   }
   
fn shifted_exp(rcut, sigma, exp)
   {
   ' sigma(r) = sigma ( 1 + (X-rcut) ^ exp )
   return sigma If( X <= rcut, 1, 1 + (X-rcut)^exp ) ;
   }
 
fn mod1(sigma, del1, del2, Qb, rcut, exp)
   {
   ' functional sigma(r) = shifted( jeong(r) )
   return shifted_exp(rcut, jeong(sigma, del1, del2, Qb), exp) ;
   }
   
fn mod2(sigma, del1, del2, Qb, rcut, exp)
   {
   ' product sigma(r) = jeong(r) shifted_exp(r)
   return jeong(sigma, del1, del2, Qb) shifted_exp(rcut, 1, exp) ;
   }

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

If you're unacquainted with Topas' syntax multiplication is implicit, so `mod1` is a functional, whereas `mod2` is a product.

I hope some of this is helpful, and attracts more discussion around this very interesting problem! I see it as a fundamental pathway forward for materials presenting multiple dissimilar bonding modes.

-Peter