# This file is part of xrayutilities.
#
# xrayutilities is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
# Copyright (c) 2017-2023 Dominik Kriegner <dominik.kriegner@gmail.com>
import numbers
from math import sqrt
import lmfit
import numpy
from scipy import interpolate
from .. import utilities
from .powder import PowderDiffraction
from .smaterials import PowderList
[docs]
class PowderModel:
"""
Class to help with powder calculations for multiple materials. For basic
calculations the Powder class together with the Fundamental parameters
approach is used.
"""
[docs]
def __init__(self, *args, **kwargs):
"""
constructor for a powder model. The arguments consist of a PowderList
or individual Powder(s). Optional parameters are specified in the
keyword arguments.
Parameters
----------
args : PowderList or Powders
either one PowderList or several Powder objects can be given
kwargs : dict
optional parameters for the simulation. supported are:
fpclass : FP_profile, optional
derived class with possible convolver mixins. (default:
FP_profile)
fpsettings : dict
settings dictionaries for the convolvers. Default settings are
loaded from the config file.
I0 : float, optional
scaling factor for the simulation result
Notes
-----
- After the end-of-use it is advisable to call the `close()` method
to cleanup the multiprocessing calculation!
- In particular interesting keys in the fpsettings dictionary are
listed in the following. Note that this a short excerpt of the full
functionality:
- 'displacement'-dictionary with keys:
- 'specimen_displacement': sample's z-displacement from the
rotation center
- 'zero_error_deg': zero error of the 2theta angle
- 'absorption'-dictionary with keys:
- 'sample_thickness': sample thickness (m),
- 'absorption_coefficient': sample's absorption (m^-1)
- 'axial'-dictionary with keys:
- 'length_sample': sample length in the axial direction (m)
"""
if len(args) == 1 and isinstance(args[0], PowderList):
self.materials = args[0]
else:
self.materials = PowderList(f'{self.__class__.__name__} List',
*args)
self.I0 = kwargs.pop('I0', 1.0)
self.pdiff = []
kwargs['enable_simulation'] = True
for mat in self.materials:
self.pdiff.append(PowderDiffraction(mat, **kwargs))
# default background
self._bckg_type = 'polynomial'
self._bckg_pol = [0, ]
[docs]
def set_parameters(self, params):
"""
set simulation parameters of all subobjects
Parameters
----------
params : dict
settings dictionaries for the convolvers.
"""
# set parameters for each convolver
for pd in self.pdiff:
pd.update_settings(params)
[docs]
def set_lmfit_parameters(self, lmparams):
"""
function to update the settings of this class during an least squares
fit
Parameters
----------
lmparams : lmfit.Parameters
lmfit Parameters list of sample and instrument parameters
"""
pv = lmparams.valuesdict()
settings = dict()
h = list(self.pdiff[0].data)[0]
fp = self.pdiff[0].data[h]['conv'].convolvers
for conv in fp:
name = conv[5:]
settings[name] = dict()
self.I0 = pv.pop('primary_beam_intensity', 1)
set_splbkg = False
spliney = {}
for p in pv:
if p.startswith('phase_'): # sample phase parameters
midx = 0
for i, name in enumerate(self.materials.namelist):
if p.find(name) > 0:
midx = i
name = self.materials.namelist[midx]
attrname = p[p.find(name) + len(name) + 1:]
setattr(self.materials[midx], attrname, pv[p])
elif p.startswith('background_coeff'):
self._bckg_pol[int(p.split('_')[-1])] = pv[p]
elif p.startswith('background_spl_coeff'):
set_splbkg = True
spliney[int(p.split('_')[-1])] = pv[p]
else: # instrument parameters
for k in settings:
if p.startswith(k):
slist = p[len(k) + 1:].split('_')
if len(slist) > 2 and slist[-2] == 'item':
name = '_'.join(slist[:-2])
if slist[-1] == '0':
settings[k][name] = []
settings[k][name].append(pv[p])
else:
name = p[len(k) + 1:]
settings[k][name] = pv[p]
break
if set_splbkg:
self._bckg_spline = interpolate.InterpolatedUnivariateSpline(
self._bckg_spline._data[0],
[spliney[k] for k in sorted(spliney)], ext=0)
self.set_parameters(settings)
[docs]
def create_fitparameters(self):
"""
function to create a fit model with all instrument and sample
parameters.
Returns
-------
lmfit.Parameters
"""
params = lmfit.Parameters()
# sample phase parameters
for mat, name in zip(self.materials, self.materials.namelist):
for k in mat.__dict__:
attr = getattr(mat, k)
if isinstance(attr, numbers.Number):
params.add('_'.join(('phase', name, k)), value=attr,
vary=False)
# instrument parameters
settings = self.pdiff[0].settings
for pg in settings:
for p in settings[pg]:
val = settings[pg][p]
if p == 'dominant_wavelength' and pg == 'global':
# wavelength must be fit using emission_emiss_wavelength
continue
if isinstance(val, numbers.Number):
params.add('_'.join((pg, p)), value=val, vary=False)
elif isinstance(val, (numpy.ndarray, tuple, list)):
for j, item in enumerate(val):
params.add('_'.join((pg, p, 'item_%d' % j)),
value=item, vary=False)
# other global parameters
params.add('primary_beam_intensity', value=self.I0, vary=False)
if self._bckg_type == 'polynomial':
for i, coeff in enumerate(self._bckg_pol):
params.add('background_coeff_%d' % i, value=coeff, vary=False)
elif self._bckg_type == 'spline':
for i, coeff in enumerate(self._bckg_spline._data[1]):
params.add('background_spl_coeff_%d' % i, value=coeff,
vary=False)
return params
[docs]
def set_background(self, btype, **kwargs):
"""
define background as spline or polynomial function
Parameters
----------
btype : {polynomial', 'spline'}
background type; Depending on this
value the expected keyword arguments differ.
kwargs : dict
optional keyword arguments
x : array-like, optional
x-values (twotheta) of the background points (if btype='spline')
y : array-like, optional
intensity values of the background (if btype='spline')
p : array-like, optional
polynomial coefficients from the highest degree to the constant
term. len of p decides about the degree of the polynomial (if
btype='polynomial')
"""
if btype == 'spline':
self._bckg_spline = interpolate.InterpolatedUnivariateSpline(
kwargs.get('x'), kwargs.get('y'), ext=0)
elif btype == 'polynomial':
self._bckg_pol = list(kwargs.get('p'))
else:
raise ValueError("btype must be either 'spline' or 'polynomial'")
self._bckg_type = btype
[docs]
def simulate(self, twotheta, **kwargs):
"""
calculate the powder diffraction pattern of all materials and sum the
results based on the relative volume of the materials.
Parameters
----------
twotheta : array-like
positions at which the powder pattern should be evaluated
kwargs : dict
optional keyword arguments
background : array-like
an array of background values (same shape as twotheta) if no
background is given then the background is calculated as previously
set by the set_background function or is 0.
further keyword arguments are passed to the Convolve function of of the
PowderDiffraction objects
Returns
-------
array-like
summed powder diffraction intensity of all materials present in the
model
"""
inte = numpy.zeros_like(twotheta)
background = kwargs.pop('background', None)
if background is None:
if self._bckg_type == 'spline':
background = self._bckg_spline(twotheta)
else:
background = numpy.polyval(self._bckg_pol, twotheta)
totalvol = sum(pd.mat.volume for pd in self.pdiff)
for pd in self.pdiff:
inte += pd.Calculate(twotheta, **kwargs) * pd.mat.volume / totalvol
return self.I0 * inte + background
[docs]
def fit(self, params, twotheta, data, std=None, maxfev=200):
"""
make least squares fit with parameters supplied by the user
Parameters
----------
params : lmfit.Parameters
object with all parameters set as intended by the user
twotheta : array-like
angular values for the fit
data : array-like
experimental intensities for the fit
std : array-like
standard deviation of the experimental data. if 'None' the sqrt of
the data will be used
maxfev: int
maximal number of simulations during the least squares refinement
Returns
-------
lmfit.MinimizerResult
"""
def residual(pars, tt, data, weight):
"""
residual function for lmfit Minimizer routine
Parameters
----------
pars : lmfit.Parameters
fit Parameters
tt : array-like
array of twotheta angles
data : array-like
experimental data, same shape as tt
eps : array-like
experimental error bars, shape as tt
"""
# set parameters in this instance
self.set_lmfit_parameters(pars)
# run simulation
model = self.simulate(tt)
return (model - data) * weight
if std is None:
weight = numpy.reciprocal(numpy.sqrt(data))
else:
weight = numpy.reciprocal(std)
weight[numpy.isinf(weight)] = 1
self.minimizer = lmfit.Minimizer(residual, params,
fcn_args=(twotheta, data, weight))
fitres = self.minimizer.minimize(max_nfev=maxfev)
self.set_lmfit_parameters(fitres.params)
return fitres
[docs]
def plot(self, twotheta, showlines=True, label='simulation', color=None,
formatspec='-', lcolors=None, ax=None, **kwargs):
"""
plot the powder diffraction pattern and indicate line positions for all
components in the model.
Parameters
----------
twotheta : array-like
positions at which the powder pattern should be evaluated
showlines : bool, optional
flag to decide if peak positions of the components will be shown on
the top of the plot
label : str
line label in the plot
color : matplotlib color or None
the color used for the line plot of the simulation
formatspec : str
format specifier of the simulation curve
lcolors : list of matplotlib colors
colors for the line indicators for the various components
ax : matplotlib.axes or None
axes object to be used for plotting, if its given no axes
decoration like labels are set.
Further keyword arguments are passed to the simulate method.
Returns
-------
matplotlib.axes object or None if matplotlib is not available
"""
plot, plt = utilities.import_matplotlib_pyplot('XU.simpack')
if not plot:
return None
if ax is None:
fig, iax = plt.subplots()
iax.set_xlabel(r"$2\theta$ (°)")
iax.set_ylabel(r"Intensity")
iax.set_xlim(twotheta.min(), twotheta.max())
else:
fig = ax.figure
iax = ax
plotkwargs = dict(label=label)
if color is not None:
plotkwargs['color'] = color
iax.plot(twotheta, self.simulate(twotheta, **kwargs), formatspec,
**plotkwargs)
if showlines:
from matplotlib.colors import is_color_like
from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(iax)
taxlist = []
lineslist = []
annotlist = []
settings = self.pdiff[0].settings
wavelengths = settings['emission']['emiss_wavelengths']
intensities = settings['emission']['emiss_intensities']
for i, pd in enumerate(self.pdiff):
tax = divider.append_axes("top", size="6%", pad=0.05,
sharex=iax)
if lcolors:
c = lcolors[i % len(lcolors)]
elif len(self.pdiff) == 1:
if is_color_like(color):
c = color
elif is_color_like(formatspec[-1]):
c = formatspec[-1]
else:
c = 'C0'
else:
c = f'C{i}'
lw = 2
wllist = []
for wl, inte in zip(wavelengths, intensities):
q = [pd.data[h]['qpos'] for h in pd.data]
tt = pd.Q2Ang(q, wl=1e10*wl)*2
lw *= inte
lines = tax.vlines(tt, 0, 1, colors=c, linewidth=lw)
wllist.append(lines)
plt.setp(tax.get_xticklabels(), visible=False)
plt.setp(tax.get_yticklabels(), visible=False)
plt.setp(tax.get_yticklines(), visible=False)
annot = tax.annotate(
"", xy=(0, 0), xytext=(20, 0), textcoords="offset points",
bbox=dict(boxstyle="round,pad=0.1", fc="w", alpha=0.8),
arrowprops=dict(arrowstyle="->"), fontsize='x-small')
annot.set_visible(False)
# next line important to avoid zorder issues
tax.figure.texts.append(tax.texts[-1])
taxlist.append(tax)
lineslist.append(wllist)
annotlist.append(annot)
def update_annot(pd, annot, lines, ind):
h = list(pd.data)[ind]
x = 2*pd.data[h]['ang']
y = 0.5
annot.xy = (x, y)
text = f"{pd.mat.name}: {h[0]} {h[1]} {h[2]}"
annot.set_text(text)
annot.get_bbox_patch().set_edgecolor(lines.get_color()[0])
annot.set_zorder(10)
def hover(event):
for pd, tax, annot, wllist in zip(self.pdiff, taxlist,
annotlist, lineslist):
vis = annot.get_visible()
if event.inaxes == tax:
for lines in wllist:
cont, ind = lines.contains(event)
if cont:
update_annot(pd, annot, lines, ind['ind'][0])
annot.set_visible(True)
fig.canvas.draw_idle()
return
if vis:
annot.set_visible(False)
fig.canvas.draw_idle()
return
def click(event):
for pd, tax, _, wllist in zip(self.pdiff, taxlist,
annotlist, lineslist):
if event.inaxes == tax:
for lines in wllist:
cont, ind = lines.contains(event)
if cont:
h = list(pd.data)[ind['ind'][0]]
text = (f'{pd.mat.name}: {h[0]} {h[1]} {h[2]};'
f' 2Theta = ')
for wl in wavelengths:
tt = 2 * pd.Q2Ang(pd.data[h]['qpos'],
wl=1e10*wl)
text += f'{tt:.4f}°, '
print(text[:-2])
return
fig.canvas.mpl_connect("motion_notify_event", hover)
fig.canvas.mpl_connect("button_press_event", click)
if ax is None:
fig.tight_layout()
return iax
[docs]
def close(self):
for pd in self.pdiff:
pd.close()
def __str__(self):
"""
string representation of the PowderModel
"""
ostr = "PowderModel {\n"
ostr += f"I0: {self.I0:f}\n"
ostr += str(self.materials)
ostr += "}"
return ostr
[docs]
def Rietveld_error_metrics(exp, sim, weight=None, std=None,
Nvar=0, disp=False):
"""
calculates common error metrics for Rietveld refinement.
Parameters
----------
exp : array-like
experimental datapoints
sim : array-like
simulated data
weight : array-like, optional
weight factor in the least squares sum. If it is None the weight is
estimated from the counting statistics of 'exp'
std : array-like, optional
standard deviation of the experimental data. alternative way of
specifying the weight factor. when both are given weight overwrites
std!
Nvar : int, optional
number of variables in the refinement
disp : bool, optional
flag to tell if a line with the calculated values should be printed.
Returns
-------
M, Rp, Rwp, Rwpexp, chi2: float
"""
if weight is None and std is None:
weight = numpy.reciprocal(exp)
elif weight is None:
weight = numpy.reciprocal(std**2)
weight[numpy.isinf(weight)] = 1
M = numpy.sum((exp - sim)**2 * weight)
Rp = numpy.sum(numpy.abs(exp - sim))/numpy.sum(exp)
Rwp = sqrt(M / numpy.sum(weight * exp**2))
chi2 = M / (len(exp) - Nvar)
Rwpexp = Rwp / sqrt(chi2)
if disp:
print(f"Rp={Rp:.4f} Rwp={Rwp:.4f} Rwpexp={Rwpexp:.4f} chi2={chi2:.4f}")
return M, Rp, Rwp, Rwpexp, chi2
[docs]
def plot_powder(twotheta, exp, sim, mask=None, scale='sqrt', fig='XU:powder',
show_diff=True, show_legend=True, labelexp='experiment',
labelsim='simulation', formatexp='.-k', formatsim='-r'):
"""
Convenience function to plot the comparison between experimental and
simulated powder diffraction data
Parameters
----------
twotheta : array-like
angle values used for the x-axis of the plot (deg)
exp : array-like
experimental data (same shape as twotheta). If None only the simulation
and no difference will be plotted
sim : array-like or PowderModel
simulated data or PowderModel instance. If a PowderModel instance is
given the plot-method of PowderModel is used.
mask : array-like, optional
mask to reduce the twotheta values to the be used as x-coordinates of
sim
scale : {'linear', 'sqrt', 'log'}, optional
string specifying the scale of the y-axis.
fig : str or int, optional
matplotlib figure name (figure will be cleared!)
show_diff : bool, optional
flag to specify if a difference curve should be shown
show_legend: bool, optional
flag to specify if a legend should be shown
labelexp : str
plot label (legend entry) for the experimental data
labelsim : str
plot label for the simulation data
formatexp : str
format specifier for the experimental data
formatsim : str
format specifier for the simulation curve
Returns
-------
List of lines in the plot. Empty list in case matplotlib can't be imported
"""
plot, plt = utilities.import_matplotlib_pyplot('XU.simpack')
if not plot:
return []
if scale == 'sqrt':
from ..mpl_helper import SqrtAllowNegScale # noqa: F401
f = plt.figure(fig, figsize=(10, 7))
f.clf()
ax = plt.subplot(111)
if exp is not None:
ax.plot(twotheta, exp, formatexp, label=labelexp)
if mask is None:
mask = numpy.ones_like(twotheta, dtype=bool)
if isinstance(sim, PowderModel):
simdata = sim.simulate(twotheta[mask])
sim.plot(twotheta[mask], label=labelsim, formatspec=formatsim, ax=ax)
else:
simdata = sim
ax.plot(twotheta[mask], simdata, formatsim, label=labelsim)
if show_diff and exp is not None:
# plot error between simulation and experiment
ax.plot(twotheta[mask], exp[mask]-simdata, '.-', color='0.5',
label='difference')
ax.set_xlabel('2Theta (deg)')
ax.set_ylabel('Intensity')
lines = ax.get_lines()
if show_legend:
plt.figlegend(lines, [line.get_label() for line in lines],
loc='upper right', frameon=True)
ax.set_yscale(scale)
f.tight_layout()
return lines