Julia in Practice: Building Scattering.jl from Scratch (2)

Re-posted from: http://www.yxliu.group/2020/03/scattering-2

In this post we will implement a submodule, peak.jl, to model the shape of scattering peaks. Essential Julia language features will be introduced along the development of the submodule. You shall learn types, constructors, functions, methods, functors, modules, testing, and benchmarking after reading this post.

Shape of the Scattering Peak

For a perfect periodic lattice, i.e. a single crystal, as will be derived in later posts, the shape of the scattering peak has the form

$\begin{equation} L(x) = \frac{\sin^2(Nx/2)}{N^2\sin^2(x/2)} \end{equation}\label{eq:L-raw}$

where $N$ is number of unit cells which also determines the grain size of a single crystal. As the number of unit cells increases, more cells interfere constructively and the scattering peak becomes sharper. The trend has been shown in Figure 1. Note that we have shifted the peak position to $x=2$ since in scattering the peak position is always positive.

The shape of the scattering peak. — Figure 1. The shape of the scattering peak described by Eq.\eqref{eq:L-raw}.

Note that $N$ in the denominator of $L(x)$ is introduced to normalize the function. Although the function oscillates strongly beyond the central lobe, these oscillations are usually insignificant and can be safely ignored for common crystal grain size ($N$ > 50). In practice, however, it is more convenient to choose other non-oscillating peak shapes to simulate the actual one. Gaussian and Lorentzian peaks are two popular choices. The Gaussian peak is defined as

$\begin{equation} L(x) = \frac{2}{\pi\delta}\exp\left( -\frac{4x^2}{\pi\delta^2} \right) \end{equation}\label{eq:L-gauss}$

where the peak width (full-width at half-maximum, $h$) depends on the parameter $\delta$ as $h = \sqrt{\pi\ln2}\delta$. Note that we use a specialized form for the Gaussian distribution to be compatible with the following introduced generalized peak. The Lorentzian peak is

$\begin{equation} L(x) = \frac{\delta/2\pi}{x^2 + (\delta/2)^2} \end{equation}\label{eq:L-loren}$

where the peak width is simply $h = \delta$. A more generalized peak shape function has been proposed in the literature¹, which can be seen as a mixture of the Gaussian and Lorentzian peaks. It has the form

$\begin{equation} L(x) = \frac{2}{\pi\delta} \prod_{n=0}^{\infty} \left( 1 + \frac{\gamma_{\nu}^2}{(n+\nu/2)^2}\frac{4x^2}{\pi^2\delta^2} \right)^{-1} \end{equation}\label{eq:L-general}$

where the parameter is a ratio of gamma functions

$\begin{equation} \gamma_{\nu} = \sqrt{\pi}\frac{\Gamma[(\nu+1)/2]}{\Gamma(\nu/2)} \end{equation}\label{eq:gamma-nu}$

Its peak width is also $h = \delta$. The general peak shape has the property that it reduces to the Gaussian peak in the form of Eq.\eqref{eq:L-gauss} in the limit of $\nu \to \infty$, while it approaches to the Lorentzian peak in the form of Eq.\eqref{eq:L-loren} as $\nu \to 0$. Therefore, the parameter $\nu$ controls the shape of the generalized peak. This versatile peak shape allows one to account for the various contributions to peak broadening (instrumental resolution, grain size, grain shape, microstrain).

Implementation

Now we will dive into the implementation details. We will use this opportunity to introduce some core language features about Julia. You will see that it is quite neat to use Julia to do the work.

Types

The shape of the scattering peak is fully specified by two parameters, $\nu$ and $\delta$. When $\nu$ is bigger than a threshold value $\nu_{max}$, we use the Gaussian peak to approximate the generalized peak. And when $\nu$ is smaller than a threshold value $\nu_{min}$, we use the Lorentzian peak instead of the generalized peak. Therefore, we have to implement three types of peaks, which is natural to define a Julia type for each one. And they shall all inherit from an abstract type that unifies the interface to produce a specific peak shape.

