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

.. only:: html

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

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

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

.. _sphx_glr_auto_examples_plot_03_custom_fitting.py:


Custom fitting example
======================

This is a short example for fitting data, which is not covered by the fitting widget.
It relies on calling lmfit manually. Therefore one needs to provide a fitting function, which calculates the numerical residual between measurement and model.
This notebook is about fitting multiple datasets from different angles of incidents simultaneously, but can be adapted to a wide range of different use cases.

.. GENERATED FROM PYTHON SOURCE LINES 12-20

.. code-block:: Python

    import numpy as np
    import elli
    from elli.fitting import ParamsHist
    from lmfit import minimize, fit_report
    import matplotlib.pyplot as plt










.. GENERATED FROM PYTHON SOURCE LINES 22-24

Data import
-----------

.. GENERATED FROM PYTHON SOURCE LINES 24-30

.. code-block:: Python

    data = elli.read_nexus_psi_delta("SiO2onSi.ellips.nxs").loc[
        (slice(None), slice(210, 800)), :
    ]
    lbda = data.loc[50].index.get_level_values("Wavelength").to_numpy()
    data






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th></th>
          <th>Ψ</th>
          <th>Δ</th>
        </tr>
        <tr>
          <th>Angle of Incidence</th>
          <th>Wavelength</th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th rowspan="5" valign="top">50</th>
          <th>210.0</th>
          <td>40.028156</td>
          <td>146.370560</td>
        </tr>
        <tr>
          <th>211.0</th>
          <td>40.012478</td>
          <td>146.607407</td>
        </tr>
        <tr>
          <th>212.0</th>
          <td>40.007420</td>
          <td>146.799438</td>
        </tr>
        <tr>
          <th>213.0</th>
          <td>40.005928</td>
          <td>147.055832</td>
        </tr>
        <tr>
          <th>214.0</th>
          <td>40.001190</td>
          <td>147.294281</td>
        </tr>
        <tr>
          <th>...</th>
          <th>...</th>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th rowspan="5" valign="top">70</th>
          <th>796.0</th>
          <td>9.051100</td>
          <td>174.423874</td>
        </tr>
        <tr>
          <th>797.0</th>
          <td>9.054414</td>
          <td>174.494553</td>
        </tr>
        <tr>
          <th>798.0</th>
          <td>9.044289</td>
          <td>174.405136</td>
        </tr>
        <tr>
          <th>799.0</th>
          <td>9.048944</td>
          <td>174.376541</td>
        </tr>
        <tr>
          <th>800.0</th>
          <td>9.021362</td>
          <td>174.418228</td>
        </tr>
      </tbody>
    </table>
    <p>1773 rows × 2 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 31-33

Setting up invariant materials and fitting parameters
-----------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 33-46

.. code-block:: Python

    rii_db = elli.db.RII()
    Si = rii_db.get_mat("Si", "Aspnes")

    params = ParamsHist()
    params.add("SiO2_n0", value=1.452, min=-100, max=100, vary=True)
    params.add("SiO2_n1", value=36.0, min=-40000, max=40000, vary=True)
    params.add("SiO2_n2", value=0, min=-40000, max=40000, vary=False)
    params.add("SiO2_k0", value=0, min=-100, max=100, vary=False)
    params.add("SiO2_k1", value=0, min=-40000, max=40000, vary=False)
    params.add("SiO2_k2", value=0, min=-40000, max=40000, vary=False)
    params.add("SiO2_d", value=20, min=0, max=40000, vary=True)









.. GENERATED FROM PYTHON SOURCE LINES 47-50

Model helper function
---------------------
This model function is not strictly needed, but simplifies the fit function, as the model only needs to be defined once.

.. GENERATED FROM PYTHON SOURCE LINES 50-69

.. code-block:: Python

    def model(lbda, angle, params):
        SiO2 = elli.Cauchy(
            params["SiO2_n0"],
            params["SiO2_n1"],
            params["SiO2_n2"],
            params["SiO2_k0"],
            params["SiO2_k1"],
            params["SiO2_k2"],
        ).get_mat()

        structure = elli.Structure(
            elli.AIR,
            [elli.Layer(SiO2, params["SiO2_d"])],
            Si,
        )

        return structure.evaluate(lbda, angle, solver=elli.Solver2x2)









.. GENERATED FROM PYTHON SOURCE LINES 70-74

Defining the fit function
-------------------------
The fit function follows the protocol defined by the lmfit package and needs the parameters dictionary as first argument.
It has to return a residual value, which will be minimized. Here psi and delta are used to calculate the residual, but could be changed to transmission or reflection data.

.. GENERATED FROM PYTHON SOURCE LINES 74-91

