#!/usr/bin/env python
#
# Copyright (C) 2015, the ixpeobssim team.
#
# This program 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 3 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, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
from __future__ import print_function, division
import numpy
import scipy.stats
from ixpeobssim.utils.math_ import fold_angle_deg, fold_angle_rad
from ixpeobssim.utils.matplotlib_ import plt, xStatBox
from ixpeobssim.utils.logging_ import abort
[docs]
class xFitModelBase:
"""Base class for a fittable model.
This base class isn't really doing anything useful, the idea being that
actual models that can be instantiated subclass xFitModelBase overloading
the relevant class members.
The class features a number of static members that derived class
should redefine as needed:
* ``PARAMETER_NAMES`` is a list containing the names of the model
parameters. (It goes without saying thet its length should match the
number of parameters in the model.)
* ``PARAMETER_DEFAULT_VALUES`` is a list containing the default values
of the model parameters. When a concrete model object is instantiated
these are the values being attached to it at creation time.
* ``PARAMETER_DEFAULT_BOUNDS`` is a tuple containing the default values
of the parameter bounds to be used for the fitting. The values in the
tuple are attached to each model object at creation time and are
intended to be passed as the ``bounds`` argument of the
``scipy.optimize.curve_fit()`` function. From the ``scipy`` documentation:
Lower and upper bounds on independent variables. Defaults to no bounds.
Each element of the tuple must be either an array with the length equal
to the number of parameters, or a scalar (in which case the bound is
taken to be the same for all parameters.) Use np.inf with an appropriate
sign to disable bounds on all or some parameters. By default
models have no built-in bounds.
* ``DEFAULT_PLOTTING_RANGE`` is a two-element list with the default support
(x-axis range) for the model. This is automatically updated at runtime
depending on the input data when the model is used in a fit.
* ``DEFAULT_STAT_BOX_POSITION`` is the default location of the stat box for
the model, see the :meth:`gpdswpy.matplotlib_.stat_box()` function for
all the details.
In addition, each derived class should override the following things:
* the ``value(x, *args)`` static method: this should return the value of
the model at a given point for a given set of values of the underlying
parameters;
* (optionally) the ``jacobian(x, *args)`` static method. (If defined, this
is passed to the underlying fit engine allowing to reduce the number of
function calls in the fit; otherwise the jacobian is calculated
numerically.)
Finally, if there is a sensible way to initialize the model parameters
based on a set of input data, derived classes should overload the
``init_parameters(xdata, ydata, sigma)`` method of the base class, as the
latter is called by fitting routines if no explicit array of initial values
are passed as an argument. The default behavior of the class method defined
in the base class is to do nothing.
See :class:`gpdswpy.modeling.xGaussian` for a working example.
"""
PARAMETER_NAMES = ()
PARAMETER_DEFAULT_VALUES = ()
PARAMETER_DEFAULT_BOUNDS = (-numpy.inf, numpy.inf)
DEFAULT_PLOTTING_RANGE = (0., 1.)
DEFAULT_STAT_BOX_POSITION = 'upper right'
def __init__(self):
"""Constructor.
Here we initialize the class members holding the best-fit parameter
values and the associated covariance matrix, see the
:pymeth:`gpdswpy.modeling.xFitModelBase.reset()` method.
We also create a (private) look-up table holding the correspondence
between the parameter names and the corresponding position in the
parameter list and we cache the default support for the model for
plotting purposes.
"""
assert len(self.PARAMETER_NAMES) == len(self.PARAMETER_DEFAULT_VALUES)
self.__parameter_dict = {}
for i, name in enumerate(self.PARAMETER_NAMES):
self.__parameter_dict[name] = i
self.reset()
[docs]
def name(self):
"""Return the model name.
"""
return self.__class__.__name__
[docs]
def reset(self):
"""Reset all the necessary stuff.
This method initializes all the things that are necessry to keep
track of a parametric fit.
* the parameter values are set to what it specified in
``PARAMETER_DEFAULT_VALUES``;
* the covatiance matrix is initialized to a matrix of the proper
dimension filled with zeroes;
* the minimum and maximum values of the independent variable
(relevant for plotting) are set to the values specified in
``DEFAULT_PLOTTING_RANGE``;
* the model bounds are set to the values specified in
``PARAMETER_DEFAULT_BOUNDS``;
* the chisquare and number of degrees of freedom are initialized to
-1.
"""
self.parameters = numpy.array(self.PARAMETER_DEFAULT_VALUES, dtype='d')
self.covariance_matrix = numpy.zeros((len(self), len(self)), dtype='d')
self.xmin, self.xmax = self.DEFAULT_PLOTTING_RANGE
self.bounds = self.PARAMETER_DEFAULT_BOUNDS
self.chisq = -1.
self.ndof = -1
[docs]
def reduced_chisquare(self):
"""Return the reduced chisquare.
"""
if self.ndof > 0:
return self.chisq / self.ndof
return -1.
def __getattr__(self, name):
"""Short-hand method to retrieve the parameter values by name.
Note that we manipulate the attribute name by capitalizing the
first letter and replacing underscores with spaces.
"""
name = name.title().replace('_', ' ')
if name in self.PARAMETER_NAMES:
return self.parameter_value(name)
else:
raise AttributeError
[docs]
@staticmethod
def value(x, *parameters):
"""Eval the model at a given point and a given set of parameter values.
Warning
-------
This needs to be overloaded in any derived class for the thing to do
something sensible.
"""
raise 'value() not implemented'
[docs]
def integral(self, edges):
"""Calculate the integral of the model within pre-defined edges.
Note that this assumes that the derived class provides a suitable
``cdf()`` method.
"""
try:
return self.cdf(edges[1:], *self.parameters) - \
self.cdf(edges[:-1], *self.parameters)
except Exception as e:
abort('%s.integral() not implemened (%s)' % (self.name, e))
[docs]
def rvs(self, size=1):
"""Return random variates from the model.
Note that this assumes that the derived class provides a suitable
``ppf()`` method.
"""
try:
return self.ppf(numpy.random.random(size), *self.parameters)
except Exception as e:
abort('%s.rvs() not implemened (%s)' % (self.name(), e))
def __call__(self, x, *parameters):
"""Return the value of the model at a given point and a given set of
parameter values.
Note that unless the proper number of parameters is passed to the
function call, the model is evaluated at the best-fit parameter values.
The function is defined with this signature because it is called
with a set of parameter values during the fit process, while
tipically we want to evaluate it with the current set of parameter
values after the fact.
"""
if len(parameters) == len(self):
return self.value(x, *parameters)
else:
return self.value(x, *self.parameters)
def __parameter_index(self, name):
"""Convenience method returning the index within the parameter vector
for a given parameter name.
"""
assert(name in self.PARAMETER_NAMES)
return self.__parameter_dict[name]
[docs]
def parameter_value(self, name):
"""Return the parameter value by name.
"""
index = self.__parameter_index(name)
return self.parameters[index]
[docs]
def parameter_error(self, name):
"""Return the parameter error by name.
"""
index = self.__parameter_index(name)
return numpy.sqrt(self.covariance_matrix[index][index])
[docs]
def parameter_values(self):
"""Return the vector of parameter values.
"""
return self.parameters
[docs]
def parameter_errors(self):
"""Return the vector of parameter errors.
"""
return numpy.sqrt(self.covariance_matrix.diagonal())
[docs]
def parameter_status(self):
"""Return the complete status of the model in the form of a tuple
of tuples (parameter_name, parameter_value, parameter_error).
Note this can be overloaded by derived classes if more information
needs to be added.
"""
return tuple(zip(self.PARAMETER_NAMES, self.parameter_values(),
self.parameter_errors()))
[docs]
def set_parameters(self, *pars):
"""Set all the parameter values.
Note that the arguments must be passed in the right order.
"""
self.parameters = numpy.array(pars, dtype='d')
[docs]
def set_parameter(self, name, value):
"""Set a parameter value.
"""
index = self.__parameter_index(name)
self.parameters[index] = value
[docs]
def init_parameters(self, xdata, ydata, sigma):
"""Assign a sensible set of values to the model parameters, based
on a data set to be fitted.
Note that in the base class the method is not doing anything, but it
can be reimplemented in derived classes to help make sure the
fit converges without too much manual intervention.
"""
pass
[docs]
def set_plotting_range(self, xmin, xmax):
"""Set the plotting range.
"""
self.xmin = xmin
self.xmax = xmax
[docs]
def plot(self, *parameters, **kwargs):
"""Plot the model.
Note that when this is called with a full set of parameters, the
self.parameters class member is overwritten so that the right values
can then be picked up if the stat box is plotted.
"""
if len(parameters) == len(self):
self.parameters = parameters
display_stat_box = kwargs.pop('display_stat_box', False)
x = numpy.linspace(self.xmin, self.xmax, 1000)
y = self(x, *parameters)
plt.plot(x, y, **kwargs)
if display_stat_box:
self.stat_box(**kwargs)
[docs]
def stat_box(self, position=None, plot=True, **kwargs):
"""Plot a ROOT-style stat box for the model.
"""
if position is None:
position = self.DEFAULT_STAT_BOX_POSITION
box = xStatBox(position)
box.add_entry('Fit model: %s' % self.name())
box.add_entry('Chisquare', '%.1f / %d' % (self.chisq, self.ndof))
for name, value, error in self.parameter_status():
box.add_entry(name, value, error)
if plot:
box.plot(**kwargs)
return box
def __len__(self):
"""Return the number of model parameters.
"""
return len(self.PARAMETER_NAMES)
def __add__(self, other):
"""Add two models.
Warning
-------
This is highly experimental and guaranteed to contain bugs. Enjoy.
"""
m1 = self
m2 = other
xmin = min(m1.DEFAULT_PLOTTING_RANGE[0], m2.DEFAULT_PLOTTING_RANGE[0])
xmax = max(m1.DEFAULT_PLOTTING_RANGE[1], m2.DEFAULT_PLOTTING_RANGE[1])
name = '%s + %s' % (m1.__class__.__name__, m2.__class__.__name__)
# In order to correctly propagate the parameter boundaries of the two
# input models to their sum, we need to handle the case where one or
# both of them do not define their PARAMETER_DEFAULT_BOUNDS, relying
# on the base class default. In that case we simply repeat those
# default values for as many parameters as there are in that model.
lower_bounds = []
upper_bounds = []
for m in (m1, m2):
if m.bounds == xFitModelBase.PARAMETER_DEFAULT_BOUNDS:
lower_bounds += (len(m.parameters) * [m.bounds[0]])
upper_bounds += (len(m.parameters) * [m.bounds[1]])
else:
lower_bounds += m.bounds[0]
upper_bounds += m.bounds[1]
class _model(xFitModelBase):
PARAMETER_NAMES = m1.PARAMETER_NAMES + m2.PARAMETER_NAMES
PARAMETER_DEFAULT_VALUES = m1.PARAMETER_DEFAULT_VALUES + \
m2.PARAMETER_DEFAULT_VALUES
DEFAULT_PLOTTING_RANGE = (xmin, xmax)
PARAMETER_DEFAULT_BOUNDS = (lower_bounds, upper_bounds)
def __init__(self):
self.__class__.__name__ = name
xFitModelBase.__init__(self)
@staticmethod
def value(x, *parameters):
return m1.value(x, *parameters[:len(m1)]) +\
m2.value(x, *parameters[len(m1):])
@staticmethod
def jacobian(x, *parameters):
return numpy.hstack((m1.jacobian(x, *parameters[:len(m1)]),
m2.jacobian(x, *parameters[len(m1):])))
return _model()
def __str__(self):
"""String formatting.
"""
text = '%s model (chisq/ndof = %.2f / %d)' % (self.__class__.__name__,
self.chisq, self.ndof)
for name, value, error in self.parameter_status():
text += '\n%15s: %.5e +- %.5e' % (name, value, error)
return text
[docs]
class xConstant(xFitModelBase):
"""Constant model.
.. math::
f(x; C) = C
"""
PARAMETER_NAMES = ('Constant',)
PARAMETER_DEFAULT_VALUES = (1.,)
DEFAULT_PLOTTING_RANGE = (0., 1.)
[docs]
@staticmethod
def value(x, constant):
"""Overloaded value() method.
"""
return numpy.full(x.shape, constant)
[docs]
@staticmethod
def jacobian(x, constant):
"""Overloaded jacobian() method.
"""
d_constant = numpy.full((len(x),), 1.)
return numpy.array([d_constant]).transpose()
[docs]
def cdf(self, x):
"""Overloaded cdf() method.
"""
return self.Constant * x
[docs]
def ppf(self, q):
"""Overloaded ppf() method.
"""
return self.xmin + q * (self.xmax - self.xmin)
[docs]
def init_parameters(self, xdata, ydata, sigma):
"""Overloaded init_parameters() method.
"""
self.set_parameter('Constant', numpy.mean(ydata))
[docs]
class xLine(xFitModelBase):
"""Straight-line model.
.. math::
f(x; m, q) = mx + q
"""
PARAMETER_NAMES = ('Intercept', 'Slope')
PARAMETER_DEFAULT_VALUES = (1., 1.)
DEFAULT_PLOTTING_RANGE = (0., 1.)
[docs]
@staticmethod
def value(x, intercept, slope):
"""Overloaded value() method.
"""
return intercept + slope * x
[docs]
@staticmethod
def jacobian(x, intercept, slope):
"""Overloaded jacobian() method.
"""
d_intercept = numpy.full((len(x),), 1.)
d_slope = x
return numpy.array([d_intercept, d_slope]).transpose()
[docs]
class xGaussian(xFitModelBase):
"""One-dimensional Gaussian model.
.. math::
f(x; A, \\mu, \\sigma) = A e^{-\\frac{(x - \\mu)^2}{2\\sigma^2}}
"""
PARAMETER_NAMES = ('Amplitude', 'Peak', 'Sigma')
PARAMETER_DEFAULT_VALUES = (1., 0., 1.)
PARAMETER_DEFAULT_BOUNDS = ([0., -numpy.inf, 0], [numpy.inf] * 3)
DEFAULT_PLOTTING_RANGE = (-5., 5.)
SIGMA_TO_FWHM = 2.3548200450309493
[docs]
@staticmethod
def value(x, amplitude, peak, sigma):
"""Overloaded value() method.
"""
return amplitude * numpy.exp(-0.5 * ((x - peak)**2. / sigma**2.))
[docs]
@staticmethod
def der_amplitude(x, amplitude, peak, sigma):
"""Return the amplitude derivative of the function, to be used in the
calculation of the Jacobian.
"""
return numpy.exp(-0.5 / sigma**2. * (x - peak)**2.)
[docs]
@staticmethod
def der_peak(x, amplitude, d_amplitude, peak, sigma):
"""Return the peak derivative of the function, to be used in the
calculation of the Jacobian.
Note that we pass the pre-calculated values of the amplitude derivatives
in order not to repeat the same calculation more times than strictly
necessary.
"""
return amplitude * d_amplitude * (x - peak) / sigma**2.
[docs]
@staticmethod
def der_sigma(x, amplitude, d_amplitude, peak, sigma):
"""Return the sigma derivative of the function, to be used in the
calculation of the Jacobian.
Note that we pass the pre-calculated values of the amplitude derivatives
in order not to repeat the same calculation more times than strictly
necessary.
"""
return amplitude * d_amplitude * (x - peak)**2. / sigma**3.
[docs]
@staticmethod
def jacobian(x, amplitude, peak, sigma):
"""Overloaded jacobian() method.
"""
d_amplitude = xGaussian.der_amplitude(x, amplitude, peak, sigma)
d_peak = xGaussian.der_peak(x, amplitude, d_amplitude, peak, sigma)
d_sigma = xGaussian.der_sigma(x, amplitude, d_amplitude, peak, sigma)
return numpy.array([d_amplitude, d_peak, d_sigma]).transpose()
[docs]
def init_parameters(self, xdata, ydata, sigma):
"""Overloaded init_parameters() method.
"""
self.set_parameter('Amplitude', numpy.max(ydata))
self.set_parameter('Peak', numpy.mean(xdata))
self.set_parameter('Sigma', numpy.std(xdata))
[docs]
def fwhm(self):
"""Return the absolute FWHM of the model.
"""
return self.SIGMA_TO_FWHM * self.sigma
[docs]
def resolution(self):
"""Return the resolution of the model, i.e., the FWHM divided by the
peak value.
"""
if self.peak > 0:
return self.fwhm() / self.peak
return 0.
[docs]
def resolution_error(self):
"""Return the error on the resolution.
"""
if self.peak > 0 and self.parameter_error('Sigma') > 0:
return self.resolution() * self.parameter_error('Sigma') /\
self.parameter_value('Sigma')
return 0.
[docs]
class xFe55(xGaussian):
"""One-dimensional double gaussian model for Fe55 with 2 lines
.. math::
f(x; A, \\mu, \\sigma) = A e^{-\\frac{(x - \\mu)^2}{2\\sigma^2}} + cose
"""
PARAMETER_NAMES = ('Amplitude0', 'Amplitude1', 'Peak', 'Sigma')
PARAMETER_DEFAULT_VALUES = (1., 0.25, 0., 1.)
PARAMETER_DEFAULT_BOUNDS = ([0., 0., -numpy.inf, 0.], [numpy.inf] * 4)
DEFAULT_PLOTTING_RANGE = (-5., 5.)
PEAK_SCALE = 6.490 / 5.895
SIGMA_SCALE = numpy.sqrt(PEAK_SCALE)
@staticmethod
def _main_gaussian(x, amplitude0, peak, sigma):
"""
"""
return xGaussian.value(x, amplitude0, peak, sigma)
@staticmethod
def _secondary_gaussian(x, amplitude1, peak, sigma):
"""
"""
return xGaussian.value(x, amplitude1, xFe55.PEAK_SCALE * peak,
xFe55.SIGMA_SCALE * sigma)
[docs]
@staticmethod
def value(x, amplitude0, amplitude1, peak, sigma):
"""Overloaded value() method.
"""
return xGaussian.value(x, amplitude0, peak, sigma) + \
xGaussian.value(x, amplitude1, xFe55.PEAK_SCALE * peak,
xFe55.SIGMA_SCALE * sigma)
[docs]
@staticmethod
def jacobian(x, amplitude0, amplitude1, peak, sigma):
"""Overloaded jacobian() method.
"""
peak1 = xFe55.PEAK_SCALE * peak
sigma1 = xFe55.SIGMA_SCALE * sigma
d_amplitude0 = xGaussian.der_amplitude(x, amplitude0, peak, sigma)
d_amplitude1 = xGaussian.der_amplitude(x, amplitude0, peak1, sigma1)
d_peak = xGaussian.der_peak(x, amplitude0, d_amplitude0, peak,
sigma) + \
xGaussian.der_peak(x, amplitude1, d_amplitude1, peak1,
sigma1)
d_sigma = xGaussian.der_sigma(x, amplitude0, d_amplitude0, peak,
sigma) +\
xGaussian.der_sigma(x, amplitude1, d_amplitude1, peak1,
sigma1)
return numpy.array([d_amplitude0, d_amplitude1,
d_peak, d_sigma]).transpose()
[docs]
def init_parameters(self, xdata, ydata, sigma):
"""
"""
pass
[docs]
def plot(self, *parameters, **kwargs):
"""Overloaded plot() method.
"""
amplitude0, amplitude1, peak, sigma = self.parameters
xFitModelBase.plot(self, *parameters, **kwargs)
x = numpy.linspace(self.xmin, self.xmax, 1000)
y = self(x, *self.parameters)
y = xGaussian.value(x, amplitude0, peak, sigma)
plt.plot(x, y, ls='dashed')
peak1 = xFe55.PEAK_SCALE * peak
sigma1 = xFe55.SIGMA_SCALE * sigma
y = xGaussian.value(x, amplitude1, peak1, sigma1)
plt.plot(x, y, ls='dashed')
[docs]
class xModulationCurveRad(xFitModelBase):
"""Modulation curve model (for fitting azimuthal distributions expressed
in radians).
.. math::
f(x; A, m, \\bar\\varphi) =
A (1 + m \\cos\\left(2(x - \\bar\\varphi)\\right))
"""
PARAMETER_NAMES = ('Amplitude', 'Modulation', 'Phase')
PARAMETER_DEFAULT_VALUES = (1., 1., 0.)
PARAMETER_DEFAULT_BOUNDS = ([0., 0., -0.5 * numpy.pi],
[numpy.inf, 1., 0.5 * numpy.pi])
DEFAULT_PLOTTING_RANGE = (-numpy.pi, numpy.pi)
DEFAULT_STAT_BOX_POSITION = 'lower right'
[docs]
@staticmethod
def value(x, amplitude, modulation, phase):
"""Overloaded value() method.
"""
return amplitude * (1 + modulation * numpy.cos(2 * (x - phase)))
[docs]
@staticmethod
def jacobian(x, amplitude, modulation, phase):
"""Overloaded jacobian() method.
"""
d_amplitude = xModulationCurveRad.value(x, 1., modulation, phase)
d_modulation = amplitude * (numpy.cos(2 * (x - phase)))
d_phase = -2. * amplitude * modulation * numpy.sin(2 * (phase - x))
return numpy.array([d_amplitude, d_modulation, d_phase]).transpose()
[docs]
def init_parameters(self, xdata, ydata, sigma):
"""Overloaded init_parameters() method.
"""
self.set_parameter('Amplitude', numpy.mean(ydata))
ymin = numpy.min(ydata)
ymax = numpy.max(ydata)
self.set_parameter('Modulation', (ymax - ymin) / (ymax + ymin))
self.set_parameter('Phase', fold_angle_rad(xdata[numpy.argmax(ydata)]))
[docs]
class xModulationCurveDeg(xFitModelBase):
"""Modulation curve model (for fitting azimuthal distributions expressed
in degrees).
.. math::
f(x; A, m, \\bar\\varphi) =
A (1 + m \\cos\\left(2(x - \\bar\\varphi)\\right))
"""
PARAMETER_NAMES = ('Amplitude', 'Modulation', 'Phase')
PARAMETER_DEFAULT_VALUES = (1., 1., 0.)
PARAMETER_DEFAULT_BOUNDS = ([0., 0., -90.], [numpy.inf, 1., 90.])
DEFAULT_PLOTTING_RANGE = (-180., 180.)
DEFAULT_STAT_BOX_POSITION = 'lower right'
DEG_TO_RAD = numpy.pi / 180.
[docs]
@staticmethod
def value(x, amplitude, modulation, phase):
"""Overloaded value() method.
Here we are essentially calling the modulation curve model expressed
in radians converting to radians the input values of the independent
variable and the phase parameter.
"""
return xModulationCurveRad.value(numpy.radians(x), amplitude,
modulation, numpy.radians(phase))
[docs]
@staticmethod
def jacobian(x, amplitude, modulation, phase):
"""Overloaded jacobian() method.
In addition to the degree-to-radians conversion of the inputs, here we
have to multiply the last column of the jacobian (i.e., the phase
derivatives) by pi / 180.
"""
_jac = xModulationCurveRad.jacobian(numpy.radians(x), amplitude,
modulation, numpy.radians(phase))
_jac[:,2] = xModulationCurveDeg.DEG_TO_RAD * _jac[:,2]
return _jac
[docs]
def init_parameters(self, xdata, ydata, sigma):
"""Overloaded init_parameters() method.
"""
self.set_parameter('Amplitude', numpy.mean(ydata))
ymin = numpy.min(ydata)
ymax = numpy.max(ydata)
self.set_parameter('Modulation', (ymax - ymin) / (ymax + ymin))
self.set_parameter('Phase', fold_angle_deg(xdata[numpy.argmax(ydata)]))
[docs]
class xPowerLaw(xFitModelBase):
"""Power law model.
.. math::
f(x; N, \\Gamma) = N x^\\Gamma
"""
PARAMETER_NAMES = ('Normalization', 'Index')
PARAMETER_DEFAULT_VALUES = (1., -1.)
PARAMETER_DEFAULT_BOUNDS = ([0., -numpy.inf], [numpy.inf, numpy.inf])
DEFAULT_PLOTTING_RANGE = (1.e-2, 1.)
[docs]
@staticmethod
def value(x, normalization, index):
"""Overloaded value() method.
"""
return normalization * (x**index)
[docs]
@staticmethod
def jacobian(x, normalization, index):
"""Overloaded jacobian() method.
"""
d_normalization = (x**index)
d_index = numpy.log(x) * normalization * (x**index)
return numpy.array([d_normalization, d_index]).transpose()
[docs]
class xPowerLawExpCutoff(xFitModelBase):
"""Power law with an exponential cutoff.
.. math::
f(x; N, \\Gamma, x_0) = N x^\\Gamma e^{-\\frac{x}{x_0}}
"""
PARAMETER_NAMES = ('Normalization', 'Index', 'Cutoff')
PARAMETER_DEFAULT_VALUES = (1., -1., 1.)
PARAMETER_DEFAULT_BOUNDS = ([0., -numpy.inf, 0.], [numpy.inf] * 3)
DEFAULT_PLOTTING_RANGE = (1.e-2, 1.)
[docs]
@staticmethod
def value(x, normalization, index, cutoff):
"""Overloaded value() method.
"""
return normalization * (x**index) * numpy.exp(-x / cutoff)
[docs]
@staticmethod
def jacobian(x, normalization, index, cutoff):
"""Overloaded jacobian() method.
"""
d_normalization = (x**index) * numpy.exp(-x / cutoff)
val = xPowerLawExpCutoff.value(x, normalization, index, cutoff)
d_index = numpy.log(x) * val
d_cutoff = val * x / (cutoff**2)
return numpy.array([d_normalization, d_index, d_cutoff]).transpose()
[docs]
class xPixelPha(xFitModelBase):
"""Composite model for fitting the single-pixel pulse-height distributions.
This is essentially the sum of a gaussian centered at zero (and modeling
the noise) and a power law with exponential cutoff (modeling the signal).
.. math::
f(x; A_n, \\sigma_n, N_s, \\Gamma_s, x_0) =
A_n e^{-\\frac{x^2}{2\\sigma_n^2}} +
N_s x^\\Gamma_s e^{-\\frac{x}{x_0}}
"""
PARAMETER_NAMES = ('Noise amplitude', 'Noise sigma',
'Signal normalization', 'Signal index', 'Signal cutoff')
PARAMETER_DEFAULT_VALUES = (1., 1., 1., -1., 1.)
DEFAULT_PLOTTING_RANGE = (0.5, 10.)
@staticmethod
def _noise(x, amplitude, sigma):
"""Noise component.
"""
return xGaussian.value(x, amplitude, 0., sigma)
@staticmethod
def _signal(x, normalization, index, cutoff):
"""Signal component.
"""
return xPowerLawExpCutoff.value(x, normalization, index, cutoff)
[docs]
@staticmethod
def value(x, noise_amplitude, noise_sigma, signal_normalization,
signal_index, signal_cutoff):
"""Overloaded value() method.
"""
return xPixelPha._noise(x, noise_amplitude, noise_sigma) + \
xPixelPha._signal(x, signal_normalization, signal_index,
signal_cutoff)
[docs]
def plot(self, *parameters, **kwargs):
"""Overloaded plot() method.
"""
xFitModelBase.plot(self, *parameters, **kwargs)
x = numpy.linspace(self.xmin, self.xmax, 1000)
y = self(x, *self.parameters)
ymin = y.min()
ymax = y.max()
y = xPixelPha._noise(x, *self.parameters[:2])
_mask = (y >= ymin) * (y <= ymax)
plt.plot(x[_mask], y[_mask], ls='dashed')
y = xPixelPha._signal(x, *self.parameters[2:])
_mask = (y >= ymin) * (y <= ymax)
plt.plot(x[_mask], y[_mask], ls='dashed')
[docs]
@staticmethod
def jacobian(x, noise_amplitude, noise_sigma, signal_normalization,
signal_index, signal_cutoff):
"""Overloaded jacobian() method.
"""
d_noise_amplitude = xGaussian.value(x, 1., 0., noise_sigma)
d_noise_sigma = noise_amplitude * (x)**2. * d_noise_amplitude /\
(noise_sigma**3.)
d_signal_normalization = (x)**signal_index *\
numpy.exp(-x / signal_cutoff)
val = xPowerLawExpCutoff.value(x, signal_normalization, signal_index,
signal_cutoff)
d_signal_index = numpy.log(x) * val
d_signal_cutoff = val * x / (signal_cutoff**2.)
return numpy.array([d_noise_amplitude, d_noise_sigma,
d_signal_normalization, d_signal_index,
d_signal_cutoff]).transpose()
[docs]
class xExponential(xFitModelBase):
"""Exponential model.
.. math::
f(x; N, \\alpha) = N e^{\\alpha x}
"""
PARAMETER_NAMES = ('Normalization', 'Index')
PARAMETER_DEFAULT_VALUES = (1., -1.)
PARAMETER_DEFAULT_BOUNDS = ([0., -numpy.inf], [numpy.inf]*2)
DEFAULT_PLOTTING_RANGE = (0., 1.)
[docs]
@staticmethod
def value(x, normalization, exponent):
"""Overloaded value() method.
"""
return normalization * numpy.exp(exponent * x)
[docs]
@staticmethod
def jacobian(x, normalization, exponent):
"""Overloaded jacobian() method.
"""
d_normalization = numpy.exp(exponent * x)
d_exponent = normalization * x * numpy.exp(exponent * x)
return numpy.array([d_normalization, d_exponent]).transpose()
[docs]
class xExponentialOffset(xFitModelBase):
"""Exponential model.
.. math::
f(x; A, N, \\alpha) = A + N e^{\\alpha x}
"""
PARAMETER_NAMES = ('Offset', 'Normalization', 'Index')
PARAMETER_DEFAULT_VALUES = (0., 1., -1.)
PARAMETER_DEFAULT_BOUNDS = ([-numpy.inf, 0., -numpy.inf], [numpy.inf]*3)
DEFAULT_PLOTTING_RANGE = (0., 1.)
[docs]
@staticmethod
def value(x, offset, normalization, exponent):
"""Overloaded value() method.
"""
return offset + normalization * numpy.exp(exponent * x)
[docs]
@staticmethod
def jacobian(x, offset, normalization, exponent):
"""Overloaded jacobian() method.
"""
d_offset = numpy.full((len(x),), 1.)
d_normalization = numpy.exp(exponent * x)
d_exponent = normalization * x * numpy.exp(exponent * x)
return numpy.array([d_offset, d_normalization, d_exponent]).transpose()
[docs]
class xLogNormal(xFitModelBase):
"""Log-normal fitting model.
See
https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.lognorm.html
for more details on the implementation.
Note that our parametrization is done in terms of the mean and sigma of the
distribution, in addition to the shape parameter.
"""
PARAMETER_NAMES = ('Normalization', 'Mean', 'Sigma', 'Shape')
PARAMETER_DEFAULT_VALUES = (1., 0., 1., 0.02)
DEFAULT_PLOTTING_RANGE = (-5., 5.)
[docs]
@staticmethod
def value(x, norm, mu, sigma, shape):
"""Overloaded value() method.
"""
p = numpy.exp(shape * shape)
scale = sigma / numpy.sqrt(p * (p - 1.))
loc = mu - scale * numpy.sqrt(p)
return norm * scipy.stats.lognorm.pdf(x, shape, loc, scale)
[docs]
class xGeneralizedGaussian(xFitModelBase):
"""Generalized gaussian fitting model.
See
https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gennorm.html
for implementation details.
"""
PARAMETER_NAMES = ('Normalization', 'Mean', 'Sigma', 'Beta')
PARAMETER_DEFAULT_VALUES = (1., 0., 1., 1.)
DEFAULT_PLOTTING_RANGE = (-5., 5.)
[docs]
@staticmethod
def value(x, norm, mean, sigma, beta):
"""Overloaded value() method.
"""
return norm * scipy.stats.gennorm.pdf(x, beta, mean, sigma)
[docs]
class xHat(xFitModelBase):
"""Fit model representing a flat distribution truncated with an error
function at both ends.
The ASYMMETRY class member controls the slope of the central part of the model.
"""
PARAMETER_NAMES = ('Normalization', 'X1', 'Sigma1', 'X2', 'Sigma2')
PARAMETER_DEFAULT_VALUES = (1., 0.2, 0.05, 0.8, 0.05)
DEFAULT_PLOTTING_RANGE = (0., 1.)
ASYMMETRY = 0.
[docs]
@staticmethod
def value(x, norm, x1, sigma1, x2, sigma2):
"""Overloaded value() method.
"""
z1 = (x - x1) / sigma1
z2 = (x - x2) / sigma2
dx = (x2 - x1)
val = norm * 0.25 * (1. + scipy.special.erf(z1)) * (1. + scipy.special.erf(-z2)) / dx
if xHat.ASYMMETRY != 0:
val *= 1. - xHat.ASYMMETRY * (x - x1) / dx
return val
[docs]
class xLorentzian(xFitModelBase):
"""Lorentzian function
"""
PARAMETER_NAMES = ('Normalization', 'Peak', 'HWHM' )
PARAMETER_DEFAULT_VALUES = (1., 0., 1.)
PARAMETER_DEFAULT_BOUNDS = ([0., -numpy.inf, 0.], [numpy.inf, numpy.inf, numpy.inf])
DEFAULT_PLOTTING_RANGE = (-1., 1.)
[docs]
@staticmethod
def value(x, normalization, peak, hwhm):
"""Overloaded value() method.
"""
return normalization * hwhm / ((x - peak)**2. + hwhm**2.)
[docs]
@staticmethod
def jacobian(x, normalization, peak, hwhm):
"""Overloaded jacobian() method.
"""
d_x = x - peak
denom = d_x**2. + hwhm**2.
d_normalization = hwhm / denom
d_peak = 2. * normalization * hwhm * d_x / denom**2.
d_hwhm = normalization * (d_x**2. - hwhm**2.)/ denom**2.
return numpy.array([d_normalization, d_peak, d_hwhm]).transpose()