# 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. # # ############################################################################*/ """ Nominal tests of the leastsq function. """ import unittest import numpy import sys from silx.test import utils 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 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.float) 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.float) 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.float).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) @utils.test_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) @utils.test_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]) test_cases = (Test_leastsq,) def suite(): loader = unittest.defaultTestLoader test_suite = unittest.TestSuite() for test_class in test_cases: tests = loader.loadTestsFromTestCase(test_class) test_suite.addTests(tests) return test_suite if __name__ == '__main__': unittest.main(defaultTest="suite")