
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples/example_fit_jacobian_benchmark.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_example_fit_jacobian_benchmark.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_example_fit_jacobian_benchmark.py:


Benchmarks of methods with and without computing the Jacobian analytically
==========================================================================

Providing a function that calculates the Jacobian matrix analytically can
reduce the time spent finding a solution. The results from benchmarks comparing
two methods (``leastsq`` and ``least_squares``) with and without a function to
calculate the Jacobian matrix analytically are presented below.

First we define the model function, the residual function, and the appropriate
Jacobian functions:

.. GENERATED FROM PYTHON SOURCE LINES 13-63

.. code-block:: Python

    from timeit import timeit
    from types import SimpleNamespace

    import matplotlib.pyplot as plt
    import numpy as np

    from lmfit import Parameters, minimize

    NUM_JACOBIAN_CALLS = 0


    def func(var, x):
        return var[0] * np.exp(-var[1]*x) + var[2]


    def residual(pars, x, data):
        a, b, c = pars['a'], pars['b'], pars['c']
        model = func((a, b, c), x)
        return model - data


    def dfunc(pars, x, data):
        global NUM_JACOBIAN_CALLS
        NUM_JACOBIAN_CALLS += 1

        a, b = pars['a'], pars['b']
        v = np.exp(-b*x)
        return np.array([v, -a*x*v, np.ones(len(x))])


    def jacfunc(pars, x, data):
        global NUM_JACOBIAN_CALLS
        NUM_JACOBIAN_CALLS += 1

        a, b = pars['a'], pars['b']
        v = np.exp(-b*x)
        jac = np.ones((len(x), 3), dtype=np.float64)
        jac[:, 0] = v
        jac[:, 1] = -a * x * v
        return jac


    a, b, c = 2.5, 1.3, 0.8

    x = np.linspace(0, 4, 50)
    y = func([a, b, c], x)

    data = y + 0.15*np.random.RandomState(seed=2021).normal(size=x.size)









.. GENERATED FROM PYTHON SOURCE LINES 64-67

Then we define the different cases to benchmark (i.e., different methods with
and without a function to calculate the Jacobian analytically) and the number
of repetitions per case:

.. GENERATED FROM PYTHON SOURCE LINES 67-114

.. code-block:: Python

    cases = (
        dict(
            method='leastsq',
        ),
        dict(
            method='leastsq',
            Dfun=dfunc,
            col_deriv=1,
        ),
        dict(
            method='least_squares',
        ),
        dict(
            method='least_squares',
            jac=jacfunc,
        ),
    )

    num_repeats = 100
    results = []

    for kwargs in cases:
        params = Parameters()
        params.add('a', value=10)
        params.add('b', value=10)
        params.add('c', value=10)

        wrapper = lambda: minimize(
            residual,
            params,
            args=(x,),
            kws={'data': data},
            **kwargs,
        )
        time = timeit(wrapper, number=num_repeats) / num_repeats

        NUM_JACOBIAN_CALLS = 0
        fit = wrapper()

        results.append(SimpleNamespace(
            time=time,
            num_jacobian_calls=NUM_JACOBIAN_CALLS,
            fit=fit,
            kwargs=kwargs,
        ))









.. GENERATED FROM PYTHON SOURCE LINES 115-116

Finally, we present the results:

.. GENERATED FROM PYTHON SOURCE LINES 116-179

.. code-block:: Python

    labels = []

    for result in results:
        label = result.kwargs['method']
        if result.num_jacobian_calls > 0:
            label += ' with Jac.'

        labels.append(label)

    label_width = max(map(len, labels))
    lines = [
        '| '
        + ' | '.join([
            'Method'.ljust(label_width),
            'Avg. time (ms)',
            '# func. (+ Jac.) calls',
            'Chi-squared',
            'a'.ljust(5),
            'b'.ljust(5),
            'c'.ljust(6),
        ])
        + '|'
    ]

    print(f'The "true" parameters are: a = {a:.3f}, b = {b:.3f}, c = {c:.3f}\n')
    fig, ax = plt.subplots()
    ax.plot(x, data, marker='.', linestyle='none', label='data')

    for (result, label) in zip(results, labels):
        linestyle = '-'
        if result.num_jacobian_calls > 0:
            linestyle = '--'

        a = result.fit.params['a'].value
        b = result.fit.params['b'].value
        c = result.fit.params['c'].value
        y = func([a, b, c], x)
        ax.plot(x, y, label=label, alpha=0.5, linestyle=linestyle)

        columns = [
            label.ljust(label_width),
            f'{result.time * 1000:.2f}'.ljust(14),
            (
                f'{result.fit.nfev}'
                + (
                    f' (+{result.num_jacobian_calls})'
                    if result.num_jacobian_calls > 0 else
                    ''
                )
            ).ljust(22),
            f'{result.fit.chisqr:.3f}'.ljust(11),
            f'{a:.3f}'.ljust(5, '0'),
            f'{b:.3f}'.ljust(5, '0'),
            f'{c:.3f}'.ljust(5, '0'),
        ]
        lines.append('| ' + ' | '.join(columns) + ' |')

    lines.insert(1, '|-' + '-|-'.join('-' * len(col) for col in columns) + '-|')
    print('\n'.join(lines))

    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.legend()



.. image-sg:: /examples/images/sphx_glr_example_fit_jacobian_benchmark_001.png
   :alt: example fit jacobian benchmark
   :srcset: /examples/images/sphx_glr_example_fit_jacobian_benchmark_001.png, /examples/images/sphx_glr_example_fit_jacobian_benchmark_001_3_00x.png 3.00x
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    The "true" parameters are: a = 2.500, b = 1.300, c = 0.800

    | Method                  | Avg. time (ms) | # func. (+ Jac.) calls | Chi-squared | a     | b     | c     |
    |-------------------------|----------------|------------------------|-------------|-------|-------|-------|
    | leastsq                 | 1.53           | 39                     | 1.092       | 2.564 | 1.359 | 0.824 |
    | leastsq with Jac.       | 1.30           | 12 (+10)               | 1.092       | 2.564 | 1.359 | 0.824 |
    | least_squares           | 3.53           | 38                     | 1.092       | 2.564 | 1.359 | 0.824 |
    | least_squares with Jac. | 2.34           | 11 (+9)                | 1.092       | 2.564 | 1.359 | 0.824 |





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.197 seconds)


.. _sphx_glr_download_examples_example_fit_jacobian_benchmark.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: example_fit_jacobian_benchmark.ipynb <example_fit_jacobian_benchmark.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: example_fit_jacobian_benchmark.py <example_fit_jacobian_benchmark.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: example_fit_jacobian_benchmark.zip <example_fit_jacobian_benchmark.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
