# 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-2021, 2023 Dominik Kriegner <dominik.kriegner@gmail.com>
import numpy
from lmfit import Model
from .. import utilities
from . import models
[docs]
class FitModel(Model):
"""
Wrapper for the lmfit Model class working for instances of LayerModel
Typically this means that after initialization of `FitModel` you want to
use make_params to get a `lmfit.Parameters` list which one customizes for
fitting.
Later on you can call `fit` and `eval` methods with those parameter list.
"""
[docs]
def __init__(self, lmodel, verbose=False, plot=False, elog=True, **kwargs):
"""
initialization of a FitModel which uses lmfit for the actual fitting,
and generates an according lmfit.Model internally for the given
pre-configured LayerModel, or subclasses thereof which includes models
for reflectivity, kinematic and dynamic diffraction.
Parameters
----------
lmodel : LayerModel
pre-configured instance of LayerModel or any subclass
verbose : bool, optional
flag to enable verbose printing during fitting
plot : bool or str, optional
flag to decide wheter an plot should be created showing the fit's
progress. If plot is a string it will be used as figure name, which
makes reusing the figures easier.
elog : bool, optional
flag to enable a logarithmic error metric between model and data.
Since the dynamic range of data is often large its often benefitial
to keep this enabled.
kwargs : dict, optional
additional keyword arguments are forwarded to the `simulate` method
of `lmodel`
"""
self.verbose = verbose
self.plot = plot
self.elog = elog
assert isinstance(lmodel, models.LayerModel)
self.lmodel = lmodel
# generate dynamic function for model evalution
funcstr = "def func(x, "
# add LayerModel parameters
for p in self.lmodel.fit_paramnames:
funcstr += f"{p}, "
# add LayerStack parameters
for layer in self.lmodel.lstack:
for param in self.lmodel.lstack_params:
funcstr += f"{layer.name}_{param}, "
if self.lmodel.lstack_structural_params:
for param in layer._structural_params:
funcstr += f"{layer.name}_{param}, "
funcstr += "lmodel=self.lmodel, **kwargs):\n"
# define modelfunc content
for p in self.lmodel.fit_paramnames:
funcstr += f" setattr(lmodel, '{p}', {p})\n"
for i, layer in enumerate(self.lmodel.lstack):
for param in self.lmodel.lstack_params:
varname = f"{layer.name}_{param}"
cmd = " setattr(lmodel.lstack[{}], '{}', {})\n"
funcstr += cmd.format(i, param, varname)
if self.lmodel.lstack_structural_params:
for param in layer._structural_params:
varname = f"{layer.name}_{param}"
cmd = " setattr(lmodel.lstack[{}], '{}', {})\n"
funcstr += cmd.format(i, param, varname)
# perform actual model calculation
funcstr += " return lmodel.simulate(x, **kwargs)"
namespace = {"self": self}
exec(funcstr, {"lmodel": self.lmodel}, namespace)
self.func = namespace["func"]
self.emetricfunc = numpy.log10 if self.elog else lambda x: x
def _residual(params, data, weights, **kwargs):
"""
Return the residual. This is a (simplified, only real values)
reimplementation of the lmfit.Model._residual function which adds
the possibility of a logarithmic error metric.
Default residual: (data-model)*weights.
"""
scale = self.emetricfunc
model = scale(self.eval(params, **kwargs))
sdata = scale(data)
mask = numpy.logical_and(
numpy.isfinite(model), numpy.isfinite(sdata)
)
diff = model[mask] - sdata[mask]
if weights is not None and scale(1) == 1:
diff *= weights
return numpy.asarray(diff).ravel()
super().__init__(
self.func,
independent_vars="x",
name=self.lmodel.__class__.__name__,
**kwargs,
)
self._residual = _residual
# set default parameter hints
self._default_hints()
self.set_fit_limits()
[docs]
def set_fit_limits(self, xmin=-numpy.inf, xmax=numpy.inf, mask=None):
"""
set fit limits. If mask is given it must have the same size as the
`data` and `x` variables given to fit. If mask is None it will be
generated from xmin and xmax.
Parameters
----------
xmin : float, optional
minimum value of x-values to include in the fit
xmax : float, optional
maximum value of x-values to include in the fit
mask : boolean array, optional
mask to be used for the data given to the fit
"""
self.mask = mask
self.xmin = xmin
self.xmax = xmax
[docs]
def fit(
self,
data,
params,
x,
weights=None,
fit_kws=None,
lmfit_kws=None,
**kwargs,
):
"""
wrapper around lmfit.Model.fit which enables plotting during the
fitting
Parameters
----------
data : ndarray
experimental values
params : lmfit.Parameters
list of parameters for the fit, use make_params for generation
x : ndarray
independent variable (incidence angle or q-position depending on
the model)
weights : ndarray, optional
values of weights for the fit, same size as data
fit_kws : dict, optional
Options to pass to the minimizer being used
lmfit_kws : dict, optional
keyword arguments which are passed to lmfit.Model.fit
kwargs : dict, optional
keyword arguments passed to lmfit.Model.eval
Returns
-------
lmfit.ModelResult
"""
if not lmfit_kws:
lmfit_kws = {}
class FitPlot:
def __init__(self, figname, logscale):
self.figname = figname
self.logscale = logscale
if not self.figname:
self.plot = False
else:
f, plt = utilities.import_matplotlib_pyplot("XU.simpack")
self.plt = plt
self.plot = f
def plot_init(self, x, data, weights, model, mask, verbose):
if not self.plot:
return
self.plt.ion()
if isinstance(self.figname, str):
self.fig = self.plt.figure(self.figname)
else:
self.fig = self.plt.figure("XU:FitModel")
self.plt.clf()
self.ax = self.plt.subplot(111)
if weights is not None:
eline = self.ax.errorbar(
x,
data,
yerr=1 / weights,
ecolor="0.3",
fmt="ok",
errorevery=int(x.size / 80),
label="data",
)[0]
else:
(eline,) = self.ax.plot(x, data, "ok", label="data")
if verbose:
self.ax.plot(x, model, "-", color="0.5", label="initial")
if eline:
self.zord = eline.zorder + 2
else:
self.zord = 1
if self.logscale:
self.ax.set_yscale("log")
self.fline = None
def showplot(self, xlab, ylab="Intensity (arb. u.)"):
if not self.plot:
return
self.plt.xlabel(xlab)
self.plt.ylabel(ylab)
self.plt.legend()
self.fig.set_tight_layout(True)
self.plt.show()
def updatemodelline(self, x, newmodel):
if not self.plot:
return
try:
self.plt.sca(self.ax)
except ValueError:
return
if self.fline is None:
(self.fline,) = self.ax.plot(
x, newmodel, "-r", lw=2, label="fit", zorder=self.zord
)
else:
self.fline.set_data(x, newmodel)
canvas = self.fig.canvas # see plt.draw function (avoid show!)
canvas.draw_idle()
canvas.start_event_loop(0.001)
def addfullmodelline(self, x, y):
if not self.plot:
return
self.ax.plot(
x, y, "-g", lw=1, label="full model", zorder=self.zord - 1
)
if self.mask:
mask = self.mask
else:
mask = numpy.logical_and(x >= self.xmin, x <= self.xmax)
mweights = weights
if mweights is not None:
mweights = weights[mask]
# create initial plot
self.fitplot = FitPlot(self.plot, self.elog)
initmodel = self.eval(params, x=x, **kwargs)
self.fitplot.plot_init(x, data, weights, initmodel, mask, self.verbose)
self.fitplot.showplot(xlab=self.lmodel.xlabelstr)
# create callback function
def cb_func(params, niter, resid, *args, **kwargs):
if self.verbose:
print(f"{niter:04d} {numpy.sum(resid**2):12.3e}")
if self.fitplot.plot and niter % 20 == 0:
self.fitplot.updatemodelline(
kwargs["x"], self.eval(params, **kwargs)
)
# perform fitting
lmfit_kws.update(kwargs) # add kwargs for function call
res = super().fit(
data[mask],
params,
x=x[mask],
weights=mweights,
fit_kws=fit_kws,
iter_cb=cb_func,
**lmfit_kws,
)
# final update of plot
if self.fitplot.plot:
try:
self.fitplot.plt.sca(self.fitplot.ax)
except ValueError:
self.fitplot.plot_init(
x, data, weights, initmodel, mask, self.verbose
)
fittedmodel = self.eval(res.params, x=x, **kwargs)
self.fitplot.addfullmodelline(x, fittedmodel)
self.fitplot.updatemodelline(x[mask], fittedmodel[mask])
self.fitplot.showplot(xlab=self.lmodel.xlabelstr)
return res
def _default_hints(self):
"""
set useful hints for parameters all LayerModels have
"""
# general parameters
for pn in self.lmodel.fit_paramnames:
self.set_param_hint(pn, value=getattr(self.lmodel, pn), vary=False)
for pn in ("I0", "background"):
self.set_param_hint(pn, vary=True, min=0)
self.set_param_hint("resolution_width", min=0, vary=False)
self.set_param_hint("energy", min=1000, vary=False)
# parameters of the layerstack
for lay in self.lmodel.lstack:
for param in self.lmodel.lstack_params:
varname = f"{lay.name}_{param}"
self.set_param_hint(varname, value=getattr(lay, param), min=0)
if param == "density":
self.set_param_hint(
varname, max=1.5 * lay.material.density
)
if param == "thickness":
self.set_param_hint(varname, max=2 * lay.thickness)
if param == "roughness":
self.set_param_hint(varname, max=50)
if self.lmodel.lstack_structural_params:
for param in lay._structural_params:
varname = f"{lay.name}_{param}"
self.set_param_hint(
varname, value=getattr(lay, param), vary=False
)
if "occupation" in param:
self.set_param_hint(varname, min=0, max=1)
if "biso" in param:
self.set_param_hint(varname, min=0, max=5)
if self.lmodel.lstack[0].thickness == numpy.inf:
varname = f"{self.lmodel.lstack[0].name}_thickness"
self.set_param_hint(varname, vary=False)