In Julia, we can define an abstract type as

1 abstract type ScatteringPeak end

Now we can define a concrete Gaussian peak type inherited from it as

1
2
3
4



struct GaussianPeak <: ScatteringPeak
    σ::Real # the variance
    μ::Real # the mean value, i.e. the position of the peak
end

The first line states that the GaussianPeak type is inherited from the ScatteringPeak denoted by the symbol <:. Note that a concrete type such as GaussianPeak should be defined using a keyword struct other than type. The type has two fields, σ and μ. Here we plan to use the standard expression of the Gaussian peak:

$\begin{equation} L(x) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left[ -\frac{(x-\mu)^2}{2\sigma^2} \right] \end{equation}\label{eq:gauss}$

Thus, we have to do a conversion between eq.\eqref{eq:L-gauss} and eq.\eqref{eq:gauss} using $\delta = \sqrt{\frac{8}{\pi}}\sigma$ and $\mu=0$.

Similarly, we can define a Lorentzian type

1
2
3
4



struct LorentzianPeak <: ScatteringPeak
    δ::Real # the peak width
    μ::Real # the center of the peak
end

and a generalized peak type

1
2
3
4
5
6



struct GeneralizedPeak <: ScatteringPeak
    δ::Real # the variance
    ν::Real # the mean value, i.e. the position of the peak
    γν::Real # = √π * Γ[(ν+1)/2] / Γ(ν/2)
    γνhalf::Real # = Γ(ν/2)
end

Note that besides δ and ν, we have two more fields γν and γνhalf which store precomputed parameters in eq.\eqref{eq:gamma-nu} to save some time for later computations.

Constructors

Just as in C++, a constructor is a special function that use the type name as its name and create an instance of that type (instantiate an instance). Julia provides default constructors for all user-defined types, which accept the equal number of arguments as the number of its fields. Sometimes, however, it is convenient to custom our own constructors. Julia allows unlimited number of constructor for a type. For example, for the GaussianPeak type, the mean value parameter $\mu$ is commonly set to 0, so it is convenient to provide a constructor which takes only one arguments which provides the value of $\sigma$:

1
2
3
4



function GaussianPeak(σ::Real, μ::Real=0)
    @assert(σ > 0, "The variance must be positive.")
    new(σ, μ)
end

Now, we can initiate an instance as

1



p = GaussianPeak(0.1)

Using default constructor, we have to do the following

1



p = GaussianPeak(0.1, 0.0)

Our own constructor also checks the value of $\sigma$, ensuring its value is positive. We can write a similar constructor for the LorentzianPeak type.

Constructors like above are put outside of the definition of its corresponding type. They are outer constructors. When only outer constructors are presented for a type, the default constructor provided by Julia is still working. However, sometimes it is better to disable the default constructor when we don’t want user to provide values for some fields explicitly. For example, the fields γν and γνhalf both depend on another field ν. To ensure the consistency in the instance, we don’t want users of our package to initialize the GeneralizedPeak type using the default constructor. To this end, we can define an internal constructor which will override the default one:

1
2
3
4
5
6
7
8
9
10
11
12
13



struct GeneralizedPeak <: ScatteringPeak
    δ::Real # the variance
    ν::Real # the mean value, i.e. the position of the peak
    γν::Real # = √π * Γ[(ν+1)/2] / Γ(ν/2)
    γνhalf::Real # = Γ(ν/2)
    function GeneralizedPeak(δ::Real, ν::Real)
        @assert(δ > 0, "The peak width must be positive.")
        @assert(ν > 0, "The peak shape parameter must be positive.")
        γνhalf = gamma(ν/2)
        γν = √π * gamma((ν+1)/2) / γνhalf
        new(δ, ν, γν, γνhalf)
    end
end

With above definition, we can only initiate a GeneralizedPeak instance using

1



p = GeneralizedPeak(0.1, 1.0)

and following attempt will fail

1