.. code-block:: Python



    def fit_function(params, lbda, data):
        residual = []

        for phi_i in [50, 60, 70]:
            model_result = model(lbda, phi_i, params)

            resid_psi = data.loc[(phi_i, "Ψ")].to_numpy() - model_result.psi
            resid_delta = data.loc[(phi_i, "Δ")].to_numpy() - model_result.delta

            residual.append(resid_psi)
            residual.append(resid_delta)

        return np.concatenate(residual)









.. GENERATED FROM PYTHON SOURCE LINES 92-96

Running the fit
---------------
The fitting is performed by calling the minimize function with the fit_function and the needed arguments.
It is possible to change the underlying algorithm by providing the method kwarg.

.. GENERATED FROM PYTHON SOURCE LINES 96-100

.. code-block:: Python


    out = minimize(fit_function, params, args=(lbda, data), method="leastsq")
    print(fit_report(out))





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

 .. code-block:: none

    [[Fit Statistics]]
        # fitting method   = leastsq
        # function evals   = 119
        # data points      = 3546
        # variables        = 3
        chi-square         = 157.748841
        reduced chi-square = 0.04452409
        Akaike info crit   = -11031.1779
        Bayesian info crit = -11012.6572
    [[Variables]]
        SiO2_n0:  1.63549687 +/- 0.03012997 (1.84%) (init = 1.452)
        SiO2_n1:  27.6408772 +/- 8.75248306 (31.66%) (init = 36)
        SiO2_n2:  0 (fixed)
        SiO2_k0:  0 (fixed)
        SiO2_k1:  0 (fixed)
        SiO2_k2:  0 (fixed)
        SiO2_d:   1.75831329 +/- 0.02404272 (1.37%) (init = 20)
    [[Correlations]] (unreported correlations are < 0.100)
        C(SiO2_n0, SiO2_d)  = -0.9916
        C(SiO2_n0, SiO2_n1) = -0.9596
        C(SiO2_n1, SiO2_d)  = +0.9259




.. GENERATED FROM PYTHON SOURCE LINES 101-103

Plotting the results
----------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 103-133

.. code-block:: Python


    fit_50 = model(lbda, 50, out.params)
    fit_60 = model(lbda, 60, out.params)
    fit_70 = model(lbda, 70, out.params)

    fig = plt.figure(dpi=100)
    ax = fig.add_subplot(1, 1, 1)
    ax.scatter(lbda, data.loc[(50, "Ψ")], s=20, alpha=0.1, label="50° Measurement")
    ax.scatter(lbda, data.loc[(60, "Ψ")], s=20, alpha=0.1, label="Psi 60° Measurement")
    ax.scatter(lbda, data.loc[(70, "Ψ")], s=20, alpha=0.1, label="Psi 70° Measurement")
    (psi50,) = ax.plot(lbda, fit_50.psi, c="tab:blue", label="Psi 50°")
    (psi60,) = ax.plot(lbda, fit_60.psi, c="tab:orange", label="Psi 60°")
    (psi70,) = ax.plot(lbda, fit_70.psi, c="tab:green", label="Psi 70°")
    ax.set_xlabel("wavelenth / nm")
    ax.set_ylabel("psi / degree")
    ax.legend(handles=[psi50, psi60, psi70], loc="lower left")
    fig.canvas.draw()

    fig = plt.figure(dpi=100)
    ax = fig.add_subplot(1, 1, 1)
    ax.scatter(lbda, data.loc[(50, "Δ")], s=20, alpha=0.1, label="Delta 50° Measurement")
    ax.scatter(lbda, data.loc[(60, "Δ")], s=20, alpha=0.1, label="Delta 60° Measurement")
    ax.scatter(lbda, data.loc[(70, "Δ")], s=20, alpha=0.1, label="Delta 70° Measurement")
    (delta50,) = ax.plot(lbda, fit_50.delta, c="tab:blue", label="Delta 50°")
    (delta60,) = ax.plot(lbda, fit_60.delta, c="tab:orange", label="Delta 60°")
    (delta70,) = ax.plot(lbda, fit_70.delta, c="tab:green", label="Delta 70°")
    ax.set_xlabel("wavelenth / nm")
    ax.set_ylabel("delta / degree")
    ax.legend(handles=[delta50, delta60, delta70], loc="lower right")
    fig.canvas.draw()



.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /auto_examples/images/sphx_glr_plot_03_custom_fitting_001.png
         :alt: plot 03 custom fitting
         :srcset: /auto_examples/images/sphx_glr_plot_03_custom_fitting_001.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_examples/images/sphx_glr_plot_03_custom_fitting_002.png
         :alt: plot 03 custom fitting
         :srcset: /auto_examples/images/sphx_glr_plot_03_custom_fitting_002.png
         :class: sphx-glr-multi-img






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

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


.. _sphx_glr_download_auto_examples_plot_03_custom_fitting.py:

.. only:: html

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

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

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

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

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

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

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