path: root/silx/math/fit/test/test_fitmanager.py

# coding: utf-8
# /*##########################################################################
# Copyright (C) 2016-2017 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.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(unittest.TestCase):
    """
    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.float)

        p = [1000, 100., 250,
             255, 650., 45,
             1500, 800.5, 95]
        linear_bg = 2.65 * x + 13
        y = linear_bg + sum_gauss(x, *p)

        # Fitting
        fit = fitmanager.FitManager()
        fit.setdata(x=x, y=y)
        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.float)

        # 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.float)

        # 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.float)

        # 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.float)

            # ('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.float)

    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)


test_cases = (TestFitmanager, TestPolynomials)


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")