p = GeneralizedPeak(0.1, 1.0, 0.0, 0.0)

Functions

It is useful to provide a unified interface to create a correct ScatteringPeak instance based on the values of $\delta$ and $\nu$. We can write a function for this purpose

1
2
3
4
5
6
7
8
9
10



function peak(δ::Real, ν::Real)
    if ν < 0.01
        LorentzianPeak(δ)
    elseif ν > 200
        σ = √(π/8) * δ
        GaussianPeak(σ)
    else
        GeneralizedPeak(δ, ν)
    end
end

It accepts two real numbers and returns a correct ScatteringPeak instance accordingly. Note that we have set $\nu_{min}=0.01$ and $\nu_{max}=200$. Also a conversion has been done for the GaussianPeak instance in line 5.

Tests

We can verify the peak in the Julia REPL

1
2
3



julia> p = peak(0.1, 1.0);
julia> typeof(p)
# GeneralizedPeak

There is also a very useful package, Test.jl, which provides functionalities to facilitate the testing work. We can perform a test as

1
2
3
4



julia> using Test
julia> p = peak(0.1, 1.0);
julia> @test typeof(p) isa GeneralizedPeak
# Test Passed

In the above, we first import the Test package use using keyword. Then @test macro is used to verify that the instance created by the peak function is a GeneralizedPeak peak as expected (since $\nu_{min} < \nu = 1.0 < \nu_{max}$).

Methods

Now it is time for us to add the essential part of each type: to actually create a peak. We shall use Julia methods. A method is a special function whose arguments are marked with explicit types. A function can have many methods. Therefore, a method can be viewed as a specific instance of a function.

1
2
3



peakshape(p::LorentzianPeak, x) = @. (p.δ/(2π)) / ((x-p.μ)^2 .+ (p.δ/2)^2)
peakshape(p::GaussianPeak, x) = @. exp(-(x-p.μ)^2/(2p.σ^2)) / (√(2π)*p.σ)
peakshape(p::GeneralizedPeak, x; tol::Real=1e-4, nmax::Integer=2000) = [compute_single_point(p, x[i], tol, nmax) for i in 1:length(x)]

Above codes define three methods all named peakshape for three types of scattering peak, respectively. Note that the keyword function can be omitted when the function body can be written in one line (one-liner). This greatly reduces the redundancy of the code and make the code looks extremely clean. The @. macro broadcast all operations to every elements of an array-like object. The implementation of compute_single_point will be presented in the following subsection.

Using these peakshape methods, we can now compute any type of scattering peaks, such as:

1
2
3
4
5
6



julia> p = peak(0.1, 1.0);
julia> q = collect(0.001:0.001:4.0);
julia> qhkl = 2.0;
julia> qs = q .- qhkl;
julia> peakshape(p, qs)
# produce a generalized peak with its peak locates at q = 2.0

The Julia compiler first deduces the type of instance p, which is GeneralizedPeak in the above example. Then it will dispatch it to the correct version of the peakshape method, which is the third one in the example.

We can compare these types of scattering peaks by plotting them together as shown in Figure 2.

Figure 2. A comparison of various scattering peaks.

It is obvious that the GeneralizedPeak is indeed a mixture of other two types.

Functors

Methods are associated with types, so it is possible to make any arbitrary Julia object “callable” by adding methods to its type. Such “callable” objects are sometimes called functors. Since the only functionality of the ScatteringPeak type is computing the peak shape, the functor fits this use case perfectly. For example, we can define functors as

1



(p::ScatteringPeak)(x) = peakshape(p, x)

The functor is a function without name and its only argument is a type instance. Now we can compute the peak simply using

1



julia> p(qs)

where p is a <:ScatteringPeak instance and qs is an array of $x$ values.

Benchmark and Performance

To compute a generalized peak, we have to evaluate an infinite product series in Eq.\eqref{eq:L-general}. By doing some numerical experiments, we learn that the series converges fast near and far away from the peak, but it converges much slower in between (around peak shoulder). So we have to allow the number of terms varying during computation. A tolerance is set to control the minimum number of terms for each $x$ value. The implementation is given below

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16



