# /*##########################################################################
#
# Copyright (c) 2004-2023 European Synchrotron Radiation Facility
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
########################################################################### */
"""This modules provides a set of fit functions and associated
estimation functions in a format that can be imported into a
:class:`silx.math.fit.FitManager` instance.
These functions are well suited for fitting multiple gaussian shaped peaks
typically found in spectroscopy data. The estimation functions are designed
to detect how many peaks are present in the data, and provide an initial
estimate for their height, their center location and their full-width
at half maximum (fwhm).
The limitation of these estimation algorithms is that only gaussians having a
similar fwhm can be detected by the peak search algorithm.
This *search fwhm* can be defined by the user, if
he knows the characteristics of his data, or can be automatically estimated
based on the fwhm of the largest peak in the data.
The source code of this module can serve as template for defining your own
fit functions.
The functions to be imported by :meth:`FitManager.loadtheories` are defined by
a dictionary :const:`THEORY`: with the following structure::
from silx.math.fit.fittheory import FitTheory
THEORY = {
'theory_name_1': FitTheory(
description='Description of theory 1',
function=fitfunction1,
parameters=('param name 1', 'param name 2', …),
estimate=estimation_function1,
configure=configuration_function1,
derivative=derivative_function1),
'theory_name_2': FitTheory(…),
}
.. note::
The order of the provided dictionary is taken into account.
Theory names can be customized (e.g. ``gauss, lorentz, splitgauss``…).
The mandatory parameters for :class:`FitTheory` are ``function`` and
``parameters``.
You can also define an ``INIT`` function that will be executed by
:meth:`FitManager.loadtheories`.
See the documentation of :class:`silx.math.fit.fittheory.FitTheory`
for more information.
Module members:
---------------
"""
import numpy
import logging
from silx.math.fit import functions
from silx.math.fit.peaks import peak_search, guess_fwhm
from silx.math.fit.filters import strip, savitsky_golay
from silx.math.fit.leastsq import leastsq
from silx.math.fit.fittheory import FitTheory
_logger = logging.getLogger(__name__)
__authors__ = ["V.A. Sole", "P. Knobel"]
__license__ = "MIT"
__date__ = "15/05/2017"
DEFAULT_CONFIG = {
"NoConstraintsFlag": False,
"PositiveFwhmFlag": True,
"PositiveHeightAreaFlag": True,
"SameFwhmFlag": False,
"QuotedPositionFlag": False, # peak not outside data range
"QuotedEtaFlag": False, # force 0 < eta < 1
# Peak detection
"AutoScaling": False,
"Yscaling": 1.0,
"FwhmPoints": 8,
"AutoFwhm": True,
"Sensitivity": 2.5,
"ForcePeakPresence": True,
# Hypermet
"HypermetTails": 15,
"QuotedFwhmFlag": 0,
"MaxFwhm2InputRatio": 1.5,
"MinFwhm2InputRatio": 0.4,
# short tail parameters
"MinGaussArea4ShortTail": 50000.0,
"InitialShortTailAreaRatio": 0.050,
"MaxShortTailAreaRatio": 0.100,
"MinShortTailAreaRatio": 0.0010,
"InitialShortTailSlopeRatio": 0.70,
"MaxShortTailSlopeRatio": 2.00,
"MinShortTailSlopeRatio": 0.50,
# long tail parameters
"MinGaussArea4LongTail": 1000.0,
"InitialLongTailAreaRatio": 0.050,
"MaxLongTailAreaRatio": 0.300,
"MinLongTailAreaRatio": 0.010,
"InitialLongTailSlopeRatio": 20.0,
"MaxLongTailSlopeRatio": 50.0,
"MinLongTailSlopeRatio": 5.0,
# step tail
"MinGaussHeight4StepTail": 5000.0,
"InitialStepTailHeightRatio": 0.002,
"MaxStepTailHeightRatio": 0.0100,
"MinStepTailHeightRatio": 0.0001,
# Hypermet constraints
# position in range [estimated position +- estimated fwhm/2]
"HypermetQuotedPositionFlag": True,
"DeltaPositionFwhmUnits": 0.5,
"SameSlopeRatioFlag": 1,
"SameAreaRatioFlag": 1,
# Strip bg removal
"StripBackgroundFlag": True,
"SmoothingFlag": True,
"SmoothingWidth": 5,
"StripWidth": 2,
"StripIterations": 5000,
"StripThresholdFactor": 1.0,
}
"""This dictionary defines default configuration parameters that have effects
on fit functions and estimation functions, mainly on fit constraints.
This dictionary is accessible as attribute :attr:`FitTheories.config`,
which can be modified by configuration functions defined in
:const:`CONFIGURE`.
"""
CFREE = 0
CPOSITIVE = 1
CQUOTED = 2
CFIXED = 3
CFACTOR = 4
CDELTA = 5
CSUM = 6
CIGNORED = 7
[docs]
class FitTheories(object):
"""Class wrapping functions from :class:`silx.math.fit.functions`
and providing estimate functions for all of these fit functions."""
def __init__(self, config=None):
if config is None:
self.config = DEFAULT_CONFIG
else:
self.config = config
[docs]
def ahypermet(self, x, *pars):
"""
Wrapping of :func:`silx.math.fit.functions.sum_ahypermet` without
the tail flags in the function signature.
Depending on the value of `self.config['HypermetTails']`, one can
activate or deactivate the various terms of the hypermet function.
`self.config['HypermetTails']` must be an integer between 0 and 15.
It is a set of 4 binary flags, one for activating each one of the
hypermet terms: *gaussian function, short tail, long tail, step*.
For example, 15 can be expressed as ``1111`` in base 2, so a flag of
15 means all terms are active.
"""
g_term = self.config["HypermetTails"] & 1
st_term = (self.config["HypermetTails"] >> 1) & 1
lt_term = (self.config["HypermetTails"] >> 2) & 1
step_term = (self.config["HypermetTails"] >> 3) & 1
return functions.sum_ahypermet(
x,
*pars,
gaussian_term=g_term,
st_term=st_term,
lt_term=lt_term,
step_term=step_term,
)
[docs]
def poly(self, x, *pars):
"""Order n polynomial.
The order of the polynomial is defined by the number of
coefficients (``*pars``).
"""
p = numpy.poly1d(pars)
return p(x)
[docs]
@staticmethod
def estimate_poly(x, y, n=2):
"""Estimate polynomial coefficients for a degree n polynomial."""
pcoeffs = numpy.polyfit(x, y, n)
constraints = numpy.zeros((n + 1, 3), numpy.float64)
return pcoeffs, constraints
[docs]
def estimate_quadratic(self, x, y):
"""Estimate quadratic coefficients"""
return self.estimate_poly(x, y, n=2)
[docs]
def estimate_cubic(self, x, y):
"""Estimate coefficients for a degree 3 polynomial"""
return self.estimate_poly(x, y, n=3)
[docs]
def estimate_quartic(self, x, y):
"""Estimate coefficients for a degree 4 polynomial"""
return self.estimate_poly(x, y, n=4)
[docs]
def estimate_quintic(self, x, y):
"""Estimate coefficients for a degree 5 polynomial"""
return self.estimate_poly(x, y, n=5)
[docs]
def strip_bg(self, y):
"""Return the strip background of y, using parameters from
:attr:`config` dictionary (*StripBackgroundFlag, StripWidth,
StripIterations, StripThresholdFactor*)"""
remove_strip_bg = self.config.get("StripBackgroundFlag", False)
if remove_strip_bg:
if self.config["SmoothingFlag"]:
y = savitsky_golay(y, self.config["SmoothingWidth"])
strip_width = self.config["StripWidth"]
strip_niterations = self.config["StripIterations"]
strip_thr_factor = self.config["StripThresholdFactor"]
return strip(
y, w=strip_width, niterations=strip_niterations, factor=strip_thr_factor
)
else:
return numpy.zeros_like(y)
[docs]
def guess_yscaling(self, y):
"""Estimate scaling for y prior to peak search.
A smoothing filter is applied to y to estimate the noise level
(chi-squared)
:param y: Data array
:return: Scaling factor
"""
# ensure y is an array
yy = numpy.asarray(y)
# smooth
convolution_kernel = numpy.ones(shape=(3,)) / 3.0
ysmooth = numpy.convolve(y, convolution_kernel, mode="same")
# remove zeros
idx_array = numpy.fabs(y) > 0.0
yy = yy[idx_array]
ysmooth = ysmooth[idx_array]
# compute scaling factor
chisq = numpy.mean((yy - ysmooth) ** 2 / numpy.fabs(yy))
if chisq > 0:
return 1.0 / chisq
else:
return 1.0
[docs]
def peak_search(self, y, fwhm, sensitivity):
"""Search for peaks in y array, after padding the array and
multiplying its value by a scaling factor.
:param y: 1-D data array
:param int fwhm: Typical full width at half maximum for peaks,
in number of points. This parameter is used for to discriminate between
true peaks and background fluctuations.
:param float sensitivity: Sensitivity parameter. This is a threshold factor
for peak detection. Only peaks larger than the standard deviation
of the noise multiplied by this sensitivity parameter are detected.
:return: List of peak indices
"""
# add padding
ysearch = numpy.ones((len(y) + 2 * fwhm,), numpy.float64)
ysearch[0:fwhm] = y[0]
ysearch[-1 : -fwhm - 1 : -1] = y[len(y) - 1]
ysearch[fwhm : fwhm + len(y)] = y[:]
scaling = (
self.guess_yscaling(y)
if self.config["AutoScaling"]
else self.config["Yscaling"]
)
if len(ysearch) > 1.5 * fwhm:
peaks = peak_search(scaling * ysearch, fwhm=fwhm, sensitivity=sensitivity)
return [
peak_index - fwhm
for peak_index in peaks
if 0 <= peak_index - fwhm < len(y)
]
else:
return []
[docs]
def estimate_height_position_fwhm(self, x, y):
"""Estimation of *Height, Position, FWHM* of peaks, for gaussian-like
curves.
This functions finds how many parameters are needed, based on the
number of peaks detected. Then it estimates the fit parameters
with a few iterations of fitting gaussian functions.
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
Parameters to be estimated for each peak are:
*Height, Position, FWHM*.
Fit constraints depend on :attr:`config`.
"""
fittedpar = []
bg = self.strip_bg(y)
if self.config["AutoFwhm"]:
search_fwhm = guess_fwhm(y)
else:
search_fwhm = int(float(self.config["FwhmPoints"]))
search_sens = float(self.config["Sensitivity"])
if search_fwhm < 3:
_logger.warning("Setting peak fwhm to 3 (lower limit)")
search_fwhm = 3
self.config["FwhmPoints"] = 3
if search_sens < 1:
_logger.warning(
"Setting peak search sensitivity to 1. "
+ "(lower limit to filter out noise peaks)"
)
search_sens = 1
self.config["Sensitivity"] = 1
npoints = len(y)
# Find indices of peaks in data array
peaks = self.peak_search(y, fwhm=search_fwhm, sensitivity=search_sens)
if not len(peaks):
forcepeak = int(float(self.config.get("ForcePeakPresence", 0)))
if forcepeak:
delta = y - bg
# get index of global maximum
# (first one if several samples are equal to this value)
peaks = [numpy.nonzero(delta == delta.max())[0][0]]
# Find index of largest peak in peaks array
index_largest_peak = 0
if len(peaks) > 0:
# estimate fwhm as 5 * sampling interval
sig = 5 * abs(x[npoints - 1] - x[0]) / npoints
peakpos = x[int(peaks[0])]
if abs(peakpos) < 1.0e-16:
peakpos = 0.0
param = numpy.array([y[int(peaks[0])] - bg[int(peaks[0])], peakpos, sig])
height_largest_peak = param[0]
peak_index = 1
for i in peaks[1:]:
param2 = numpy.array([y[int(i)] - bg[int(i)], x[int(i)], sig])
param = numpy.concatenate((param, param2))
if param2[0] > height_largest_peak:
height_largest_peak = param2[0]
index_largest_peak = peak_index
peak_index += 1
# Subtract background
xw = x
yw = y - bg
cons = numpy.zeros((len(param), 3), numpy.float64)
# peak height must be positive
cons[0 : len(param) : 3, 0] = CPOSITIVE
# force peaks to stay around their position
cons[1 : len(param) : 3, 0] = CQUOTED
# set possible peak range to estimated peak +- guessed fwhm
if len(xw) > search_fwhm:
fwhmx = numpy.fabs(xw[int(search_fwhm)] - xw[0])
cons[1 : len(param) : 3, 1] = param[1 : len(param) : 3] - 0.5 * fwhmx
cons[1 : len(param) : 3, 2] = param[1 : len(param) : 3] + 0.5 * fwhmx
else:
shape = [max(1, int(x)) for x in (param[1 : len(param) : 3])]
cons[1 : len(param) : 3, 1] = min(xw) * numpy.ones(shape, numpy.float64)
cons[1 : len(param) : 3, 2] = max(xw) * numpy.ones(shape, numpy.float64)
# ensure fwhm is positive
cons[2 : len(param) : 3, 0] = CPOSITIVE
# run a quick iterative fit (4 iterations) to improve
# estimations
fittedpar, _, _ = leastsq(
functions.sum_gauss,
xw,
yw,
param,
max_iter=4,
constraints=cons.tolist(),
full_output=True,
)
# set final constraints based on config parameters
cons = numpy.zeros((len(fittedpar), 3), numpy.float64)
peak_index = 0
for i in range(len(peaks)):
# Setup height area constrains
if not self.config["NoConstraintsFlag"]:
if self.config["PositiveHeightAreaFlag"]:
cons[peak_index, 0] = CPOSITIVE
cons[peak_index, 1] = 0
cons[peak_index, 2] = 0
peak_index += 1
# Setup position constrains
if not self.config["NoConstraintsFlag"]:
if self.config["QuotedPositionFlag"]:
cons[peak_index, 0] = CQUOTED
cons[peak_index, 1] = min(x)
cons[peak_index, 2] = max(x)
peak_index += 1
# Setup positive FWHM constrains
if not self.config["NoConstraintsFlag"]:
if self.config["PositiveFwhmFlag"]:
cons[peak_index, 0] = CPOSITIVE
cons[peak_index, 1] = 0
cons[peak_index, 2] = 0
if self.config["SameFwhmFlag"]:
if i != index_largest_peak:
cons[peak_index, 0] = CFACTOR
cons[peak_index, 1] = 3 * index_largest_peak + 2
cons[peak_index, 2] = 1.0
peak_index += 1
return fittedpar, cons
[docs]
def estimate_agauss(self, x, y):
"""Estimation of *Area, Position, FWHM* of peaks, for gaussian-like
curves.
This functions uses :meth:`estimate_height_position_fwhm`, then
converts the height parameters to area under the curve with the
formula ``area = sqrt(2*pi) * height * fwhm / (2 * sqrt(2 * log(2))``
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
Parameters to be estimated for each peak are:
*Area, Position, FWHM*.
Fit constraints depend on :attr:`config`.
"""
fittedpar, cons = self.estimate_height_position_fwhm(x, y)
# get the number of found peaks
npeaks = len(fittedpar) // 3
for i in range(npeaks):
height = fittedpar[3 * i]
fwhm = fittedpar[3 * i + 2]
# Replace height with area in fittedpar
fittedpar[3 * i] = (
numpy.sqrt(2 * numpy.pi)
* height
* fwhm
/ (2.0 * numpy.sqrt(2 * numpy.log(2)))
)
return fittedpar, cons
[docs]
def estimate_alorentz(self, x, y):
"""Estimation of *Area, Position, FWHM* of peaks, for Lorentzian
curves.
This functions uses :meth:`estimate_height_position_fwhm`, then
converts the height parameters to area under the curve with the
formula ``area = height * fwhm * 0.5 * pi``
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
Parameters to be estimated for each peak are:
*Area, Position, FWHM*.
Fit constraints depend on :attr:`config`.
"""
fittedpar, cons = self.estimate_height_position_fwhm(x, y)
# get the number of found peaks
npeaks = len(fittedpar) // 3
for i in range(npeaks):
height = fittedpar[3 * i]
fwhm = fittedpar[3 * i + 2]
# Replace height with area in fittedpar
fittedpar[3 * i] = height * fwhm * 0.5 * numpy.pi
return fittedpar, cons
[docs]
def estimate_splitgauss(self, x, y):
"""Estimation of *Height, Position, FWHM1, FWHM2* of peaks, for
asymmetric gaussian-like curves.
This functions uses :meth:`estimate_height_position_fwhm`, then
adds a second (identical) estimation of FWHM to the fit parameters
for each peak, and the corresponding constraint.
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
Parameters to be estimated for each peak are:
*Height, Position, FWHM1, FWHM2*.
Fit constraints depend on :attr:`config`.
"""
fittedpar, cons = self.estimate_height_position_fwhm(x, y)
# get the number of found peaks
npeaks = len(fittedpar) // 3
estimated_parameters = []
estimated_constraints = numpy.zeros((4 * npeaks, 3), numpy.float64)
for i in range(npeaks):
for j in range(3):
estimated_parameters.append(fittedpar[3 * i + j])
# fwhm2 estimate = fwhm1
estimated_parameters.append(fittedpar[3 * i + 2])
# height
estimated_constraints[4 * i, 0] = cons[3 * i, 0]
estimated_constraints[4 * i, 1] = cons[3 * i, 1]
estimated_constraints[4 * i, 2] = cons[3 * i, 2]
# position
estimated_constraints[4 * i + 1, 0] = cons[3 * i + 1, 0]
estimated_constraints[4 * i + 1, 1] = cons[3 * i + 1, 1]
estimated_constraints[4 * i + 1, 2] = cons[3 * i + 1, 2]
# fwhm1
estimated_constraints[4 * i + 2, 0] = cons[3 * i + 2, 0]
estimated_constraints[4 * i + 2, 1] = cons[3 * i + 2, 1]
estimated_constraints[4 * i + 2, 2] = cons[3 * i + 2, 2]
# fwhm2
estimated_constraints[4 * i + 3, 0] = cons[3 * i + 2, 0]
estimated_constraints[4 * i + 3, 1] = cons[3 * i + 2, 1]
estimated_constraints[4 * i + 3, 2] = cons[3 * i + 2, 2]
if cons[3 * i + 2, 0] == CFACTOR:
# convert indices of related parameters
# (this happens if SameFwhmFlag == True)
estimated_constraints[4 * i + 2, 1] = (
int(cons[3 * i + 2, 1] / 3) * 4 + 2
)
estimated_constraints[4 * i + 3, 1] = (
int(cons[3 * i + 2, 1] / 3) * 4 + 3
)
return estimated_parameters, estimated_constraints
[docs]
def estimate_pvoigt(self, x, y):
"""Estimation of *Height, Position, FWHM, eta* of peaks, for
pseudo-Voigt curves.
Pseudo-Voigt are a sum of a gaussian curve *G(x)* and a lorentzian
curve *L(x)* with the same height, center, fwhm parameters:
``y(x) = eta * G(x) + (1-eta) * L(x)``
This functions uses :meth:`estimate_height_position_fwhm`, then
adds a constant estimation of *eta* (0.5) to the fit parameters
for each peak, and the corresponding constraint.
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
Parameters to be estimated for each peak are:
*Height, Position, FWHM, eta*.
Constraint for the eta parameter can be set to QUOTED (0.--1.)
by setting :attr:`config`['QuotedEtaFlag'] to ``True``.
If this is not the case, the constraint code is set to FREE.
"""
fittedpar, cons = self.estimate_height_position_fwhm(x, y)
npeaks = len(fittedpar) // 3
newpar = []
newcons = numpy.zeros((4 * npeaks, 3), numpy.float64)
# find out related parameters proper index
if not self.config["NoConstraintsFlag"]:
if self.config["SameFwhmFlag"]:
j = 0
# get the index of the free FWHM
for i in range(npeaks):
if cons[3 * i + 2, 0] != 4:
j = i
for i in range(npeaks):
if i != j:
cons[3 * i + 2, 1] = 4 * j + 2
for i in range(npeaks):
newpar.append(fittedpar[3 * i])
newpar.append(fittedpar[3 * i + 1])
newpar.append(fittedpar[3 * i + 2])
newpar.append(0.5)
# height
newcons[4 * i, 0] = cons[3 * i, 0]
newcons[4 * i, 1] = cons[3 * i, 1]
newcons[4 * i, 2] = cons[3 * i, 2]
# position
newcons[4 * i + 1, 0] = cons[3 * i + 1, 0]
newcons[4 * i + 1, 1] = cons[3 * i + 1, 1]
newcons[4 * i + 1, 2] = cons[3 * i + 1, 2]
# fwhm
newcons[4 * i + 2, 0] = cons[3 * i + 2, 0]
newcons[4 * i + 2, 1] = cons[3 * i + 2, 1]
newcons[4 * i + 2, 2] = cons[3 * i + 2, 2]
# Eta constrains
newcons[4 * i + 3, 0] = CFREE
newcons[4 * i + 3, 1] = 0
newcons[4 * i + 3, 2] = 0
if self.config["QuotedEtaFlag"]:
newcons[4 * i + 3, 0] = CQUOTED
newcons[4 * i + 3, 1] = 0.0
newcons[4 * i + 3, 2] = 1.0
return newpar, newcons
[docs]
def estimate_splitpvoigt(self, x, y):
"""Estimation of *Height, Position, FWHM1, FWHM2, eta* of peaks, for
asymmetric pseudo-Voigt curves.
This functions uses :meth:`estimate_height_position_fwhm`, then
adds an identical FWHM2 parameter and a constant estimation of
*eta* (0.5) to the fit parameters for each peak, and the corresponding
constraints.
Constraint for the eta parameter can be set to QUOTED (0.--1.)
by setting :attr:`config`['QuotedEtaFlag'] to ``True``.
If this is not the case, the constraint code is set to FREE.
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
Parameters to be estimated for each peak are:
*Height, Position, FWHM1, FWHM2, eta*.
"""
fittedpar, cons = self.estimate_height_position_fwhm(x, y)
npeaks = len(fittedpar) // 3
newpar = []
newcons = numpy.zeros((5 * npeaks, 3), numpy.float64)
# find out related parameters proper index
if not self.config["NoConstraintsFlag"]:
if self.config["SameFwhmFlag"]:
j = 0
# get the index of the free FWHM
for i in range(npeaks):
if cons[3 * i + 2, 0] != 4:
j = i
for i in range(npeaks):
if i != j:
cons[3 * i + 2, 1] = 4 * j + 2
for i in range(npeaks):
# height
newpar.append(fittedpar[3 * i])
# position
newpar.append(fittedpar[3 * i + 1])
# fwhm1
newpar.append(fittedpar[3 * i + 2])
# fwhm2 estimate equal to fwhm1
newpar.append(fittedpar[3 * i + 2])
# eta
newpar.append(0.5)
# constraint codes
# ----------------
# height
newcons[5 * i, 0] = cons[3 * i, 0]
# position
newcons[5 * i + 1, 0] = cons[3 * i + 1, 0]
# fwhm1
newcons[5 * i + 2, 0] = cons[3 * i + 2, 0]
# fwhm2
newcons[5 * i + 3, 0] = cons[3 * i + 2, 0]
# cons 1
# ------
newcons[5 * i, 1] = cons[3 * i, 1]
newcons[5 * i + 1, 1] = cons[3 * i + 1, 1]
newcons[5 * i + 2, 1] = cons[3 * i + 2, 1]
newcons[5 * i + 3, 1] = cons[3 * i + 2, 1]
# cons 2
# ------
newcons[5 * i, 2] = cons[3 * i, 2]
newcons[5 * i + 1, 2] = cons[3 * i + 1, 2]
newcons[5 * i + 2, 2] = cons[3 * i + 2, 2]
newcons[5 * i + 3, 2] = cons[3 * i + 2, 2]
if cons[3 * i + 2, 0] == CFACTOR:
# fwhm2 connstraint depends on fwhm1
newcons[5 * i + 3, 1] = newcons[5 * i + 2, 1] + 1
# eta constraints
newcons[5 * i + 4, 0] = CFREE
newcons[5 * i + 4, 1] = 0
newcons[5 * i + 4, 2] = 0
if self.config["QuotedEtaFlag"]:
newcons[5 * i + 4, 0] = CQUOTED
newcons[5 * i + 4, 1] = 0.0
newcons[5 * i + 4, 2] = 1.0
return newpar, newcons
[docs]
def estimate_splitpvoigt2(self, x, y):
"""Estimation of *Height, Position, FWHM1, FWHM2, eta1, eta2* of peaks, for
asymmetric pseudo-Voigt curves.
This functions uses :meth:`estimate_height_position_fwhm`, then
adds an identical FWHM2 parameter and a constant estimation of
*eta1* and *eta2* (0.5) to the fit parameters for each peak, and the corresponding
constraints.
Constraint for the eta parameter can be set to QUOTED (0.--1.)
by setting :attr:`config`['QuotedEtaFlag'] to ``True``.
If this is not the case, the constraint code is set to FREE.
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
Parameters to be estimated for each peak are:
*Height, Position, FWHM1, FWHM2, eta1, eta2*.
"""
fittedpar, cons = self.estimate_height_position_fwhm(x, y)
npeaks = len(fittedpar) // 3
newpar = []
newcons = numpy.zeros((5 * npeaks, 3), numpy.float64)
# find out related parameters proper index
if not self.config["NoConstraintsFlag"]:
if self.config["SameFwhmFlag"]:
j = 0
# get the index of the free FWHM
for i in range(npeaks):
if cons[3 * i + 2, 0] != 4:
j = i
for i in range(npeaks):
if i != j:
cons[3 * i + 2, 1] = 4 * j + 2
for i in range(npeaks):
# height
newpar.append(fittedpar[3 * i])
# position
newpar.append(fittedpar[3 * i + 1])
# fwhm1
newpar.append(fittedpar[3 * i + 2])
# fwhm2 estimate equal to fwhm1
newpar.append(fittedpar[3 * i + 2])
# eta1
newpar.append(0.5)
# eta2
newpar.append(0.5)
# constraint codes
# ----------------
# height
newcons[6 * i, 0] = cons[3 * i, 0]
# position
newcons[6 * i + 1, 0] = cons[3 * i + 1, 0]
# fwhm1
newcons[6 * i + 2, 0] = cons[3 * i + 2, 0]
# fwhm2
newcons[6 * i + 3, 0] = cons[3 * i + 2, 0]
# cons 1
# ------
newcons[6 * i, 1] = cons[3 * i, 1]
newcons[6 * i + 1, 1] = cons[3 * i + 1, 1]
newcons[6 * i + 2, 1] = cons[3 * i + 2, 1]
newcons[6 * i + 3, 1] = cons[3 * i + 2, 1]
# cons 2
# ------
newcons[6 * i, 2] = cons[3 * i, 2]
newcons[6 * i + 1, 2] = cons[3 * i + 1, 2]
newcons[6 * i + 2, 2] = cons[3 * i + 2, 2]
newcons[6 * i + 3, 2] = cons[3 * i + 2, 2]
if cons[3 * i + 2, 0] == CFACTOR:
# fwhm2 constraint depends on fwhm1
newcons[6 * i + 3, 1] = newcons[6 * i + 2, 1] + 1
# eta1 constraints
newcons[6 * i + 4, 0] = CFREE
newcons[6 * i + 4, 1] = 0
newcons[6 * i + 4, 2] = 0
if self.config["QuotedEtaFlag"]:
newcons[6 * i + 4, 0] = CQUOTED
newcons[6 * i + 4, 1] = 0.0
newcons[6 * i + 4, 2] = 1.0
# eta2 constraints
newcons[6 * i + 5, 0] = CFREE
newcons[6 * i + 5, 1] = 0
newcons[6 * i + 5, 2] = 0
if self.config["QuotedEtaFlag"]:
newcons[6 * i + 5, 0] = CQUOTED
newcons[6 * i + 5, 1] = 0.0
newcons[6 * i + 5, 2] = 1.0
return newpar, newcons
[docs]
def estimate_apvoigt(self, x, y):
"""Estimation of *Area, Position, FWHM1, eta* of peaks, for
pseudo-Voigt curves.
This functions uses :meth:`estimate_pvoigt`, then converts the height
parameter to area.
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
Parameters to be estimated for each peak are:
*Area, Position, FWHM, eta*.
"""
fittedpar, cons = self.estimate_pvoigt(x, y)
npeaks = len(fittedpar) // 4
# Assume 50% of the area is determined by the gaussian and 50% by
# the Lorentzian.
for i in range(npeaks):
height = fittedpar[4 * i]
fwhm = fittedpar[4 * i + 2]
fittedpar[4 * i] = 0.5 * (height * fwhm * 0.5 * numpy.pi) + 0.5 * (
height * fwhm / (2.0 * numpy.sqrt(2 * numpy.log(2)))
) * numpy.sqrt(2 * numpy.pi)
return fittedpar, cons
[docs]
def estimate_ahypermet(self, x, y):
"""Estimation of *area, position, fwhm, st_area_r, st_slope_r,
lt_area_r, lt_slope_r, step_height_r* of peaks, for hypermet curves.
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
Parameters to be estimated for each peak are:
*area, position, fwhm, st_area_r, st_slope_r,
lt_area_r, lt_slope_r, step_height_r* .
"""
yscaling = self.config.get("Yscaling", 1.0)
if yscaling == 0:
yscaling = 1.0
fittedpar, cons = self.estimate_height_position_fwhm(x, y)
npeaks = len(fittedpar) // 3
newpar = []
newcons = numpy.zeros((8 * npeaks, 3), numpy.float64)
main_peak = 0
# find out related parameters proper index
if not self.config["NoConstraintsFlag"]:
if self.config["SameFwhmFlag"]:
j = 0
# get the index of the free FWHM
for i in range(npeaks):
if cons[3 * i + 2, 0] != 4:
j = i
for i in range(npeaks):
if i != j:
cons[3 * i + 2, 1] = 8 * j + 2
main_peak = j
for i in range(npeaks):
if fittedpar[3 * i] > fittedpar[3 * main_peak]:
main_peak = i
for i in range(npeaks):
height = fittedpar[3 * i]
position = fittedpar[3 * i + 1]
fwhm = fittedpar[3 * i + 2]
area = (height * fwhm / (2.0 * numpy.sqrt(2 * numpy.log(2)))) * numpy.sqrt(
2 * numpy.pi
)
# the gaussian parameters
newpar.append(area)
newpar.append(position)
newpar.append(fwhm)
# print "area, pos , fwhm = ",area,position,fwhm
# Avoid zero derivatives because of not calculating contribution
g_term = 1
st_term = 1
lt_term = 1
step_term = 1
if self.config["HypermetTails"] != 0:
g_term = self.config["HypermetTails"] & 1
st_term = (self.config["HypermetTails"] >> 1) & 1
lt_term = (self.config["HypermetTails"] >> 2) & 1
step_term = (self.config["HypermetTails"] >> 3) & 1
if g_term == 0:
# fix the gaussian parameters
newcons[8 * i, 0] = CFIXED
newcons[8 * i + 1, 0] = CFIXED
newcons[8 * i + 2, 0] = CFIXED
# the short tail parameters
if ((area * yscaling) < self.config["MinGaussArea4ShortTail"]) | (
st_term == 0
):
newpar.append(0.0)
newpar.append(0.0)
newcons[8 * i + 3, 0] = CFIXED
newcons[8 * i + 3, 1] = 0.0
newcons[8 * i + 3, 2] = 0.0
newcons[8 * i + 4, 0] = CFIXED
newcons[8 * i + 4, 1] = 0.0
newcons[8 * i + 4, 2] = 0.0
else:
newpar.append(self.config["InitialShortTailAreaRatio"])
newpar.append(self.config["InitialShortTailSlopeRatio"])
newcons[8 * i + 3, 0] = CQUOTED
newcons[8 * i + 3, 1] = self.config["MinShortTailAreaRatio"]
newcons[8 * i + 3, 2] = self.config["MaxShortTailAreaRatio"]
newcons[8 * i + 4, 0] = CQUOTED
newcons[8 * i + 4, 1] = self.config["MinShortTailSlopeRatio"]
newcons[8 * i + 4, 2] = self.config["MaxShortTailSlopeRatio"]
# the long tail parameters
if ((area * yscaling) < self.config["MinGaussArea4LongTail"]) | (
lt_term == 0
):
newpar.append(0.0)
newpar.append(0.0)
newcons[8 * i + 5, 0] = CFIXED
newcons[8 * i + 5, 1] = 0.0
newcons[8 * i + 5, 2] = 0.0
newcons[8 * i + 6, 0] = CFIXED
newcons[8 * i + 6, 1] = 0.0
newcons[8 * i + 6, 2] = 0.0
else:
newpar.append(self.config["InitialLongTailAreaRatio"])
newpar.append(self.config["InitialLongTailSlopeRatio"])
newcons[8 * i + 5, 0] = CQUOTED
newcons[8 * i + 5, 1] = self.config["MinLongTailAreaRatio"]
newcons[8 * i + 5, 2] = self.config["MaxLongTailAreaRatio"]
newcons[8 * i + 6, 0] = CQUOTED
newcons[8 * i + 6, 1] = self.config["MinLongTailSlopeRatio"]
newcons[8 * i + 6, 2] = self.config["MaxLongTailSlopeRatio"]
# the step parameters
if ((height * yscaling) < self.config["MinGaussHeight4StepTail"]) | (
step_term == 0
):
newpar.append(0.0)
newcons[8 * i + 7, 0] = CFIXED
newcons[8 * i + 7, 1] = 0.0
newcons[8 * i + 7, 2] = 0.0
else:
newpar.append(self.config["InitialStepTailHeightRatio"])
newcons[8 * i + 7, 0] = CQUOTED
newcons[8 * i + 7, 1] = self.config["MinStepTailHeightRatio"]
newcons[8 * i + 7, 2] = self.config["MaxStepTailHeightRatio"]
# if self.config['NoConstraintsFlag'] == 1:
# newcons=numpy.zeros((8*npeaks, 3),numpy.float64)
if npeaks > 0:
if g_term:
if self.config["PositiveHeightAreaFlag"]:
for i in range(npeaks):
newcons[8 * i, 0] = CPOSITIVE
if self.config["PositiveFwhmFlag"]:
for i in range(npeaks):
newcons[8 * i + 2, 0] = CPOSITIVE
if self.config["SameFwhmFlag"]:
for i in range(npeaks):
if i != main_peak:
newcons[8 * i + 2, 0] = CFACTOR
newcons[8 * i + 2, 1] = 8 * main_peak + 2
newcons[8 * i + 2, 2] = 1.0
if self.config["HypermetQuotedPositionFlag"]:
for i in range(npeaks):
delta = self.config["DeltaPositionFwhmUnits"] * fwhm
newcons[8 * i + 1, 0] = CQUOTED
newcons[8 * i + 1, 1] = newpar[8 * i + 1] - delta
newcons[8 * i + 1, 2] = newpar[8 * i + 1] + delta
if self.config["SameSlopeRatioFlag"]:
for i in range(npeaks):
if i != main_peak:
newcons[8 * i + 4, 0] = CFACTOR
newcons[8 * i + 4, 1] = 8 * main_peak + 4
newcons[8 * i + 4, 2] = 1.0
newcons[8 * i + 6, 0] = CFACTOR
newcons[8 * i + 6, 1] = 8 * main_peak + 6
newcons[8 * i + 6, 2] = 1.0
if self.config["SameAreaRatioFlag"]:
for i in range(npeaks):
if i != main_peak:
newcons[8 * i + 3, 0] = CFACTOR
newcons[8 * i + 3, 1] = 8 * main_peak + 3
newcons[8 * i + 3, 2] = 1.0
newcons[8 * i + 5, 0] = CFACTOR
newcons[8 * i + 5, 1] = 8 * main_peak + 5
newcons[8 * i + 5, 2] = 1.0
return newpar, newcons
[docs]
def estimate_stepdown(self, x, y):
"""Estimation of parameters for stepdown curves.
The functions estimates gaussian parameters for the derivative of
the data, takes the largest gaussian peak and uses its estimated
parameters to define the center of the step and its fwhm. The
estimated amplitude returned is simply ``max(y) - min(y)``.
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit newconstraints.
Parameters to be estimated for each stepdown are:
*height, centroid, fwhm* .
"""
crappyfilter = [-0.25, -0.75, 0.0, 0.75, 0.25]
cutoff = len(crappyfilter) // 2
y_deriv = numpy.convolve(y, crappyfilter, mode="valid")
# make the derivative's peak have the same amplitude as the step
if max(y_deriv) > 0:
y_deriv = y_deriv * max(y) / max(y_deriv)
fittedpar, newcons = self.estimate_height_position_fwhm(
x[cutoff:-cutoff], y_deriv
)
data_amplitude = max(y) - min(y)
# use parameters from largest gaussian found
if len(fittedpar):
npeaks = len(fittedpar) // 3
largest_index = 0
largest = [
data_amplitude,
fittedpar[3 * largest_index + 1],
fittedpar[3 * largest_index + 2],
]
for i in range(npeaks):
if fittedpar[3 * i] > largest[0]:
largest_index = i
largest = [
data_amplitude,
fittedpar[3 * largest_index + 1],
fittedpar[3 * largest_index + 2],
]
else:
# no peak was found
largest = [
data_amplitude, # height
x[len(x) // 2], # center: middle of x range
self.config["FwhmPoints"] * (x[1] - x[0]),
] # fwhm: default value
# Setup constrains
newcons = numpy.zeros((3, 3), numpy.float64)
if not self.config["NoConstraintsFlag"]:
# Setup height constrains
if self.config["PositiveHeightAreaFlag"]:
newcons[0, 0] = CPOSITIVE
newcons[0, 1] = 0
newcons[0, 2] = 0
# Setup position constrains
if self.config["QuotedPositionFlag"]:
newcons[1, 0] = CQUOTED
newcons[1, 1] = min(x)
newcons[1, 2] = max(x)
# Setup positive FWHM constrains
if self.config["PositiveFwhmFlag"]:
newcons[2, 0] = CPOSITIVE
newcons[2, 1] = 0
newcons[2, 2] = 0
return largest, newcons
[docs]
def estimate_slit(self, x, y):
"""Estimation of parameters for slit curves.
The functions estimates stepup and stepdown parameters for the largest
steps, and uses them for calculating the center (middle between stepup
and stepdown), the height (maximum amplitude in data), the fwhm
(distance between the up- and down-step centers) and the beamfwhm
(average of FWHM for up- and down-step).
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
Parameters to be estimated for each slit are:
*height, position, fwhm, beamfwhm* .
"""
largestup, cons = self.estimate_stepup(x, y)
largestdown, cons = self.estimate_stepdown(x, y)
fwhm = numpy.fabs(largestdown[1] - largestup[1])
beamfwhm = 0.5 * (largestup[2] + largestdown[1])
beamfwhm = min(beamfwhm, fwhm / 10.0)
beamfwhm = max(beamfwhm, (max(x) - min(x)) * 3.0 / len(x))
y_minus_bg = y - self.strip_bg(y)
height = max(y_minus_bg)
i1 = numpy.nonzero(y_minus_bg >= 0.5 * height)[0]
xx = numpy.take(x, i1)
position = (xx[0] + xx[-1]) / 2.0
fwhm = xx[-1] - xx[0]
largest = [height, position, fwhm, beamfwhm]
cons = numpy.zeros((4, 3), numpy.float64)
# Setup constrains
if not self.config["NoConstraintsFlag"]:
# Setup height constrains
if self.config["PositiveHeightAreaFlag"]:
cons[0, 0] = CPOSITIVE
cons[0, 1] = 0
cons[0, 2] = 0
# Setup position constrains
if self.config["QuotedPositionFlag"]:
cons[1, 0] = CQUOTED
cons[1, 1] = min(x)
cons[1, 2] = max(x)
# Setup positive FWHM constrains
if self.config["PositiveFwhmFlag"]:
cons[2, 0] = CPOSITIVE
cons[2, 1] = 0
cons[2, 2] = 0
# Setup positive FWHM constrains
if self.config["PositiveFwhmFlag"]:
cons[3, 0] = CPOSITIVE
cons[3, 1] = 0
cons[3, 2] = 0
return largest, cons
[docs]
def estimate_stepup(self, x, y):
"""Estimation of parameters for a single step up curve.
The functions estimates gaussian parameters for the derivative of
the data, takes the largest gaussian peak and uses its estimated
parameters to define the center of the step and its fwhm. The
estimated amplitude returned is simply ``max(y) - min(y)``.
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
Parameters to be estimated for each stepup are:
*height, centroid, fwhm* .
"""
crappyfilter = [0.25, 0.75, 0.0, -0.75, -0.25]
cutoff = len(crappyfilter) // 2
y_deriv = numpy.convolve(y, crappyfilter, mode="valid")
if max(y_deriv) > 0:
y_deriv = y_deriv * max(y) / max(y_deriv)
fittedpar, cons = self.estimate_height_position_fwhm(x[cutoff:-cutoff], y_deriv)
# for height, use the data amplitude after removing the background
data_amplitude = max(y) - min(y)
# find params of the largest gaussian found
if len(fittedpar):
npeaks = len(fittedpar) // 3
largest_index = 0
largest = [
data_amplitude,
fittedpar[3 * largest_index + 1],
fittedpar[3 * largest_index + 2],
]
for i in range(npeaks):
if fittedpar[3 * i] > largest[0]:
largest_index = i
largest = [
fittedpar[3 * largest_index],
fittedpar[3 * largest_index + 1],
fittedpar[3 * largest_index + 2],
]
else:
# no peak was found
largest = [
data_amplitude, # height
x[len(x) // 2], # center: middle of x range
self.config["FwhmPoints"] * (x[1] - x[0]),
] # fwhm: default value
newcons = numpy.zeros((3, 3), numpy.float64)
# Setup constrains
if not self.config["NoConstraintsFlag"]:
# Setup height constraints
if self.config["PositiveHeightAreaFlag"]:
newcons[0, 0] = CPOSITIVE
newcons[0, 1] = 0
newcons[0, 2] = 0
# Setup position constraints
if self.config["QuotedPositionFlag"]:
newcons[1, 0] = CQUOTED
newcons[1, 1] = min(x)
newcons[1, 2] = max(x)
# Setup positive FWHM constraints
if self.config["PositiveFwhmFlag"]:
newcons[2, 0] = CPOSITIVE
newcons[2, 1] = 0
newcons[2, 2] = 0
return largest, newcons
[docs]
def estimate_periodic_gauss(self, x, y):
"""Estimation of parameters for periodic gaussian curves:
*number of peaks, distance between peaks, height, position of the
first peak, fwhm*
The functions detects all peaks, then computes the parameters the
following way:
- *distance*: average of distances between detected peaks
- *height*: average height of detected peaks
- *fwhm*: fwhm of the highest peak (in number of samples) if
field ``'AutoFwhm'`` in :attr:`config` is ``True``, else take
the default value (field ``'FwhmPoints'`` in :attr:`config`)
:param x: Array of abscissa values
:param y: Array of ordinate values (``y = f(x)``)
:return: Tuple of estimated fit parameters and fit constraints.
"""
yscaling = self.config.get("Yscaling", 1.0)
if yscaling == 0:
yscaling = 1.0
bg = self.strip_bg(y)
if self.config["AutoFwhm"]:
search_fwhm = guess_fwhm(y)
else:
search_fwhm = int(float(self.config["FwhmPoints"]))
search_sens = float(self.config["Sensitivity"])
if search_fwhm < 3:
search_fwhm = 3
if search_sens < 1:
search_sens = 1
if len(y) > 1.5 * search_fwhm:
peaks = peak_search(yscaling * y, fwhm=search_fwhm, sensitivity=search_sens)
else:
peaks = []
npeaks = len(peaks)
if not npeaks:
fittedpar = []
cons = numpy.zeros((len(fittedpar), 3), numpy.float64)
return fittedpar, cons
fittedpar = [0.0, 0.0, 0.0, 0.0, 0.0]
# The number of peaks
fittedpar[0] = npeaks
# The separation between peaks in x units
delta = 0.0
height = 0.0
for i in range(npeaks):
height += y[int(peaks[i])] - bg[int(peaks[i])]
if i != npeaks - 1:
delta += x[int(peaks[i + 1])] - x[int(peaks[i])]
# delta between peaks
if npeaks > 1:
fittedpar[1] = delta / (npeaks - 1)
# starting height
fittedpar[2] = height / npeaks
# position of the first peak
fittedpar[3] = x[int(peaks[0])]
# Estimate the fwhm
fittedpar[4] = search_fwhm
# setup constraints
cons = numpy.zeros((5, 3), numpy.float64)
cons[0, 0] = CFIXED # the number of gaussians
if npeaks == 1:
cons[1, 0] = CFIXED # the delta between peaks
else:
cons[1, 0] = CFREE
j = 2
# Setup height area constrains
if not self.config["NoConstraintsFlag"]:
if self.config["PositiveHeightAreaFlag"]:
# POSITIVE = 1
cons[j, 0] = CPOSITIVE
cons[j, 1] = 0
cons[j, 2] = 0
j += 1
# Setup position constrains
if not self.config["NoConstraintsFlag"]:
if self.config["QuotedPositionFlag"]:
# QUOTED = 2
cons[j, 0] = CQUOTED
cons[j, 1] = min(x)
cons[j, 2] = max(x)
j += 1
# Setup positive FWHM constrains
if not self.config["NoConstraintsFlag"]:
if self.config["PositiveFwhmFlag"]:
# POSITIVE=1
cons[j, 0] = CPOSITIVE
cons[j, 1] = 0
cons[j, 2] = 0
j += 1
return fittedpar, cons
fitfuns = FitTheories()
THEORY = dict(
(
(
"Gaussians",
FitTheory(
description="Gaussian functions",
function=functions.sum_gauss,
parameters=("Height", "Position", "FWHM"),
estimate=fitfuns.estimate_height_position_fwhm,
configure=fitfuns.configure,
),
),
(
"Lorentz",
FitTheory(
description="Lorentzian functions",
function=functions.sum_lorentz,
parameters=("Height", "Position", "FWHM"),
estimate=fitfuns.estimate_height_position_fwhm,
configure=fitfuns.configure,
),
),
(
"Area Gaussians",
FitTheory(
description="Gaussian functions (area)",
function=functions.sum_agauss,
parameters=("Area", "Position", "FWHM"),
estimate=fitfuns.estimate_agauss,
configure=fitfuns.configure,
),
),
(
"Area Lorentz",
FitTheory(
description="Lorentzian functions (area)",
function=functions.sum_alorentz,
parameters=("Area", "Position", "FWHM"),
estimate=fitfuns.estimate_alorentz,
configure=fitfuns.configure,
),
),
(
"Pseudo-Voigt Line",
FitTheory(
description="Pseudo-Voigt functions",
function=functions.sum_pvoigt,
parameters=("Height", "Position", "FWHM", "Eta"),
estimate=fitfuns.estimate_pvoigt,
configure=fitfuns.configure,
),
),
(
"Area Pseudo-Voigt",
FitTheory(
description="Pseudo-Voigt functions (area)",
function=functions.sum_apvoigt,
parameters=("Area", "Position", "FWHM", "Eta"),
estimate=fitfuns.estimate_apvoigt,
configure=fitfuns.configure,
),
),
(
"Split Gaussian",
FitTheory(
description="Asymmetric gaussian functions",
function=functions.sum_splitgauss,
parameters=("Height", "Position", "LowFWHM", "HighFWHM"),
estimate=fitfuns.estimate_splitgauss,
configure=fitfuns.configure,
),
),
(
"Split Lorentz",
FitTheory(
description="Asymmetric lorentzian functions",
function=functions.sum_splitlorentz,
parameters=("Height", "Position", "LowFWHM", "HighFWHM"),
estimate=fitfuns.estimate_splitgauss,
configure=fitfuns.configure,
),
),
(
"Split Pseudo-Voigt",
FitTheory(
description="Asymmetric pseudo-Voigt functions",
function=functions.sum_splitpvoigt,
parameters=("Height", "Position", "LowFWHM", "HighFWHM", "Eta"),
estimate=fitfuns.estimate_splitpvoigt,
configure=fitfuns.configure,
),
),
(
"Split Pseudo-Voigt 2",
FitTheory(
description="Asymmetric pseudo-Voigt functions",
function=functions.sum_splitpvoigt2,
parameters=(
"Height",
"Position",
"LowFWHM",
"HighFWHM",
"LowEta",
"HighEta",
),
estimate=fitfuns.estimate_splitpvoigt2,
configure=fitfuns.configure,
),
),
(
"Step Down",
FitTheory(
description="Step down function",
function=functions.sum_stepdown,
parameters=("Height", "Position", "FWHM"),
estimate=fitfuns.estimate_stepdown,
configure=fitfuns.configure,
),
),
(
"Step Up",
FitTheory(
description="Step up function",
function=functions.sum_stepup,
parameters=("Height", "Position", "FWHM"),
estimate=fitfuns.estimate_stepup,
configure=fitfuns.configure,
),
),
(
"Slit",
FitTheory(
description="Slit function",
function=functions.sum_slit,
parameters=("Height", "Position", "FWHM", "BeamFWHM"),
estimate=fitfuns.estimate_slit,
configure=fitfuns.configure,
),
),
(
"Atan",
FitTheory(
description="Arctan step up function",
function=functions.atan_stepup,
parameters=("Height", "Position", "Width"),
estimate=fitfuns.estimate_stepup,
configure=fitfuns.configure,
),
),
(
"Hypermet",
FitTheory(
description="Hypermet functions",
function=fitfuns.ahypermet, # customized version of functions.sum_ahypermet
parameters=(
"G_Area",
"Position",
"FWHM",
"ST_Area",
"ST_Slope",
"LT_Area",
"LT_Slope",
"Step_H",
),
estimate=fitfuns.estimate_ahypermet,
configure=fitfuns.configure,
),
),
# ('Periodic Gaussians',
# FitTheory(description='Periodic gaussian functions',
# function=functions.periodic_gauss,
# parameters=('N', 'Delta', 'Height', 'Position', 'FWHM'),
# estimate=fitfuns.estimate_periodic_gauss,
# configure=fitfuns.configure))
(
"Degree 2 Polynomial",
FitTheory(
description="Degree 2 polynomial" "\ny = a*x^2 + b*x +c",
function=fitfuns.poly,
parameters=["a", "b", "c"],
estimate=fitfuns.estimate_quadratic,
),
),
(
"Degree 3 Polynomial",
FitTheory(
description="Degree 3 polynomial" "\ny = a*x^3 + b*x^2 + c*x + d",
function=fitfuns.poly,
parameters=["a", "b", "c", "d"],
estimate=fitfuns.estimate_cubic,
),
),
(
"Degree 4 Polynomial",
FitTheory(
description="Degree 4 polynomial"
"\ny = a*x^4 + b*x^3 + c*x^2 + d*x + e",
function=fitfuns.poly,
parameters=["a", "b", "c", "d", "e"],
estimate=fitfuns.estimate_quartic,
),
),
(
"Degree 5 Polynomial",
FitTheory(
description="Degree 5 polynomial"
"\ny = a*x^5 + b*x^4 + c*x^3 + d*x^2 + e*x + f",
function=fitfuns.poly,
parameters=["a", "b", "c", "d", "e", "f"],
estimate=fitfuns.estimate_quintic,
),
),
)
)
"""Dictionary of fit theories: fit functions and their associated estimation
function, parameters list, configuration function and description.
"""
def test(a):
from silx.math.fit import fitmanager
x = numpy.arange(1000).astype(numpy.float64)
p = [1500, 100.0, 50.0, 1500, 700.0, 50.0]
y_synthetic = functions.sum_gauss(x, *p) + 1
fit = fitmanager.FitManager(x, y_synthetic)
fit.addtheory(
"Gaussians",
functions.sum_gauss,
["Height", "Position", "FWHM"],
a.estimate_height_position_fwhm,
)
fit.settheory("Gaussians")
fit.setbackground("Linear")
fit.estimate()
fit.runfit()
y_fit = fit.gendata()
print("Fit parameter names: %s" % str(fit.get_names()))
print("Theoretical parameters: %s" % str(numpy.append([1, 0], p)))
print("Fitted parameters: %s" % str(fit.get_fitted_parameters()))
try:
from silx.gui import qt
from silx.gui.plot import plot1D
app = qt.QApplication([])
# Offset of 1 to see the difference in log scale
plot1D(x, (y_synthetic + 1, y_fit), "Input data + 1, Fit")
app.exec()
except ImportError:
_logger.warning("Unable to load qt binding, can't plot results.")
if __name__ == "__main__":
test(fitfuns)
