diff options
Diffstat (limited to 'src/silx/math/fit/bgtheories.py')
| -rw-r--r-- | src/silx/math/fit/bgtheories.py | 459 |
1 files changed, 459 insertions, 0 deletions
diff --git a/src/silx/math/fit/bgtheories.py b/src/silx/math/fit/bgtheories.py new file mode 100644 index 0000000..e698927 --- /dev/null +++ b/src/silx/math/fit/bgtheories.py @@ -0,0 +1,459 @@ +# /*########################################################################## +# +# 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 defines a set of background model functions and associated +estimation functions in a format that can be imported into a +:class:`silx.math.fit.FitManager` object. + +A background function is a function that you want to add to a regular fit +function prior to fitting the sum of both functions. This is useful, for +instance, if you need to fit multiple gaussian peaks in an array of +measured data points when the measurement is polluted by a background signal. + +The models include common background models such as a constant value or a +linear background. + +It also includes background computation filters - *strip* and *snip* - that +can extract a more complex low-curvature background signal from a signal with +peaks having higher curvatures. + +The source code of this module can serve as a template for defining your +own fit background theories. The minimal skeleton of such a theory definition +file is:: + + from silx.math.fit.fittheory import FitTheory + + def bgfunction1(x, y0, …): + bg_signal = … + return bg_signal + + def estimation_function1(x, y): + … + estimated_params = … + constraints = … + return estimated_params, constraints + + THEORY = { + 'bg_theory_name1': FitTheory( + description='Description of theory 1', + function=bgfunction1, + parameters=('param name 1', 'param name 2', …), + estimate=estimation_function1, + configure=configuration_function1, + derivative=derivative_function1, + is_background=True), + 'theory_name_2': …, + } +""" +__authors__ = ["P. Knobel"] +__license__ = "MIT" +__date__ = "16/01/2017" + +import numpy +from silx.math.fit.filters import strip, snip1d, savitsky_golay +from silx.math.fit.fittheory import FitTheory + +CONFIG = { + "SmoothingFlag": False, + "SmoothingWidth": 5, + "AnchorsFlag": False, + "AnchorsList": [], + "StripWidth": 2, + "StripIterations": 5000, + "StripThresholdFactor": 1.0, + "SnipWidth": 16, + "EstimatePolyOnStrip": True, +} + +# to avoid costly computations when parameters stay the same +_BG_STRIP_OLDY = numpy.array([]) +_BG_STRIP_OLDPARS = [0, 0] +_BG_STRIP_OLDBG = numpy.array([]) + +_BG_SNIP_OLDY = numpy.array([]) +_BG_SNIP_OLDWIDTH = None +_BG_SNIP_OLDBG = numpy.array([]) + + +_BG_OLD_ANCHORS = [] +_BG_OLD_ANCHORS_FLAG = None + +_BG_SMOOTH_OLDWIDTH = None +_BG_SMOOTH_OLDFLAG = None + + +def _convert_anchors_to_indices(x): + """Anchors stored in CONFIG["AnchorsList"] are abscissa. + Convert then to indices (take first index where x >= anchor), + then return the list of indices. + + :param x: Original array of abscissa + :return: List of indices of anchors in x array. + If CONFIG['AnchorsFlag'] is False or None, or if the list + of indices is empty, return None. + """ + # convert anchor X abscissa to index + if CONFIG["AnchorsFlag"] and CONFIG["AnchorsList"] is not None: + anchors_indices = [] + for anchor_x in CONFIG["AnchorsList"]: + if anchor_x <= x[0]: + continue + # take the first index where x > anchor_x + indices = numpy.nonzero(x >= anchor_x)[0] + if len(indices): + anchors_indices.append(min(indices)) + if not len(anchors_indices): + anchors_indices = None + else: + anchors_indices = None + + return anchors_indices + + +def strip_bg(x, y0, width, niter): + """Extract and return the strip bg from y0. + + Use anchors coordinates in CONFIG["AnchorsList"] if flag + CONFIG["AnchorsFlag"] is True. Convert anchors from x coordinate + to array index prior to passing it to silx.math.fit.filters.strip + + :param x: Abscissa array + :param x: Ordinate array (data values at x positions) + :param width: strip width + :param niter: strip niter + """ + global _BG_STRIP_OLDY + global _BG_STRIP_OLDPARS + global _BG_STRIP_OLDBG + global _BG_SMOOTH_OLDWIDTH + global _BG_SMOOTH_OLDFLAG + global _BG_OLD_ANCHORS + global _BG_OLD_ANCHORS_FLAG + + parameters_changed = ( + _BG_STRIP_OLDPARS != [width, niter] + or _BG_SMOOTH_OLDWIDTH != CONFIG["SmoothingWidth"] + or _BG_SMOOTH_OLDFLAG != CONFIG["SmoothingFlag"] + or _BG_OLD_ANCHORS_FLAG != CONFIG["AnchorsFlag"] + or _BG_OLD_ANCHORS != CONFIG["AnchorsList"] + ) + + # same parameters + if not parameters_changed: + # same data + if numpy.array_equal(_BG_STRIP_OLDY, y0): + # same result + return _BG_STRIP_OLDBG + + _BG_STRIP_OLDY = y0 + _BG_STRIP_OLDPARS = [width, niter] + _BG_SMOOTH_OLDWIDTH = CONFIG["SmoothingWidth"] + _BG_SMOOTH_OLDFLAG = CONFIG["SmoothingFlag"] + _BG_OLD_ANCHORS = CONFIG["AnchorsList"] + _BG_OLD_ANCHORS_FLAG = CONFIG["AnchorsFlag"] + + y1 = savitsky_golay(y0, CONFIG["SmoothingWidth"]) if CONFIG["SmoothingFlag"] else y0 + + anchors_indices = _convert_anchors_to_indices(x) + + background = strip( + y1, + w=width, + niterations=niter, + factor=CONFIG["StripThresholdFactor"], + anchors=anchors_indices, + ) + + _BG_STRIP_OLDBG = background + + return background + + +def snip_bg(x, y0, width): + """Compute the snip bg for y0""" + global _BG_SNIP_OLDY + global _BG_SNIP_OLDWIDTH + global _BG_SNIP_OLDBG + global _BG_SMOOTH_OLDWIDTH + global _BG_SMOOTH_OLDFLAG + global _BG_OLD_ANCHORS + global _BG_OLD_ANCHORS_FLAG + + parameters_changed = ( + _BG_SNIP_OLDWIDTH != width + or _BG_SMOOTH_OLDWIDTH != CONFIG["SmoothingWidth"] + or _BG_SMOOTH_OLDFLAG != CONFIG["SmoothingFlag"] + or _BG_OLD_ANCHORS_FLAG != CONFIG["AnchorsFlag"] + or _BG_OLD_ANCHORS != CONFIG["AnchorsList"] + ) + + # same parameters + if not parameters_changed: + # same data + if numpy.sum(_BG_SNIP_OLDY == y0) == len(y0): + # same result + return _BG_SNIP_OLDBG + + _BG_SNIP_OLDY = y0 + _BG_SNIP_OLDWIDTH = width + _BG_SMOOTH_OLDWIDTH = CONFIG["SmoothingWidth"] + _BG_SMOOTH_OLDFLAG = CONFIG["SmoothingFlag"] + _BG_OLD_ANCHORS = CONFIG["AnchorsList"] + _BG_OLD_ANCHORS_FLAG = CONFIG["AnchorsFlag"] + + y1 = savitsky_golay(y0, CONFIG["SmoothingWidth"]) if CONFIG["SmoothingFlag"] else y0 + + anchors_indices = _convert_anchors_to_indices(x) + + if anchors_indices is None or not len(anchors_indices): + anchors_indices = [0, len(y1) - 1] + + background = numpy.zeros_like(y1) + previous_anchor = 0 + for anchor_index in anchors_indices: + if (anchor_index > previous_anchor) and (anchor_index < len(y1)): + background[previous_anchor:anchor_index] = snip1d( + y1[previous_anchor:anchor_index], width + ) + previous_anchor = anchor_index + + if previous_anchor < len(y1): + background[previous_anchor:] = snip1d(y1[previous_anchor:], width) + + _BG_SNIP_OLDBG = background + + return background + + +def estimate_linear(x, y): + """ + Estimate the linear parameters (constant, slope) of a y signal. + + Strip peaks, then perform a linear regression. + """ + bg = strip_bg(x, y, width=CONFIG["StripWidth"], niter=CONFIG["StripIterations"]) + n = float(len(bg)) + Sy = numpy.sum(bg) + Sx = float(numpy.sum(x)) + Sxx = float(numpy.sum(x * x)) + Sxy = float(numpy.sum(x * bg)) + + deno = n * Sxx - (Sx * Sx) + if deno != 0: + bg = (Sxx * Sy - Sx * Sxy) / deno + slope = (n * Sxy - Sx * Sy) / deno + else: + bg = 0.0 + slope = 0.0 + estimated_par = [bg, slope] + # code = 0: FREE + constraints = [[0, 0, 0], [0, 0, 0]] + return estimated_par, constraints + + +def estimate_strip(x, y): + """Estimation function for strip parameters. + + Return parameters as defined in CONFIG dict, + set constraints to FIXED. + """ + estimated_par = [CONFIG["StripWidth"], CONFIG["StripIterations"]] + constraints = numpy.zeros((len(estimated_par), 3), numpy.float64) + # code = 3: FIXED + constraints[0][0] = 3 + constraints[1][0] = 3 + return estimated_par, constraints + + +def estimate_snip(x, y): + """Estimation function for snip parameters. + + Return parameters as defined in CONFIG dict, + set constraints to FIXED. + """ + estimated_par = [CONFIG["SnipWidth"]] + constraints = numpy.zeros((len(estimated_par), 3), numpy.float64) + # code = 3: FIXED + constraints[0][0] = 3 + return estimated_par, constraints + + +def poly(x, y, *pars): + """Order n polynomial. + The order of the polynomial is defined by the number of + coefficients (``*pars``). + + """ + p = numpy.poly1d(pars) + return p(x) + + +def estimate_poly(x, y, deg=2): + """Estimate polynomial coefficients.""" + # extract bg signal with strip, to estimate polynomial on background + if CONFIG["EstimatePolyOnStrip"]: + y = strip_bg(x, y, CONFIG["StripWidth"], CONFIG["StripIterations"]) + pcoeffs = numpy.polyfit(x, y, deg) + cons = numpy.zeros((deg + 1, 3), numpy.float64) + return pcoeffs, cons + + +def estimate_quadratic_poly(x, y): + """Estimate quadratic polynomial coefficients.""" + return estimate_poly(x, y, deg=2) + + +def estimate_cubic_poly(x, y): + """Estimate cubic polynomial coefficients.""" + return estimate_poly(x, y, deg=3) + + +def estimate_quartic_poly(x, y): + """Estimate degree 4 polynomial coefficients.""" + return estimate_poly(x, y, deg=4) + + +def estimate_quintic_poly(x, y): + """Estimate degree 5 polynomial coefficients.""" + return estimate_poly(x, y, deg=5) + + +def configure(**kw): + """Update the CONFIG dict""" + # inspect **kw to find known keys, update them in CONFIG + for key in CONFIG: + if key in kw: + CONFIG[key] = kw[key] + + return CONFIG + + +THEORY = dict( + ( + ( + "No Background", + FitTheory( + description="No background function", + function=lambda x, y0: numpy.zeros_like(x), + parameters=[], + is_background=True, + ), + ), + ( + "Constant", + FitTheory( + description="Constant background", + function=lambda x, y0, c: c * numpy.ones_like(x), + parameters=[ + "Constant", + ], + estimate=lambda x, y: ([min(y)], [[0, 0, 0]]), + is_background=True, + ), + ), + ( + "Linear", + FitTheory( + description="Linear background, parameters 'Constant' and" " 'Slope'", + function=lambda x, y0, a, b: a + b * x, + parameters=["Constant", "Slope"], + estimate=estimate_linear, + configure=configure, + is_background=True, + ), + ), + ( + "Strip", + FitTheory( + description="Compute background using a strip filter\n" + "Parameters 'StripWidth', 'StripIterations'", + function=strip_bg, + parameters=["StripWidth", "StripIterations"], + estimate=estimate_strip, + configure=configure, + is_background=True, + ), + ), + ( + "Snip", + FitTheory( + description="Compute background using a snip filter\n" + "Parameter 'SnipWidth'", + function=snip_bg, + parameters=["SnipWidth"], + estimate=estimate_snip, + configure=configure, + is_background=True, + ), + ), + ( + "Degree 2 Polynomial", + FitTheory( + description="Quadratic polynomial background, Parameters " + "'a', 'b' and 'c'\ny = a*x^2 + b*x +c", + function=poly, + parameters=["a", "b", "c"], + estimate=estimate_quadratic_poly, + configure=configure, + is_background=True, + ), + ), + ( + "Degree 3 Polynomial", + FitTheory( + description="Cubic polynomial background, Parameters " + "'a', 'b', 'c' and 'd'\n" + "y = a*x^3 + b*x^2 + c*x + d", + function=poly, + parameters=["a", "b", "c", "d"], + estimate=estimate_cubic_poly, + configure=configure, + is_background=True, + ), + ), + ( + "Degree 4 Polynomial", + FitTheory( + description="Quartic polynomial background\n" + "y = a*x^4 + b*x^3 + c*x^2 + d*x + e", + function=poly, + parameters=["a", "b", "c", "d", "e"], + estimate=estimate_quartic_poly, + configure=configure, + is_background=True, + ), + ), + ( + "Degree 5 Polynomial", + FitTheory( + description="Quaintic polynomial background\n" + "y = a*x^5 + b*x^4 + c*x^3 + d*x^2 + e*x + f", + function=poly, + parameters=["a", "b", "c", "d", "e", "f"], + estimate=estimate_quintic_poly, + configure=configure, + is_background=True, + ), + ), + ) +) |