function compute_single_point(p::GeneralizedPeak, q::Real; tol::Real=1e-4, nmax::Integer=2000)
    c = 2/(p.δ*π)
    s = c
    s_prev = 0.0
    ds = 1.0
    n = 0
    while ds > tol && n < nmax
        s_prev = s
        t1 = (p.γν/(n+p.ν/2))^2
        t2 = (c * q)^2
        s /= (1 + t1*t2)
        ds = abs(s - s_prev)
        n += 1
    end
    s
end

This is the only computation intensive method in this submodule. To ensure its performance, we can run some benchmarks. Julia core provides a @time macro which can estimate the computation time of a method. However, its results are known to be disturbed by environment settings. So it is not quite reliable. The BenchmarkTools.jl package, on the other hand, provide a reliable macro @btime that isolates the environment. We can check out the performance of our code by running

1
2
3
4



julia> using BenchmarkTools
julia> p = peak(0.1, 1.0);
julia> @btime compute_single_point($p, 0.2)
# 11.xx μs (xx allocations: xx bytes)

Julia also provides a macro @code_warntype which helps us to identify the type instability issue. Using @code_warntype, it shows that our current implementation of compute_single_point is not optimized.

1
2



julia> @code_warntype compute_single_point(p, 0.2)
# pseudo codes presented here.

From the output of above @code_warntype, we can observe that Julia compiler cannot infer the types of c, s, s_prev, and ds due to the lack of type information of p.δ. To solve this issue, we should re-define GaussianPeak, LorentzianPeak and GeneralizedPeak to accept a type parameter which defines the type of their internal fields. For example,

1
2
3
4



struct GaussianPeak{T<:Real} <: ScatteringPeak
    σ::T # the variance
    μ::T # the mean value, i.e. the position of the peak
end

We should also update its constructor

1
2
3
4



function GaussianPeak(σ::T, μ::T=zero(T)) where {T<:Real}
    @assert(σ > 0, "The variance must be positive.")
    new{T}(σ, μ)
end

After all refactoring work is done, we re-run @code_warntype and will notice that all types are inferred by the compiler successfully. Using @btime to benchmark again, we will notice that the running time is only 207 ns which is 50 times faster!

Modules

The above codes are organized into a submodule and put in a single file named peak.jl

1 2 3 4 5 6 7 module Peak using SpecialFunctions # for gamma function # all other codes ... end

We use the gamma function in the SpecialFunctions.jl package to compute the $\Gamma$ function in Eq.\eqref{eq:gamma-nu}.

We will create a new file Scattering.jl to organize all its submodules and export Julia objects.

1 2 3 4 5 6 7 8 9 module ScatteringPeak include("peak.jl") using .Peak: ScatteringPeak, GaussianPeak, LorentzianPeak, GeneralizedPeak, peak export ScatteringPeak, GaussianPeak, LorentzianPeak, GeneralizedPeak export peak end

Currently, only the Peak submodule is included and exported.

Usage

Now we can generate a peak using the Scattering.jl package as follows

1
2
3
4



julia> using Scattering
julia> p = peak(0.1, 1.0);
julia> qs = collect(0.001:0.001:4.0) .- 2.0
julia> p(qs) # peak is generated here

Acknowledgements

This work is partially supported by the General Program of the National Natural Science Foundation of China (No. 21873021).

References

Yager, K. G.; Zhang, Y.; Lu, F.; Gang, O. Periodic Lattices of Arbitrary Nano-Objects: Modeling and Applications for Self-Assembled Systems. J. Appl. Crystallogr. 2013, 47, 118–129. ↩

Julia in Practice: Building Scattering.jl from Scratch (2) was originally published by Yi-Xin Liu at Yi-Xin Liu on March 23, 2020.

juliabloggers.com

A Julia Language Blog Aggregator