<jupyter_start><jupyter_text>Disclaimer: The main objective of this jupyter notebook is to show how to fit simple data by    - defining a fitting model- creating the reference data to which the model will be fitted to. In standard cases, this step should be replaced by loading your experimental data- setting and running the fit   - extracting and displaying information about the resultsThe syntax to fit data is minimizer-dependent. Here we focus on a fitting routing provided by the `scipy` python package: scipy.optimize.curve_fit.       Please also note that the initial guessed parameters might not be optimal, resulting in a poor fit of the reference data.Table of Contents- Importing libraries- Definition of the fitting model- Plot of the fitting model- Creating reference data- Setting and fitting- Displaying the results  - Exercises (Top)Importing librariesPlease refer to the notebooks about `numpy` and `matplotlib` for more details.<jupyter_code># import python modules for plotting, fitting
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# for interactive widgets
import ipywidgets as ipw
%matplotlib notebook<jupyter_output><empty_output><jupyter_text>(Top)Definition of the fitting modelHere we define a simple lorentzian function.  Please note that there are other options, such as importing [`lmfit.LorentzianModel`](https://lmfit.github.io/lmfit-py/builtin_models.htmllmfit.models.LorentzianModel),  [`QENSmodels.lorentzian`](https://github.com/QENSlibrary/QENSmodels/blob/master/QENSmodels/lorentzian.py) or [`pyspec.fitfuncs.lor`](https://pythonhosted.org/pyspec/fit.html).<jupyter_code>def lorentzian(xx, scale=1.0, center=1.0, hwhm=3.0):
    if hwhm == 0:
        raise ValueError('hwhm of the lorentzian is equal to zero.')
    return scale * hwhm / ((xx-center)**2 + hwhm**2) / np.pi<jupyter_output><empty_output><jupyter_text>(Top)Plot of the fitting modelThe widget below shows how the variations of the lorentzian's parameters, *Scale*, *Center* and *FWHM*,  influence its shape.<jupyter_code># define values of the variable of the functions to be used in this tutorial
xx = np.linspace(-10, 10, 500)

initial_params = {"scale": 5.5,
                  "center": 5.0,
                  "hwhm": 3.0}

fig1, ax1 = plt.subplots()
lines = ax1.plot(xx, lorentzian(xx, *list(initial_params.values())))
ax1.set_ylabel('lorentzian(x,{},{},{})'.format(*list(initial_params.values())))
ax1.set_xlabel('x')
ax1.grid()

# define slider to interactively modify the parameters
sliders = {key: ipw.FloatSlider(min=1.0, max=10.0, value=value, description=key)
           for key, value in initial_params.items()}
sliders["center"].min = -sliders["center"].max

# define function to be plotted
def interactive_plot(change):
    params = [sliders[key].value for key in sliders]
    lines[0].set_ydata(lorentzian(xx, *params))
    ax1.set_ylabel('lorentzian(x,{},{},{})'.format(*params))

# add observers to the sliders to update the plot
for sl in sliders.values():
    sl.observe(interactive_plot, names="value")

# Define function to reset all parameters' values to the initial ones
def reset_values(b):
    """Reset the interactive plots to inital values."""
    for key, value in initial_params.items():
        sliders[key].value = value

# Define reset button and occurring action when clicking on it
reset_button = ipw.Button(description = "Reset")
reset_button.on_click(reset_values)

# display the interactive plot
display(ipw.VBox(list(sliders.values())), reset_button)<jupyter_output><empty_output><jupyter_text>(Top)Creating reference data**Input:** the reference data for this simple example correspond to a Lorentzian with added noise.<jupyter_code># Creation of reference data
# lorentzian with noise added as background and modulation in its amplitude
added_noise = 0.5*np.random.normal(0, 1, 500)
lorentzian_noisy = lorentzian(xx, 0.89, -0.025, 0.45)*(1. + 0.1*added_noise) + 0.01*added_noise

# Display reference data
fig2, ax2 = plt.subplots()
lines = ax2.plot(xx, lorentzian_noisy, label='reference data')
ax2.set_xlabel('x')
ax2.grid()
ax2.legend()<jupyter_output><empty_output><jupyter_text>(Top)Setting and fittingThe fit is performed using `scipy.optimize.curve_fit`.  The example is based on implementations from https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html<jupyter_code># define set of initial values for the parameters to be refined 
initial_parameters_values = [1, 0.2, 0.5]

# plot initial model and reference data
fig3, ax3 = plt.subplots()
ax3.plot(xx, lorentzian_noisy, label='reference data')
ax3.plot(xx, lorentzian(xx, *initial_parameters_values), '.', label='model with initial guess')
ax3.set_xlabel('x')
ax3.grid()
ax3.legend()<jupyter_output><empty_output><jupyter_text>`scipy.optimize_curve_fit` offers the possibility of choosing the method for optimization. Here we are going exploit this option by providing a widget to select your method.<jupyter_code>chosen_method_optim = ipw.RadioButtons(
    options=['lm', 'trf', 'dogbox'],
    value='lm', # Defaults to 'lm'
    description='Method for optimization',
    style={'description_width': 'initial'},
    disabled=False
)

chosen_method_optim
# perform fit
params, params_covariance = curve_fit(lorentzian, xx, lorentzian_noisy,
                                      method=chosen_method_optim.value,
                                      p0=initial_parameters_values)<jupyter_output><empty_output><jupyter_text>(Top)Displaying the results<jupyter_code># Calculation of the errors on the refined parameters:
params_error = np.sqrt(np.diag(params_covariance))

print('Values of refined parameters:')
print('scale:', params[0],'+/-', params_error[0])
print('center :', params[1],'+/-', params_error[1])
print('HWHM', params[2],'+/-', params_error[2])
# Comparison of reference data with fitting result
fig4 = plt.figure()
gs = fig4.add_gridspec(3, 1)
f4_ax1 = fig4.add_subplot(gs[0:2, :])
f4_ax2 = fig4.add_subplot(gs[2, :])

f4_ax1.plot(xx, lorentzian_noisy, label='reference data')
f4_ax1.plot(xx, lorentzian(xx, *params), '--', label='fit: %5.3f, %5.3f, %5.3f' % tuple(params))
f4_ax2.plot(xx, lorentzian_noisy - lorentzian(xx, *params), label='residuals')
f4_ax1.set_xlabel('x')
f4_ax1.grid()
f4_ax1.legend()
f4_ax2.set_xlabel('x')
f4_ax2.grid()
f4_ax2.legend()<jupyter_output><empty_output><jupyter_text>(Top)ExerciseThis exercise uses the same procedure as above but with slightly different data.  The objective is to determine the new fitting results.  **(1)** Determine the initial guesses for the amplitude, center and width, *Scale*, *Center* and *FWHM* to fit the following data based on a Lorentzian profile with added noise by running the cell below and using the generated widget.<jupyter_code># Define the new data
new_data = lorentzian(xx, 3, 4, 0.5) * (
    1. + 0.1*np.random.normal(0,1,500)) + 0.01*np.random.normal(0,1,500)

initial_params = {"scale": 5.5,
                  "center": 1.0,
                  "hwhm": 2.5}

fig5, ax5 = plt.subplots()
ax5.plot(xx, new_data, label="reference data for exercise")
lines = ax5.plot(xx, lorentzian(xx, *list(initial_params.values())))
ax5.set_ylabel('lorentzian(x,{},{},{})'.format(*list(initial_params.values())))
ax5.set_xlabel('x')
ax5.grid()
ax5.legend()

# define slider to interactively modify the parameters
sliders = {key: ipw.FloatSlider(min=0.1, max=10.0, value=value, description=key)
           for key, value in initial_params.items()}
sliders["center"].min = -sliders["center"].max

# define function to be plotted
def interactive_plot(change):
    params = [sliders[key].value for key in sliders]
    lines[0].set_ydata(lorentzian(xx, *params))
    ax1.set_ylabel('lorentzian(x,{},{},{})'.format(*params))
#     fig1.canvas.draw_idle()
# add observers to the sliders to update the plot
for sl in sliders.values():
    sl.observe(interactive_plot, names="value")

# Define function to reset all parameters' values to the initial ones
def reset_values(b):
    """Reset the interactive plots to inital values."""
    for key, value in initial_params.items():
        sliders[key].value = value

# Define reset button and occurring action when clicking on it
reset_button = ipw.Button(description = "Reset")
reset_button.on_click(reset_values)

# display the interactive plot
display(ipw.VBox(list(sliders.values())), reset_button)<jupyter_output><empty_output><jupyter_text>**(2)** Running the following cells will:- assign the values you determined just above to your model- plot the initial model and the experimental data- run the fit<jupyter_code>initial_parameters_new_data = [sl.value for sl in sliders.values()]
initial_parameters_new_data<jupyter_output><empty_output><jupyter_text>The syntax above uses list comprehensions (as described in another tutorial). Options for such inputs are:- list comprehension```pythonx = [i for i in range(10)]```- with a standard `for` loop:```pythonx = []for i in range(10):    x[i] = i```- entered by hand```python x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]```<jupyter_code>fig6, ax6 = plt.subplots()
ax6.plot(xx, new_data, label="reference data for exercise")
ax6.plot(xx, lorentzian(xx, *initial_parameters_new_data), '.', label='model with initial guesses')
ax6.set_ylabel('lorentzian(x,{},{},{})'.format(*list(initial_params.values())))
ax6.set_title("Exercise: plot before fit")
ax6.set_xlabel('x')
ax6.grid()
ax6.legend()<jupyter_output><empty_output><jupyter_text>Now run the fit<jupyter_code>params_new, pcov_new = curve_fit(lorentzian, xx, new_data, 
                                 p0=initial_parameters_new_data)<jupyter_output><empty_output><jupyter_text>**(3)** Display the final values of the refined parameters<jupyter_code>print('Values of refined parameters:')
print('scale:', params_new[0])
print('center :', params_new[1])
print('HWHM', params_new[2])
print('to be compared to 3, 4, 0.5')<jupyter_output><empty_output><jupyter_text>**(4)** Plot the reference data and the model  In the code below, replace:     - `{A}` with name of reference data   - `{B}` with name of list of values of fitted parameters Refer to the first example for help<jupyter_code>fig7 = plt.figure()
gs = fig7.add_gridspec(3, 1)
f7_ax1 = fig7.add_subplot(gs[0:2, :])
f7_ax2 = fig7.add_subplot(gs[2, :])

f7_ax1.plot(xx, {A}, label='reference data for exercise')
f7_ax1.plot(xx, lorentzian(xx, *{B}), '--', label='fit: %5.3f, %5.3f, %5.3f' % tuple({B}))
f7_ax2.plot(xx, {A} - lorentzian(xx, *{B}), label='residuals')
f7_ax1.set_xlabel('x')
f7_ax1.grid()
f7_ax1.legend()
f7_ax2.set_xlabel('x')
f7_ax2.grid()
f7_ax2.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13016917293233082, 'fraction_numerical': 0.02199248120300752, 'mean_word_length': 5.289007092198582, 'pattern_counts': {'":': 6, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '272293', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3681, 'n_tokens_neox': 3522, 'n_words': 980, 'path': 'notebooks/3_using_external_libraries/4_fitting-scipy.ipynb', 'repo_name': 'celinedurniak/python-course-ikon', 'size': '10640'}