# 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) 2016-2017 Dominik Kriegner <dominik.kriegner@gmail.com>
import numbers
from math import sqrt
import numpy
from scipy import interpolate
from .. import utilities
from .powder import PowderDiffraction
from .smaterials import PowderList
[docs]class PowderModel(object):
"""
Class to help with powder calculations for multiple materials. For basic
calculations the Powder class together with the Fundamental parameters
approach is used.
"""
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.
Note:
After the end-of-use it is advisable to call the `close()` method to
cleanup the multiprocessing calculation!
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
In particular interesting in fpsettings might be:
{'displacement': {'specimen_displacement': z-displacement of the sample
from the rotation center
'zero_error_deg': zero error of the 2theta angle}
'absorption': {'sample_thickness': sample thickness (m),
'absorption_coefficient': sample's absorption (m^-1)}
'axial': {'length_sample': the length of the sample in the axial
direction (m)}
}
"""
if len(args) == 1 and isinstance(args[0], PowderList):
self.materials = args[0]
else:
self.materials = PowderList('%s List' % self.__class__.__name__,
*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
"""
lmfit = utilities.import_lmfit('XU.PowderModel')
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
"""
lmfit = utilities.import_lmfit('XU.PowderModel')
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(maxfev=maxfev)
self.set_lmfit_parameters(fitres.params)
return fitres
[docs] def close(self):
for pd in self.pdiff:
pd.close()
def __str__(self):
"""
string representation of the PowderModel
"""
ostr = "PowderModel {\n"
ostr += "I0: %f\n" % self.I0
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('Rp=%.4f Rwp=%.4f Rwpexp=%.4f chi2=%.4f'
% (Rp, Rwp, Rwpexp, chi2))
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='simulate', 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
simulated data
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
"""
plot, plt = utilities.import_matplotlib_pyplot('XU.simpack')
if not plot:
return
plt.figure(fig, figsize=(10, 7))
plt.clf()
ax = plt.subplot(111)
lines = []
if exp is not None:
lines.append(ax.plot(twotheta, exp, formatexp, label=labelexp)[0])
if mask is None:
mask = numpy.ones_like(twotheta, dtype=numpy.bool)
lines.append(ax.plot(twotheta[mask], sim, formatsim, label=labelsim)[0])
if show_diff:
# plot error between simulation and experiment
if exp is not None:
lines.append(ax.plot(twotheta[mask], exp[mask]-sim, '.-',
color='0.5', label='difference')[0])
plt.xlabel('2Theta (deg)')
plt.ylabel('Intensity')
plt.figlegend(lines, [l.get_label() for l in lines], loc='upper right',
frameon=True)
ax.set_yscale(scale)
plt.tight_layout()
return lines