Add custom fitting example

docs/requirements.txt

@@ -17,8 +17,6 @@ sphinx-plotly-directive
 sphinxcontrib-mermaid 
 matplotlib
 h5py
-importlib-resources
-h5py
 pyyaml
 importlib-resources
 rapidfuzz

examples/Basic Usage/Basic Usage.ipynb

@@ -185,7 +185,7 @@
     "## Building the model\n",
     "\n",
     "For simple parameter estimation,\n",
-    "the fit decorator (**@fit**) in conjuction with the model definition is used.\n",
+    "the fit decorator (**@fit**) in conjunction with the model definition is used.\n",
     "The fitting decorator takes a pandas dataframe containing\n",
     "the psi/delta measurement data (**psi_delta**) and the model parameters (**params**) as an input.\n",
     "It then passes the wavelength from measurement dataframe (**lbda**)\n",
@@ -204,7 +204,7 @@
     "This is done by calling the :code:`elli.IsotropicMaterial(...)` function\n",
     "with a dispersion model as a parameter\n",
     "or simply calling :code:`.get_mat()` on a dispersion model.\n",
-    "These two approaches are equivalent.From these materials the layer is build,\n",
+    "These two approaches are equivalent. From these materials the layer is build,\n",
     "which only consists of the SiO2 layer in this example.\n",
     "The final structure consists of an incoming half-space,\n",
     "the layers and an outgoing half space. Specifically,\n",
@@ -223,7 +223,7 @@
     "Executing the cell below in a jupyter notebook displays a comparison of the simulated Ψ / Δ values\n",
     "at the current parameter values with their measured counterparts.\n",
     "Additionally, input fields for each model parameter are shown.\n",
-    "You may change the parameters and the calcualted data will change accordingly.\n",
+    "You may change the parameters and the calculated data will change accordingly.\n",
     "For clarification the modeled data is shown with `_calc` postfix in the legend.\n",
     "\n"
    ]

examples/gallery/plot_01_basic_usage.py

@@ -36,7 +36,7 @@
 # Prior to defining our model, we have to set the parameters we want to use.
 # We're going to use a :ref:`Cauchy model <cauchy>` for SiO2 and load the Si values from
 # `literature values <https://refractiveindex.info/?shelf=main&book=Si&page=Aspnes>`_.
-# The parameter names can be choosen freely,
+# The parameter names can be chosen freely,
 # but you have to use the exact same name in the later model definition.
 # The package uses lmfit as fitting tool and you may refer to their
 # `documentation <https://lmfit.github.io/lmfit-py/parameters.html#lmfit.parameter.Parameters.add>`_
@@ -71,7 +71,7 @@
 # ------------------------
 #
 # For simple parameter estimation,
-# the fit decorator (**@fit**) in conjuction with the model definition is used.
+# the fit decorator (**@fit**) in conjunction with the model definition is used.
 # The fitting decorator takes a pandas dataframe containing
 # the psi/delta measurement data (**psi_delta**) and the model parameters (**params**) as an input.
 # It then passes the wavelength from measurement dataframe (**lbda**)
@@ -90,7 +90,7 @@
 # This is done by calling the :code:`elli.IsotropicMaterial(...)` function
 # with a dispersion model as a parameter
 # or simply calling :code:`.get_mat()` on a dispersion model.
-# These two approaches are equivalent.From these materials the layer is build,
+# These two approaches are equivalent. From these materials the layer is build,
 # which only consists of the SiO2 layer in this example.
 # The final structure consists of an incoming half-space,
 # the layers and an outgoing half space. Specifically,
@@ -109,7 +109,7 @@
 # Executing the cell below in a jupyter notebook displays a comparison of the simulated Ψ / Δ values
 # at the current parameter values with their measured counterparts.
 # Additionally, input fields for each model parameter are shown.
-# You may change the parameters and the calcualted data will change accordingly.
+# You may change the parameters and the calculated data will change accordingly.
 # For clarification the modeled data is shown with `_calc` postfix in the legend.
 @fit(psi_delta, params)
 def model(lbda, params):

examples/gallery/plot_03_custom_fitting.py

@@ -0,0 +1,132 @@
+"""
+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.
+
+"""
+
+# %%
+import numpy as np
+import elli
+from elli.fitting import ParamsHist
+from lmfit import minimize, fit_report
+import matplotlib.pyplot as plt
+
+# sphinx_gallery_thumbnail_path = '_static/custom_fitting.png'
+
+
+# %%
+# Data import
+# -----------
+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
+
+# %%
+# Setting up invariant materials and fitting parameters
+# -----------------------------------------------------
+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)
+
+
+# %%
+# Model helper function
+# ---------------------
+# This model function is not strictly needed, but simplifies the fit function, as the model only needs to be defined once.
+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)
+
+
+# %%
+# 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.
+
+
+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)
+
+
+# %%
+# 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.
+
+out = minimize(fit_function, params, args=(lbda, data), method="leastsq")
+print(fit_report(out))
+
+# %%
+# Plotting the results
+# ----------------------------------------------
+
+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()