diff options
Diffstat (limited to 'src/silx/math/fit/test')
| -rw-r--r-- | src/silx/math/fit/test/__init__.py | 23 | ||||
| -rw-r--r-- | src/silx/math/fit/test/test_bgtheories.py | 154 | ||||
| -rw-r--r-- | src/silx/math/fit/test/test_filters.py | 122 | ||||
| -rw-r--r-- | src/silx/math/fit/test/test_fit.py | 373 | ||||
| -rw-r--r-- | src/silx/math/fit/test/test_fitmanager.py | 498 | ||||
| -rw-r--r-- | src/silx/math/fit/test/test_functions.py | 259 | ||||
| -rw-r--r-- | src/silx/math/fit/test/test_peaks.py | 132 |
7 files changed, 1561 insertions, 0 deletions
diff --git a/src/silx/math/fit/test/__init__.py b/src/silx/math/fit/test/__init__.py new file mode 100644 index 0000000..745efe3 --- /dev/null +++ b/src/silx/math/fit/test/__init__.py @@ -0,0 +1,23 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016 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. +# +# ############################################################################*/ diff --git a/src/silx/math/fit/test/test_bgtheories.py b/src/silx/math/fit/test/test_bgtheories.py new file mode 100644 index 0000000..6620d38 --- /dev/null +++ b/src/silx/math/fit/test/test_bgtheories.py @@ -0,0 +1,154 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016 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. +# +# ############################################################################*/ +import copy +import unittest +import numpy +import random + +from silx.math.fit import bgtheories +from silx.math.fit.functions import sum_gauss + + +class TestBgTheories(unittest.TestCase): + """ + """ + def setUp(self): + self.x = numpy.arange(100) + self.y = 10 + 0.05 * self.x + sum_gauss(self.x, 10., 45., 15.) + # add a very narrow high amplitude peak to test strip and snip + self.y += sum_gauss(self.x, 100., 75., 2.) + self.narrow_peak_index = list(self.x).index(75) + random.seed() + + def tearDown(self): + pass + + def testTheoriesAttrs(self): + for theory_name in bgtheories.THEORY: + self.assertIsInstance(theory_name, str) + self.assertTrue(hasattr(bgtheories.THEORY[theory_name], + "function")) + self.assertTrue(hasattr(bgtheories.THEORY[theory_name].function, + "__call__")) + # Ensure legacy functions are not renamed accidentally + self.assertTrue( + {"No Background", "Constant", "Linear", "Strip", "Snip"}.issubset( + set(bgtheories.THEORY))) + + def testNoBg(self): + nobgfun = bgtheories.THEORY["No Background"].function + self.assertTrue(numpy.array_equal(nobgfun(self.x, self.y), + numpy.zeros_like(self.x))) + # default estimate + self.assertEqual(bgtheories.THEORY["No Background"].estimate(self.x, self.y), + ([], [])) + + def testConstant(self): + consfun = bgtheories.THEORY["Constant"].function + c = random.random() * 100 + self.assertTrue(numpy.array_equal(consfun(self.x, self.y, c), + c * numpy.ones_like(self.x))) + # default estimate + esti_par, cons = bgtheories.THEORY["Constant"].estimate(self.x, self.y) + self.assertEqual(cons, + [[0, 0, 0]]) + self.assertAlmostEqual(esti_par, + min(self.y)) + + def testLinear(self): + linfun = bgtheories.THEORY["Linear"].function + a = random.random() * 100 + b = random.random() * 100 + self.assertTrue(numpy.array_equal(linfun(self.x, self.y, a, b), + a + b * self.x)) + # default estimate + esti_par, cons = bgtheories.THEORY["Linear"].estimate(self.x, self.y) + + self.assertEqual(cons, + [[0, 0, 0], [0, 0, 0]]) + self.assertAlmostEqual(esti_par[0], 10, places=3) + self.assertAlmostEqual(esti_par[1], 0.05, places=3) + + def testStrip(self): + stripfun = bgtheories.THEORY["Strip"].function + anchors = sorted(random.sample(list(self.x), 4)) + anchors_indices = [list(self.x).index(a) for a in anchors] + + # we really want to strip away the narrow peak + anchors_indices_copy = copy.deepcopy(anchors_indices) + for idx in anchors_indices_copy: + if abs(idx - self.narrow_peak_index) < 5: + anchors_indices.remove(idx) + anchors.remove(self.x[idx]) + + width = 2 + niter = 1000 + bgtheories.THEORY["Strip"].configure(AnchorsList=anchors, AnchorsFlag=True) + + bg = stripfun(self.x, self.y, width, niter) + + # assert peak amplitude has been decreased + self.assertLess(bg[self.narrow_peak_index], + self.y[self.narrow_peak_index]) + + # default estimate + for i in anchors_indices: + self.assertEqual(bg[i], self.y[i]) + + # estimated parameters are equal to the default ones in the config dict + bgtheories.THEORY["Strip"].configure(StripWidth=7, StripIterations=8) + esti_par, cons = bgtheories.THEORY["Strip"].estimate(self.x, self.y) + self.assertTrue(numpy.array_equal(cons, [[3, 0, 0], [3, 0, 0]])) + self.assertEqual(esti_par, [7, 8]) + + def testSnip(self): + snipfun = bgtheories.THEORY["Snip"].function + anchors = sorted(random.sample(list(self.x), 4)) + anchors_indices = [list(self.x).index(a) for a in anchors] + + # we want to strip away the narrow peak, so remove nearby anchors + anchors_indices_copy = copy.deepcopy(anchors_indices) + for idx in anchors_indices_copy: + if abs(idx - self.narrow_peak_index) < 5: + anchors_indices.remove(idx) + anchors.remove(self.x[idx]) + + width = 16 + bgtheories.THEORY["Snip"].configure(AnchorsList=anchors, AnchorsFlag=True) + bg = snipfun(self.x, self.y, width) + + # assert peak amplitude has been decreased + self.assertLess(bg[self.narrow_peak_index], + self.y[self.narrow_peak_index], + "Snip didn't decrease the peak amplitude.") + + # anchored data must remain fixed + for i in anchors_indices: + self.assertEqual(bg[i], self.y[i]) + + # estimated parameters are equal to the default ones in the config dict + bgtheories.THEORY["Snip"].configure(SnipWidth=7) + esti_par, cons = bgtheories.THEORY["Snip"].estimate(self.x, self.y) + self.assertTrue(numpy.array_equal(cons, [[3, 0, 0]])) + self.assertEqual(esti_par, [7]) diff --git a/src/silx/math/fit/test/test_filters.py b/src/silx/math/fit/test/test_filters.py new file mode 100644 index 0000000..8314bdc --- /dev/null +++ b/src/silx/math/fit/test/test_filters.py @@ -0,0 +1,122 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016 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. +# +# ############################################################################*/ +import numpy +import unittest +from silx.math.fit import filters +from silx.math.fit import functions +from silx.test.utils import add_relative_noise + + +class TestSmooth(unittest.TestCase): + """ + Unit tests of smoothing functions. + + Test that the difference between a synthetic curve with 5% added random + noise and the result of smoothing that signal is less than 5%. We compare + the sum of all samples in each curve. + """ + def setUp(self): + x = numpy.arange(5000) + # (height1, center1, fwhm1, beamfwhm...) + slit_params = (50, 500, 200, 100, + 50, 600, 80, 30, + 20, 2000, 150, 150, + 50, 2250, 110, 100, + 40, 3000, 50, 10, + 23, 4980, 250, 20) + + self.y1 = functions.sum_slit(x, *slit_params) + # 5% noise + self.y1 = add_relative_noise(self.y1, 5.) + + # (height1, center1, fwhm1...) + step_params = (50, 500, 200, + 50, 600, 80, + 20, 2000, 150, + 50, 2250, 110, + 40, 3000, 50, + 23, 4980, 250,) + + self.y2 = functions.sum_stepup(x, *step_params) + # 5% noise + self.y2 = add_relative_noise(self.y2, 5.) + + self.y3 = functions.sum_stepdown(x, *step_params) + # 5% noise + self.y3 = add_relative_noise(self.y3, 5.) + + def tearDown(self): + pass + + def testSavitskyGolay(self): + npts = 25 + for y in [self.y1, self.y2, self.y3]: + smoothed_y = filters.savitsky_golay(y, npoints=npts) + + # we added +-5% of random noise. The difference must be much lower + # than 5%. + diff = abs(sum(smoothed_y) - sum(y)) / sum(y) + self.assertLess(diff, 0.05, + "Difference between data with 5%% noise and " + + "smoothed data is > 5%% (%f %%)" % (diff * 100)) + + # Try various smoothing levels + npts += 25 + + def testSmooth1d(self): + """Test the 1D smoothing against the formula + ys[i] = (y[i-1] + 2 * y[i] + y[i+1]) / 4 (for 1 < i < n-1)""" + smoothed_y = filters.smooth1d(self.y1) + + for i in range(1, len(self.y1) - 1): + self.assertAlmostEqual(4 * smoothed_y[i], + self.y1[i-1] + 2 * self.y1[i] + self.y1[i+1]) + + def testSmooth2d(self): + """Test that a 2D smoothing is the same as two successive and + orthogonal 1D smoothings""" + x = numpy.arange(10000) + + noise = 2 * numpy.random.random(10000) - 1 + noise *= 0.05 + y = x * (1 + noise) + + y.shape = (100, 100) + + smoothed_y = filters.smooth2d(y) + + intermediate_smooth = numpy.zeros_like(y) + expected_smooth = numpy.zeros_like(y) + # smooth along first dimension + for i in range(0, y.shape[0]): + intermediate_smooth[i, :] = filters.smooth1d(y[i, :]) + + # smooth along second dimension + for j in range(0, y.shape[1]): + expected_smooth[:, j] = filters.smooth1d(intermediate_smooth[:, j]) + + for i in range(0, y.shape[0]): + for j in range(0, y.shape[1]): + self.assertAlmostEqual(smoothed_y[i, j], + expected_smooth[i, j]) diff --git a/src/silx/math/fit/test/test_fit.py b/src/silx/math/fit/test/test_fit.py new file mode 100644 index 0000000..00f04e2 --- /dev/null +++ b/src/silx/math/fit/test/test_fit.py @@ -0,0 +1,373 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016-2021 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. +# +# ############################################################################*/ +""" +Nominal tests of the leastsq function. +""" + +import unittest + +import numpy +import sys + +from silx.utils import testutils +from silx.math.fit.leastsq import _logger as fitlogger + + +class Test_leastsq(unittest.TestCase): + """ + Unit tests of the leastsq function. + """ + + ndims = None + + def setUp(self): + try: + from silx.math.fit import leastsq + self.instance = leastsq + except ImportError: + self.instance = None + + def myexp(x): + # put a (bad) filter to avoid over/underflows + # with no python looping + with numpy.errstate(invalid='ignore'): + return numpy.exp(x*numpy.less(abs(x), 250)) - \ + 1.0 * numpy.greater_equal(abs(x), 250) + + self.my_exp = myexp + + def gauss(x, *params): + params = numpy.array(params, copy=False, dtype=numpy.float64) + result = params[0] + params[1] * x + for i in range(2, len(params), 3): + p = params[i:(i+3)] + dummy = 2.3548200450309493*(x - p[1])/p[2] + result += p[0] * self.my_exp(-0.5 * dummy * dummy) + return result + + self.gauss = gauss + + def gauss_derivative(x, params, idx): + if idx == 0: + return numpy.ones(len(x), numpy.float64) + if idx == 1: + return x + gaussian_peak = (idx - 2) // 3 + gaussian_parameter = (idx - 2) % 3 + actual_idx = 2 + 3 * gaussian_peak + p = params[actual_idx:(actual_idx+3)] + if gaussian_parameter == 0: + return self.gauss(x, *[0, 0, 1.0, p[1], p[2]]) + if gaussian_parameter == 1: + tmp = self.gauss(x, *[0, 0, p[0], p[1], p[2]]) + tmp *= 2.3548200450309493*(x - p[1])/p[2] + return tmp * 2.3548200450309493/p[2] + if gaussian_parameter == 2: + tmp = self.gauss(x, *[0, 0, p[0], p[1], p[2]]) + tmp *= 2.3548200450309493*(x - p[1])/p[2] + return tmp * 2.3548200450309493*(x - p[1])/(p[2]*p[2]) + + self.gauss_derivative = gauss_derivative + + def tearDown(self): + self.instance = None + self.gauss = None + self.gauss_derivative = None + self.my_exp = None + self.model_function = None + self.model_derivative = None + + def testImport(self): + self.assertTrue(self.instance is not None, + "Cannot import leastsq from silx.math.fit") + + def testUnconstrainedFitNoWeight(self): + parameters_actual = [10.5, 2, 1000.0, 20., 15] + x = numpy.arange(10000.) + y = self.gauss(x, *parameters_actual) + parameters_estimate = [0.0, 1.0, 900.0, 25., 10] + model_function = self.gauss + + fittedpar, cov = self.instance(model_function, x, y, parameters_estimate) + test_condition = numpy.allclose(parameters_actual, fittedpar) + if not test_condition: + msg = "Unsuccessfull fit\n" + for i in range(len(fittedpar)): + msg += "Expected %g obtained %g\n" % (parameters_actual[i], + fittedpar[i]) + self.assertTrue(test_condition, msg) + + def testUnconstrainedFitWeight(self): + parameters_actual = [10.5,2,1000.0,20.,15] + x = numpy.arange(10000.) + y = self.gauss(x, *parameters_actual) + sigma = numpy.sqrt(y) + parameters_estimate = [0.0, 1.0, 900.0, 25., 10] + model_function = self.gauss + + fittedpar, cov = self.instance(model_function, x, y, + parameters_estimate, + sigma=sigma) + test_condition = numpy.allclose(parameters_actual, fittedpar) + if not test_condition: + msg = "Unsuccessfull fit\n" + for i in range(len(fittedpar)): + msg += "Expected %g obtained %g\n" % (parameters_actual[i], + fittedpar[i]) + self.assertTrue(test_condition, msg) + + def testDerivativeFunction(self): + parameters_actual = [10.5, 2, 10000.0, 20., 150, 5000, 900., 300] + x = numpy.arange(10000.) + y = self.gauss(x, *parameters_actual) + delta = numpy.sqrt(numpy.finfo(numpy.float64).eps) + for i in range(len(parameters_actual)): + p = parameters_actual * 1 + if p[i] == 0: + delta_par = delta + else: + delta_par = p[i] * delta + if i > 2: + p[0] = 0.0 + p[1] = 0.0 + p[i] += delta_par + yPlus = self.gauss(x, *p) + p[i] = parameters_actual[i] - delta_par + yMinus = self.gauss(x, *p) + numerical_derivative = (yPlus - yMinus) / (2 * delta_par) + #numerical_derivative = (self.gauss(x, *p) - y) / delta_par + p[i] = parameters_actual[i] + derivative = self.gauss_derivative(x, p, i) + diff = numerical_derivative - derivative + test_condition = numpy.allclose(numerical_derivative, + derivative, atol=5.0e-6) + if not test_condition: + msg = "Error calculating derivative of parameter %d." % i + msg += "\n diff min = %g diff max = %g" % (diff.min(), diff.max()) + self.assertTrue(test_condition, msg) + + def testConstrainedFit(self): + CFREE = 0 + CPOSITIVE = 1 + CQUOTED = 2 + CFIXED = 3 + CFACTOR = 4 + CDELTA = 5 + CSUM = 6 + parameters_actual = [10.5, 2, 10000.0, 20., 150, 5000, 900., 300] + x = numpy.arange(10000.) + y = self.gauss(x, *parameters_actual) + parameters_estimate = [0.0, 1.0, 900.0, 25., 10, 400, 850, 200] + model_function = self.gauss + model_deriv = self.gauss_derivative + constraints_all_free = [[0, 0, 0]] * len(parameters_actual) + constraints_all_positive = [[1, 0, 0]] * len(parameters_actual) + constraints_delta_position = [[0, 0, 0]] * len(parameters_actual) + constraints_delta_position[6] = [CDELTA, 3, 880] + constraints_sum_position = constraints_all_positive * 1 + constraints_sum_position[6] = [CSUM, 3, 920] + constraints_factor = constraints_delta_position * 1 + constraints_factor[2] = [CFACTOR, 5, 2] + constraints_list = [None, + constraints_all_free, + constraints_all_positive, + constraints_delta_position, + constraints_sum_position] + + # for better code coverage, the warning recommending to set full_output + # to True when using constraints should be shown at least once + full_output = True + for index, constraints in enumerate(constraints_list): + if index == 2: + full_output = None + elif index == 3: + full_output = 0 + for model_deriv in [None, self.gauss_derivative]: + for sigma in [None, numpy.sqrt(y)]: + fittedpar, cov = self.instance(model_function, x, y, + parameters_estimate, + sigma=sigma, + constraints=constraints, + model_deriv=model_deriv, + full_output=full_output)[:2] + full_output = True + + test_condition = numpy.allclose(parameters_actual, fittedpar) + if not test_condition: + msg = "Unsuccessfull fit\n" + for i in range(len(fittedpar)): + msg += "Expected %g obtained %g\n" % (parameters_actual[i], + fittedpar[i]) + self.assertTrue(test_condition, msg) + + def testUnconstrainedFitAnalyticalDerivative(self): + parameters_actual = [10.5, 2, 1000.0, 20., 15] + x = numpy.arange(10000.) + y = self.gauss(x, *parameters_actual) + sigma = numpy.sqrt(y) + parameters_estimate = [0.0, 1.0, 900.0, 25., 10] + model_function = self.gauss + model_deriv = self.gauss_derivative + + fittedpar, cov = self.instance(model_function, x, y, + parameters_estimate, + sigma=sigma, + model_deriv=model_deriv) + test_condition = numpy.allclose(parameters_actual, fittedpar) + if not test_condition: + msg = "Unsuccessfull fit\n" + for i in range(len(fittedpar)): + msg += "Expected %g obtained %g\n" % (parameters_actual[i], + fittedpar[i]) + self.assertTrue(test_condition, msg) + + @testutils.validate_logging(fitlogger.name, warning=2) + def testBadlyShapedData(self): + parameters_actual = [10.5, 2, 1000.0, 20., 15] + x = numpy.arange(10000.).reshape(1000, 10) + y = self.gauss(x, *parameters_actual) + sigma = numpy.sqrt(y) + parameters_estimate = [0.0, 1.0, 900.0, 25., 10] + model_function = self.gauss + + for check_finite in [True, False]: + fittedpar, cov = self.instance(model_function, x, y, + parameters_estimate, + sigma=sigma, + check_finite=check_finite) + test_condition = numpy.allclose(parameters_actual, fittedpar) + if not test_condition: + msg = "Unsuccessfull fit\n" + for i in range(len(fittedpar)): + msg += "Expected %g obtained %g\n" % (parameters_actual[i], + fittedpar[i]) + self.assertTrue(test_condition, msg) + + @testutils.validate_logging(fitlogger.name, warning=3) + def testDataWithNaN(self): + parameters_actual = [10.5, 2, 1000.0, 20., 15] + x = numpy.arange(10000.).reshape(1000, 10) + y = self.gauss(x, *parameters_actual) + sigma = numpy.sqrt(y) + parameters_estimate = [0.0, 1.0, 900.0, 25., 10] + model_function = self.gauss + x[500] = numpy.inf + # check default behavior + try: + self.instance(model_function, x, y, + parameters_estimate, + sigma=sigma) + except ValueError: + info = "%s" % sys.exc_info()[1] + self.assertTrue("array must not contain inf" in info) + + # check requested behavior + try: + self.instance(model_function, x, y, + parameters_estimate, + sigma=sigma, + check_finite=True) + except ValueError: + info = "%s" % sys.exc_info()[1] + self.assertTrue("array must not contain inf" in info) + + fittedpar, cov = self.instance(model_function, x, y, + parameters_estimate, + sigma=sigma, + check_finite=False) + test_condition = numpy.allclose(parameters_actual, fittedpar) + if not test_condition: + msg = "Unsuccessfull fit\n" + for i in range(len(fittedpar)): + msg += "Expected %g obtained %g\n" % (parameters_actual[i], + fittedpar[i]) + self.assertTrue(test_condition, msg) + + # testing now with ydata containing NaN + x = numpy.arange(10000.).reshape(1000, 10) + y[500] = numpy.nan + fittedpar, cov = self.instance(model_function, x, y, + parameters_estimate, + sigma=sigma, + check_finite=False) + + test_condition = numpy.allclose(parameters_actual, fittedpar) + if not test_condition: + msg = "Unsuccessfull fit\n" + for i in range(len(fittedpar)): + msg += "Expected %g obtained %g\n" % (parameters_actual[i], + fittedpar[i]) + self.assertTrue(test_condition, msg) + + # testing now with sigma containing NaN + sigma[300] = numpy.nan + fittedpar, cov = self.instance(model_function, x, y, + parameters_estimate, + sigma=sigma, + check_finite=False) + test_condition = numpy.allclose(parameters_actual, fittedpar) + if not test_condition: + msg = "Unsuccessfull fit\n" + for i in range(len(fittedpar)): + msg += "Expected %g obtained %g\n" % (parameters_actual[i], + fittedpar[i]) + self.assertTrue(test_condition, msg) + + def testUncertainties(self): + """Test for validity of uncertainties in returned full-output + dictionary. This is a non-regression test for pull request #197""" + parameters_actual = [10.5, 2, 1000.0, 20., 15, 2001.0, 30.1, 16] + x = numpy.arange(10000.) + y = self.gauss(x, *parameters_actual) + parameters_estimate = [0.0, 1.0, 900.0, 25., 10., 1500., 20., 2.0] + + # test that uncertainties are not 0. + fittedpar, cov, infodict = self.instance(self.gauss, x, y, parameters_estimate, + full_output=True) + uncertainties = infodict["uncertainties"] + self.assertEqual(len(uncertainties), len(parameters_actual)) + self.assertEqual(len(uncertainties), len(fittedpar)) + for uncertainty in uncertainties: + self.assertNotAlmostEqual(uncertainty, 0.) + + # set constraint FIXED for half the parameters. + # This should cause leastsq to return 100% uncertainty. + parameters_estimate = [10.6, 2.1, 1000.1, 20.1, 15.1, 2001.1, 30.2, 16.1] + CFIXED = 3 + CFREE = 0 + constraints = [] + for i in range(len(parameters_estimate)): + if i % 2: + constraints.append([CFIXED, 0, 0]) + else: + constraints.append([CFREE, 0, 0]) + fittedpar, cov, infodict = self.instance(self.gauss, x, y, parameters_estimate, + constraints=constraints, + full_output=True) + uncertainties = infodict["uncertainties"] + for i in range(len(parameters_estimate)): + if i % 2: + # test that all FIXED parameters have 100% uncertainty + self.assertAlmostEqual(uncertainties[i], + parameters_estimate[i]) diff --git a/src/silx/math/fit/test/test_fitmanager.py b/src/silx/math/fit/test/test_fitmanager.py new file mode 100644 index 0000000..4ab56a5 --- /dev/null +++ b/src/silx/math/fit/test/test_fitmanager.py @@ -0,0 +1,498 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016-2020 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. +# +# ############################################################################*/ +""" +Tests for fitmanager module +""" + +import unittest +import numpy +import os.path + +from silx.math.fit import fitmanager +from silx.math.fit import fittheories +from silx.math.fit import bgtheories +from silx.math.fit.fittheory import FitTheory +from silx.math.fit.functions import sum_gauss, sum_stepdown, sum_stepup + +from silx.utils.testutils import ParametricTestCase +from silx.test.utils import temp_dir + +custom_function_definition = """ +import copy +from silx.math.fit.fittheory import FitTheory + +CONFIG = {'d': 1.} + +def myfun(x, a, b, c): + "Model function" + return (a * x**2 + b * x + c) / CONFIG['d'] + +def myesti(x, y): + "Initial parameters for iterative fit (a, b, c) = (1, 1, 1)" + return (1., 1., 1.), ((0, 0, 0), (0, 0, 0), (0, 0, 0)) + +def myconfig(d=1., **kw): + "This function can modify CONFIG" + CONFIG["d"] = d + return CONFIG + +def myderiv(x, parameters, index): + "Custom derivative (does not work, causes singular matrix)" + pars_plus = copy.copy(parameters) + pars_plus[index] *= 1.0001 + + pars_minus = parameters + pars_minus[index] *= copy.copy(0.9999) + + delta_fun = myfun(x, *pars_plus) - myfun(x, *pars_minus) + delta_par = parameters[index] * 0.0001 * 2 + + return delta_fun / delta_par + +THEORY = { + 'my fit theory': + FitTheory(function=myfun, + parameters=('A', 'B', 'C'), + estimate=myesti, + configure=myconfig, + derivative=myderiv) +} + +""" + +old_custom_function_definition = """ +CONFIG = {'d': 1.0} + +def myfun(x, a, b, c): + "Model function" + return (a * x**2 + b * x + c) / CONFIG['d'] + +def myesti(x, y, bg, xscalinq, yscaling): + "Initial parameters for iterative fit (a, b, c) = (1, 1, 1)" + return (1., 1., 1.), ((0, 0, 0), (0, 0, 0), (0, 0, 0)) + +def myconfig(**kw): + "Update or complete CONFIG dictionary" + for key in kw: + CONFIG[key] = kw[key] + return CONFIG + +THEORY = ['my fit theory'] +PARAMETERS = [('A', 'B', 'C')] +FUNCTION = [myfun] +ESTIMATE = [myesti] +CONFIGURE = [myconfig] + +""" + + +def _order_of_magnitude(x): + return numpy.log10(x).round() + + +class TestFitmanager(ParametricTestCase): + """ + Unit tests of multi-peak functions. + """ + def setUp(self): + pass + + def tearDown(self): + pass + + def testFitManager(self): + """Test fit manager on synthetic data using a gaussian function + and a linear background""" + # Create synthetic data with a sum of gaussian functions + x = numpy.arange(1000).astype(numpy.float64) + + p = [1000, 100., 250, + 255, 650., 45, + 1500, 800.5, 95] + linear_bg = 2.65 * x + 13 + y = linear_bg + sum_gauss(x, *p) + + y_with_nans = numpy.array(y) + y_with_nans[::10] = numpy.nan + + x_with_nans = numpy.array(x) + x_with_nans[5::15] = numpy.nan + + tests = { + 'all finite': (x, y), + 'y with NaNs': (x, y_with_nans), + 'x with NaNs': (x_with_nans, y), + } + + for name, (xdata, ydata) in tests.items(): + with self.subTest(name=name): + # Fitting + fit = fitmanager.FitManager() + fit.setdata(x=xdata, y=ydata) + fit.loadtheories(fittheories) + # Use one of the default fit functions + fit.settheory('Gaussians') + fit.setbackground('Linear') + fit.estimate() + fit.runfit() + + # fit.fit_results[] + + # first 2 parameters are related to the linear background + self.assertEqual(fit.fit_results[0]["name"], "Constant") + self.assertAlmostEqual(fit.fit_results[0]["fitresult"], 13) + self.assertEqual(fit.fit_results[1]["name"], "Slope") + self.assertAlmostEqual(fit.fit_results[1]["fitresult"], 2.65) + + for i, param in enumerate(fit.fit_results[2:]): + param_number = i // 3 + 1 + if i % 3 == 0: + self.assertEqual(param["name"], + "Height%d" % param_number) + elif i % 3 == 1: + self.assertEqual(param["name"], + "Position%d" % param_number) + elif i % 3 == 2: + self.assertEqual(param["name"], + "FWHM%d" % param_number) + + self.assertAlmostEqual(param["fitresult"], + p[i]) + self.assertAlmostEqual(_order_of_magnitude(param["estimation"]), + _order_of_magnitude(p[i])) + + def testLoadCustomFitFunction(self): + """Test FitManager using a custom fit function defined in an external + file and imported with FitManager.loadtheories""" + # Create synthetic data with a sum of gaussian functions + x = numpy.arange(100).astype(numpy.float64) + + # a, b, c are the fit parameters + # d is a known scaling parameter that is set using configure() + a, b, c, d = 1.5, 2.5, 3.5, 4.5 + y = (a * x**2 + b * x + c) / d + + # Fitting + fit = fitmanager.FitManager() + fit.setdata(x=x, y=y) + + # Create a temporary function definition file, and import it + with temp_dir() as tmpDir: + tmpfile = os.path.join(tmpDir, 'customfun.py') + # custom_function_definition + fd = open(tmpfile, "w") + fd.write(custom_function_definition) + fd.close() + fit.loadtheories(tmpfile) + tmpfile_pyc = os.path.join(tmpDir, 'customfun.pyc') + if os.path.exists(tmpfile_pyc): + os.unlink(tmpfile_pyc) + os.unlink(tmpfile) + + fit.settheory('my fit theory') + # Test configure + fit.configure(d=4.5) + fit.estimate() + fit.runfit() + + self.assertEqual(fit.fit_results[0]["name"], + "A1") + self.assertAlmostEqual(fit.fit_results[0]["fitresult"], + 1.5) + self.assertEqual(fit.fit_results[1]["name"], + "B1") + self.assertAlmostEqual(fit.fit_results[1]["fitresult"], + 2.5) + self.assertEqual(fit.fit_results[2]["name"], + "C1") + self.assertAlmostEqual(fit.fit_results[2]["fitresult"], + 3.5) + + def testLoadOldCustomFitFunction(self): + """Test FitManager using a custom fit function defined in an external + file and imported with FitManager.loadtheories (legacy PyMca format)""" + # Create synthetic data with a sum of gaussian functions + x = numpy.arange(100).astype(numpy.float64) + + # a, b, c are the fit parameters + # d is a known scaling parameter that is set using configure() + a, b, c, d = 1.5, 2.5, 3.5, 4.5 + y = (a * x**2 + b * x + c) / d + + # Fitting + fit = fitmanager.FitManager() + fit.setdata(x=x, y=y) + + # Create a temporary function definition file, and import it + with temp_dir() as tmpDir: + tmpfile = os.path.join(tmpDir, 'oldcustomfun.py') + # custom_function_definition + fd = open(tmpfile, "w") + fd.write(old_custom_function_definition) + fd.close() + fit.loadtheories(tmpfile) + tmpfile_pyc = os.path.join(tmpDir, 'oldcustomfun.pyc') + if os.path.exists(tmpfile_pyc): + os.unlink(tmpfile_pyc) + os.unlink(tmpfile) + + fit.settheory('my fit theory') + fit.configure(d=4.5) + fit.estimate() + fit.runfit() + + self.assertEqual(fit.fit_results[0]["name"], + "A1") + self.assertAlmostEqual(fit.fit_results[0]["fitresult"], + 1.5) + self.assertEqual(fit.fit_results[1]["name"], + "B1") + self.assertAlmostEqual(fit.fit_results[1]["fitresult"], + 2.5) + self.assertEqual(fit.fit_results[2]["name"], + "C1") + self.assertAlmostEqual(fit.fit_results[2]["fitresult"], + 3.5) + + def testAddTheory(self, estimate=True): + """Test FitManager using a custom fit function imported with + FitManager.addtheory""" + # Create synthetic data with a sum of gaussian functions + x = numpy.arange(100).astype(numpy.float64) + + # a, b, c are the fit parameters + # d is a known scaling parameter that is set using configure() + a, b, c, d = -3.14, 1234.5, 10000, 4.5 + y = (a * x**2 + b * x + c) / d + + # Fitting + fit = fitmanager.FitManager() + fit.setdata(x=x, y=y) + + # Define and add the fit theory + CONFIG = {'d': 1.} + + def myfun(x_, a_, b_, c_): + """"Model function""" + return (a_ * x_**2 + b_ * x_ + c_) / CONFIG['d'] + + def myesti(x_, y_): + """"Initial parameters for iterative fit: + (a, b, c) = (1, 1, 1) + Constraints all set to 0 (FREE)""" + return (1., 1., 1.), ((0, 0, 0), (0, 0, 0), (0, 0, 0)) + + def myconfig(d_=1., **kw): + """This function can modify CONFIG""" + CONFIG["d"] = d_ + return CONFIG + + def myderiv(x_, parameters, index): + """Custom derivative""" + pars_plus = numpy.array(parameters, copy=True) + pars_plus[index] *= 1.001 + + pars_minus = numpy.array(parameters, copy=True) + pars_minus[index] *= 0.999 + + delta_fun = myfun(x_, *pars_plus) - myfun(x_, *pars_minus) + delta_par = parameters[index] * 0.001 * 2 + + return delta_fun / delta_par + + fit.addtheory("polynomial", + FitTheory(function=myfun, + parameters=["A", "B", "C"], + estimate=myesti if estimate else None, + configure=myconfig, + derivative=myderiv)) + + fit.settheory('polynomial') + fit.configure(d_=4.5) + fit.estimate() + params1, sigmas, infodict = fit.runfit() + + self.assertEqual(fit.fit_results[0]["name"], + "A1") + self.assertAlmostEqual(fit.fit_results[0]["fitresult"], + -3.14) + self.assertEqual(fit.fit_results[1]["name"], + "B1") + # params1[1] is the same as fit.fit_results[1]["fitresult"] + self.assertAlmostEqual(params1[1], + 1234.5) + self.assertEqual(fit.fit_results[2]["name"], + "C1") + self.assertAlmostEqual(params1[2], + 10000) + + # change configuration scaling factor and check that the fit returns + # different values + fit.configure(d_=5.) + fit.estimate() + params2, sigmas, infodict = fit.runfit() + for p1, p2 in zip(params1, params2): + self.assertFalse(numpy.array_equal(p1, p2), + "Fit parameters are equal even though the " + + "configuration has been changed") + + def testNoEstimate(self): + """Ensure that the in the absence of the estimation function, + the default estimation function :meth:`FitTheory.default_estimate` + is used.""" + self.testAddTheory(estimate=False) + + def testStep(self): + """Test fit manager on a step function with a more complex estimate + function than the gaussian (convolution filter)""" + for theory_name, theory_fun in (('Step Down', sum_stepdown), + ('Step Up', sum_stepup)): + # Create synthetic data with a sum of gaussian functions + x = numpy.arange(1000).astype(numpy.float64) + + # ('Height', 'Position', 'FWHM') + p = [1000, 439, 250] + + constantbg = 13 + y = theory_fun(x, *p) + constantbg + + # Fitting + fit = fitmanager.FitManager() + fit.setdata(x=x, y=y) + fit.loadtheories(fittheories) + fit.settheory(theory_name) + fit.setbackground('Constant') + + fit.estimate() + + params, sigmas, infodict = fit.runfit() + + # first parameter is the constant background + self.assertAlmostEqual(params[0], 13, places=5) + for i, param in enumerate(params[1:]): + self.assertAlmostEqual(param, p[i], places=5) + self.assertAlmostEqual(_order_of_magnitude(fit.fit_results[i+1]["estimation"]), + _order_of_magnitude(p[i])) + + +def quadratic(x, a, b, c): + return a * x**2 + b * x + c + + +def cubic(x, a, b, c, d): + return a * x**3 + b * x**2 + c * x + d + + +class TestPolynomials(unittest.TestCase): + """Test polynomial fit theories and fit background""" + def setUp(self): + self.x = numpy.arange(100).astype(numpy.float64) + + def testQuadraticBg(self): + gaussian_params = [100, 45, 8] + poly_params = [0.05, -2, 3] + p = numpy.poly1d(poly_params) + + y = p(self.x) + sum_gauss(self.x, *gaussian_params) + + fm = fitmanager.FitManager(self.x, y) + fm.loadbgtheories(bgtheories) + fm.loadtheories(fittheories) + fm.settheory("Gaussians") + fm.setbackground("Degree 2 Polynomial") + esti_params = fm.estimate() + fit_params = fm.runfit()[0] + + for p, pfit in zip(poly_params + gaussian_params, fit_params): + self.assertAlmostEqual(p, + pfit) + + def testCubicBg(self): + gaussian_params = [1000, 45, 8] + poly_params = [0.0005, -0.05, 3, -4] + p = numpy.poly1d(poly_params) + + y = p(self.x) + sum_gauss(self.x, *gaussian_params) + + fm = fitmanager.FitManager(self.x, y) + fm.loadtheories(fittheories) + fm.settheory("Gaussians") + fm.setbackground("Degree 3 Polynomial") + esti_params = fm.estimate() + fit_params = fm.runfit()[0] + + for p, pfit in zip(poly_params + gaussian_params, fit_params): + self.assertAlmostEqual(p, + pfit) + + def testQuarticcBg(self): + gaussian_params = [10000, 69, 25] + poly_params = [5e-10, 0.0005, 0.005, 2, 4] + p = numpy.poly1d(poly_params) + + y = p(self.x) + sum_gauss(self.x, *gaussian_params) + + fm = fitmanager.FitManager(self.x, y) + fm.loadtheories(fittheories) + fm.settheory("Gaussians") + fm.setbackground("Degree 4 Polynomial") + esti_params = fm.estimate() + fit_params = fm.runfit()[0] + + for p, pfit in zip(poly_params + gaussian_params, fit_params): + self.assertAlmostEqual(p, + pfit, + places=5) + + def _testPoly(self, poly_params, theory, places=5): + p = numpy.poly1d(poly_params) + + y = p(self.x) + + fm = fitmanager.FitManager(self.x, y) + fm.loadbgtheories(bgtheories) + fm.loadtheories(fittheories) + fm.settheory(theory) + esti_params = fm.estimate() + fit_params = fm.runfit()[0] + + for p, pfit in zip(poly_params, fit_params): + self.assertAlmostEqual(p, pfit, places=places) + + def testQuadratic(self): + self._testPoly([0.05, -2, 3], + "Degree 2 Polynomial") + + def testCubic(self): + self._testPoly([0.0005, -0.05, 3, -4], + "Degree 3 Polynomial") + + def testQuartic(self): + self._testPoly([1, -2, 3, -4, -5], + "Degree 4 Polynomial") + + def testQuintic(self): + self._testPoly([1, -2, 3, -4, -5, 6], + "Degree 5 Polynomial", + places=4) diff --git a/src/silx/math/fit/test/test_functions.py b/src/silx/math/fit/test/test_functions.py new file mode 100644 index 0000000..7e3ff63 --- /dev/null +++ b/src/silx/math/fit/test/test_functions.py @@ -0,0 +1,259 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016 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. +# +# ############################################################################*/ +""" +Tests for functions module +""" + +import unittest +import numpy +import math + +from silx.math.fit import functions + +__authors__ = ["P. Knobel"] +__license__ = "MIT" +__date__ = "21/07/2016" + +class Test_functions(unittest.TestCase): + """ + Unit tests of multi-peak functions. + """ + def setUp(self): + self.x = numpy.arange(11) + + # height, center, sigma1, sigma2 + (h, c, s1, s2) = (7., 5., 3., 2.1) + self.g_params = { + "height": h, + "center": c, + #"sigma": s, + "fwhm1": 2 * math.sqrt(2 * math.log(2)) * s1, + "fwhm2": 2 * math.sqrt(2 * math.log(2)) * s2, + "area1": h * s1 * math.sqrt(2 * math.pi) + } + # result of `7 * scipy.signal.gaussian(11, 3)` + self.scipy_gaussian = numpy.array( + [1.74546546, 2.87778603, 4.24571462, 5.60516182, 6.62171628, + 7., 6.62171628, 5.60516182, 4.24571462, 2.87778603, + 1.74546546] + ) + + # result of: + # numpy.concatenate((7 * scipy.signal.gaussian(11, 3)[0:5], + # 7 * scipy.signal.gaussian(11, 2.1)[5:11])) + self.scipy_asym_gaussian = numpy.array( + [1.74546546, 2.87778603, 4.24571462, 5.60516182, 6.62171628, + 7., 6.24968751, 4.44773692, 2.52313452, 1.14093853, 0.41124877] + ) + + def tearDown(self): + pass + + def testGauss(self): + """Compare sum_gauss with scipy.signals.gaussian""" + y = functions.sum_gauss(self.x, + self.g_params["height"], + self.g_params["center"], + self.g_params["fwhm1"]) + + for i in range(11): + self.assertAlmostEqual(y[i], self.scipy_gaussian[i]) + + def testAGauss(self): + """Compare sum_agauss with scipy.signals.gaussian""" + y = functions.sum_agauss(self.x, + self.g_params["area1"], + self.g_params["center"], + self.g_params["fwhm1"]) + for i in range(11): + self.assertAlmostEqual(y[i], self.scipy_gaussian[i]) + + def testFastAGauss(self): + """Compare sum_fastagauss with scipy.signals.gaussian + Limit precision to 3 decimal places.""" + y = functions.sum_fastagauss(self.x, + self.g_params["area1"], + self.g_params["center"], + self.g_params["fwhm1"]) + for i in range(11): + self.assertAlmostEqual(y[i], self.scipy_gaussian[i], 3) + + + def testSplitGauss(self): + """Compare sum_splitgauss with scipy.signals.gaussian""" + y = functions.sum_splitgauss(self.x, + self.g_params["height"], + self.g_params["center"], + self.g_params["fwhm1"], + self.g_params["fwhm2"]) + for i in range(11): + self.assertAlmostEqual(y[i], self.scipy_asym_gaussian[i]) + + def testErf(self): + """Compare erf with math.erf""" + # scalars + self.assertAlmostEqual(functions.erf(0.14), math.erf(0.14), places=5) + self.assertAlmostEqual(functions.erf(0), math.erf(0), places=5) + self.assertAlmostEqual(functions.erf(-0.74), math.erf(-0.74), places=5) + + # lists + x = [-5, -2, -1.5, -0.6, 0, 0.1, 2, 3] + erfx = functions.erf(x) + for i in range(len(x)): + self.assertAlmostEqual(erfx[i], + math.erf(x[i]), + places=5) + + # ndarray + x = numpy.array([[0.1, 0.2, 0.3], [0.4, 0.5, 0.6]]) + erfx = functions.erf(x) + for i in range(x.shape[0]): + for j in range(x.shape[1]): + self.assertAlmostEqual(erfx[i, j], + math.erf(x[i, j]), + places=5) + + def testErfc(self): + """Compare erf with math.erf""" + # scalars + self.assertAlmostEqual(functions.erfc(0.14), math.erfc(0.14), places=5) + self.assertAlmostEqual(functions.erfc(0), math.erfc(0), places=5) + self.assertAlmostEqual(functions.erfc(-0.74), math.erfc(-0.74), places=5) + + # lists + x = [-5, -2, -1.5, -0.6, 0, 0.1, 2, 3] + erfcx = functions.erfc(x) + for i in range(len(x)): + self.assertAlmostEqual(erfcx[i], math.erfc(x[i]), places=5) + + # ndarray + x = numpy.array([[0.1, 0.2, 0.3], [0.4, 0.5, 0.6]]) + erfcx = functions.erfc(x) + for i in range(x.shape[0]): + for j in range(x.shape[1]): + self.assertAlmostEqual(erfcx[i, j], math.erfc(x[i, j]), places=5) + + def testAtanStepUp(self): + """Compare atan_stepup with math.atan + + atan_stepup(x, a, b, c) = a * (0.5 + (arctan((x - b) / c) / pi))""" + x0 = numpy.arange(100) / 6.33 + y0 = functions.atan_stepup(x0, 11.1, 22.2, 3.33) + + for x, y in zip(x0, y0): + self.assertAlmostEqual( + 11.1 * (0.5 + math.atan((x - 22.2) / 3.33) / math.pi), + y + ) + + def testStepUp(self): + """sanity check for step up: + + - derivative must be largest around the step center + - max value must be close to height parameter + + """ + x0 = numpy.arange(1000) + center = 444 + height = 1234 + fwhm = 210 + y0 = functions.sum_stepup(x0, height, center, fwhm) + + self.assertLess(max(y0), height) + self.assertAlmostEqual(max(y0), height, places=1) + self.assertAlmostEqual(min(y0), 0, places=1) + + deriv0 = _numerical_derivative(functions.sum_stepup, x0, [height, center, fwhm]) + + # Test center position within +- 1 sample of max derivative + index_max_deriv = numpy.argmax(deriv0) + self.assertLess(abs(index_max_deriv - center), + 1) + + def testStepDown(self): + """sanity check for step down: + + - absolute value of derivative must be largest around the step center + - max value must be close to height parameter + + """ + x0 = numpy.arange(1000) + center = 444 + height = 1234 + fwhm = 210 + y0 = functions.sum_stepdown(x0, height, center, fwhm) + + self.assertLess(max(y0), height) + self.assertAlmostEqual(max(y0), height, places=1) + self.assertAlmostEqual(min(y0), 0, places=1) + + deriv0 = _numerical_derivative(functions.sum_stepdown, x0, [height, center, fwhm]) + + # Test center position within +- 1 sample of max derivative + index_min_deriv = numpy.argmax(-deriv0) + self.assertLess(abs(index_min_deriv - center), + 1) + + def testSlit(self): + """sanity check for slit: + + - absolute value of derivative must be largest around the step center + - max value must be close to height parameter + + """ + x0 = numpy.arange(1000) + center = 444 + height = 1234 + fwhm = 210 + beamfwhm = 30 + y0 = functions.sum_slit(x0, height, center, fwhm, beamfwhm) + + self.assertAlmostEqual(max(y0), height, places=1) + self.assertAlmostEqual(min(y0), 0, places=1) + + deriv0 = _numerical_derivative(functions.sum_slit, x0, [height, center, fwhm, beamfwhm]) + + # Test step up center position (center - fwhm/2) within +- 1 sample of max derivative + index_max_deriv = numpy.argmax(deriv0) + self.assertLess(abs(index_max_deriv - (center - fwhm/2)), + 1) + # Test step down center position (center + fwhm/2) within +- 1 sample of min derivative + index_min_deriv = numpy.argmin(deriv0) + self.assertLess(abs(index_min_deriv - (center + fwhm/2)), + 1) + + +def _numerical_derivative(f, x, params=[], delta_factor=0.0001): + """Compute the numerical derivative of ``f`` for all values of ``x``. + + :param f: function + :param x: Array of evenly spaced abscissa values + :param params: list of additional parameters + :return: Array of derivative values + """ + deltax = (x[1] - x[0]) * delta_factor + y_plus = f(x + deltax, *params) + y_minus = f(x - deltax, *params) + + return (y_plus - y_minus) / (2 * deltax) diff --git a/src/silx/math/fit/test/test_peaks.py b/src/silx/math/fit/test/test_peaks.py new file mode 100644 index 0000000..495c70d --- /dev/null +++ b/src/silx/math/fit/test/test_peaks.py @@ -0,0 +1,132 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016 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. +# +# ############################################################################*/ +""" +Tests for peaks module +""" + +import unittest +import numpy +import math + +from silx.math.fit import functions +from silx.math.fit import peaks + +class Test_peak_search(unittest.TestCase): + """ + Unit tests of peak_search on various types of multi-peak functions. + """ + def setUp(self): + self.x = numpy.arange(5000) + # (height1, center1, fwhm1, ...) + self.h_c_fwhm = (50, 500, 100, + 50, 600, 80, + 20, 2000, 100, + 50, 2250, 110, + 40, 3000, 99, + 23, 4980, 80) + # (height1, center1, fwhm1, eta1 ...) + self.h_c_fwhm_eta = (50, 500, 100, 0.4, + 50, 600, 80, 0.5, + 20, 2000, 100, 0.6, + 50, 2250, 110, 0.7, + 40, 3000, 99, 0.8, + 23, 4980, 80, 0.3,) + # (height1, center1, fwhm11, fwhm21, ...) + self.h_c_fwhm_fwhm = (50, 500, 100, 85, + 50, 600, 80, 110, + 20, 2000, 100, 100, + 50, 2250, 110, 99, + 40, 3000, 99, 110, + 23, 4980, 80, 80,) + # (height1, center1, fwhm11, fwhm21, eta1 ...) + self.h_c_fwhm_fwhm_eta = (50, 500, 100, 85, 0.4, + 50, 600, 80, 110, 0.5, + 20, 2000, 100, 100, 0.6, + 50, 2250, 110, 99, 0.7, + 40, 3000, 99, 110, 0.8, + 23, 4980, 80, 80, 0.3,) + # (area1, center1, fwhm1, ...) + self.a_c_fwhm = (2550, 500, 100, + 2000, 600, 80, + 500, 2000, 100, + 4000, 2250, 110, + 2300, 3000, 99, + 3333, 4980, 80) + # (area1, center1, fwhm1, eta1 ...) + self.a_c_fwhm_eta = (500, 500, 100, 0.4, + 500, 600, 80, 0.5, + 200, 2000, 100, 0.6, + 500, 2250, 110, 0.7, + 400, 3000, 99, 0.8, + 230, 4980, 80, 0.3,) + # (area, position, fwhm, st_area_r, st_slope_r, lt_area_r, lt_slope_r, step_height_r) + self.hypermet_params = (1000, 500, 200, 0.2, 100, 0.3, 100, 0.05, + 1000, 1000, 200, 0.2, 100, 0.3, 100, 0.05, + 1000, 2000, 200, 0.2, 100, 0.3, 100, 0.05, + 1000, 2350, 200, 0.2, 100, 0.3, 100, 0.05, + 1000, 3000, 200, 0.2, 100, 0.3, 100, 0.05, + 1000, 4900, 200, 0.2, 100, 0.3, 100, 0.05,) + + + def tearDown(self): + pass + + def get_peaks(self, function, params): + """ + + :param function: Multi-peak function + :param params: Parameter for this function + :return: list of (peak, relevance) tuples + """ + y = function(self.x, *params) + return peaks.peak_search(y=y, fwhm=100, relevance_info=True) + + def testPeakSearch_various_functions(self): + """Run peak search on a variety of synthetic functions, and + check that result falls within +-25 samples of the actual peak + (reasonable delta considering a fwhm of ~100 samples) and effects + of overlapping peaks).""" + f_p = ((functions.sum_gauss, self.h_c_fwhm ), + (functions.sum_lorentz, self.h_c_fwhm), + (functions.sum_pvoigt, self.h_c_fwhm_eta), + (functions.sum_splitgauss, self.h_c_fwhm_fwhm), + (functions.sum_splitlorentz, self.h_c_fwhm_fwhm), + (functions.sum_splitpvoigt, self.h_c_fwhm_fwhm_eta), + (functions.sum_agauss, self.a_c_fwhm), + (functions.sum_fastagauss, self.a_c_fwhm), + (functions.sum_alorentz, self.a_c_fwhm), + (functions.sum_apvoigt, self.a_c_fwhm_eta), + (functions.sum_ahypermet, self.hypermet_params), + (functions.sum_fastahypermet, self.hypermet_params),) + + for function, params in f_p: + peaks = self.get_peaks(function, params) + + self.assertEqual(len(peaks), 6, + "Wrong number of peaks detected") + + for i in range(6): + theoretical_peak_index = params[i*(len(params)//6) + 1] + found_peak_index = peaks[i][0] + self.assertLess(abs(found_peak_index - theoretical_peak_index), 25) |