From fb9280db5fbc5b9656ce7dc21557b186dde975f9 Mon Sep 17 00:00:00 2001 From: mstekiel <64463131+mstekiel@users.noreply.github.com> Date: Tue, 20 Aug 2024 18:14:07 +0200 Subject: [PATCH 1/2] DOC: #904 --- examples/example_fit_multi_datasets.py | 80 ++++++++++++-------------- 1 file changed, 37 insertions(+), 43 deletions(-) diff --git a/examples/example_fit_multi_datasets.py b/examples/example_fit_multi_datasets.py index beab3dd9..a509da6b 100644 --- a/examples/example_fit_multi_datasets.py +++ b/examples/example_fit_multi_datasets.py @@ -7,76 +7,70 @@ All minimizers require the residual array to be one-dimensional. Therefore, in the ``objective`` function we need to ``flatten`` the array before returning it. -TODO: this could/should be using the Model interface / built-in models! - """ import matplotlib.pyplot as plt import numpy as np from lmfit import Parameters, minimize, report_fit +from lmfit.models import GaussianModel +############################################################################### +# Create N simulated Gaussian data sets +N = 5 +np.random.seed(2021) +x = np.linspace(-1, 2, 151) +data = [] +for _ in np.arange(N): + params = Parameters() + params.add('amplitude', value=0.60 + 9.50*np.random.rand()) + params.add('center', value=-0.20 + 1.20*np.random.rand()) + params.add('sigma', value=0.25 + 0.03*np.random.rand()) + dat = GaussianModel().eval(x=x, params=params) + np.random.normal(size=x.size, scale=0.1) + data.append(dat) +data = np.array(data) -def gauss(x, amp, cen, sigma): - """Gaussian lineshape.""" - return amp * np.exp(-(x-cen)**2 / (2.*sigma**2)) - - -def gauss_dataset(params, i, x): - """Calculate Gaussian lineshape from parameters for data set.""" - amp = params[f'amp_{i+1}'] - cen = params[f'cen_{i+1}'] - sig = params[f'sig_{i+1}'] - return gauss(x, amp, cen, sig) - - -def objective(params, x, data): +############################################################################### +# The objective function will extract and evaluate a Gaussian from the compound model +def objective(params, x, data, model): """Calculate total residual for fits of Gaussians to several data sets.""" ndata, _ = data.shape resid = 0.0*data[:] # make residual per data set for i in range(ndata): - resid[i, :] = data[i, :] - gauss_dataset(params, i, x) + components = model.components[i].eval(params=params, x=x) + resid[i, :] = data[i, :] - components # now flatten this to a 1D array, as minimize() needs return resid.flatten() - ############################################################################### -# Create five simulated Gaussian data sets -np.random.seed(2021) -x = np.linspace(-1, 2, 151) -data = [] -for _ in np.arange(5): - amp = 0.60 + 9.50*np.random.rand() - cen = -0.20 + 1.20*np.random.rand() - sig = 0.25 + 0.03*np.random.rand() - dat = gauss(x, amp, cen, sig) + np.random.normal(size=x.size, scale=0.1) - data.append(dat) -data = np.array(data) +# Create a composite model by adding Gaussians +model_arr = [GaussianModel(prefix=f'n{i+1}_') for i,_ in enumerate(data)] +model = sum(model_arr[1:], start=model_arr[0]) ############################################################################### -# Create five sets of fitting parameters, one per data set -fit_params = Parameters() +# Prepare the fitting parameters and constrain n2_sigma, ..., nN_sigma to be equal to n1_sigma +fit_params = model.make_params() for iy, y in enumerate(data): - fit_params.add(f'amp_{iy+1}', value=0.5, min=0.0, max=200) - fit_params.add(f'cen_{iy+1}', value=0.4, min=-2.0, max=2.0) - fit_params.add(f'sig_{iy+1}', value=0.3, min=0.01, max=3.0) + fit_params.add(f'n{iy+1}_amplitude', value=0.5, min=0.0, max=200) + fit_params.add(f'n{iy+1}_center', value=0.4, min=-2.0, max=2.0) + fit_params.add(f'n{iy+1}_sigma', value=0.3, min=0.01, max=3.0) -############################################################################### -# Constrain the values of sigma to be the same for all peaks by assigning -# sig_2, ..., sig_5 to be equal to sig_1. -for iy in (2, 3, 4, 5): - fit_params[f'sig_{iy}'].expr = 'sig_1' + if iy>0: + fit_params[f'n{iy+1}_sigma'].expr = 'n1_sigma' ############################################################################### # Run the global fit and show the fitting result -out = minimize(objective, fit_params, args=(x, data)) +out = minimize(objective, fit_params, args=(x, data, model)) report_fit(out.params) + ############################################################################### # Plot the data sets and fits plt.figure() -for i in range(5): - y_fit = gauss_dataset(out.params, i, x) - plt.plot(x, data[i, :], 'o', x, y_fit, '-') +for i, y in enumerate(data): + components = model.eval_components(params=out.params, x=x) + plt.plot(x, y, 'o', x, components[f'n{i+1}_'], '-') + +plt.show() \ No newline at end of file From 3dedabdfe241e17232a7f3abb79be18d2988f055 Mon Sep 17 00:00:00 2001 From: mstekiel <64463131+mstekiel@users.noreply.github.com> Date: Wed, 21 Aug 2024 10:50:35 +0200 Subject: [PATCH 2/2] Revert "DOC: #904" This reverts commit fb9280db5fbc5b9656ce7dc21557b186dde975f9. --- examples/example_fit_multi_datasets.py | 80 ++++++++++++++------------ 1 file changed, 43 insertions(+), 37 deletions(-) diff --git a/examples/example_fit_multi_datasets.py b/examples/example_fit_multi_datasets.py index a509da6b..beab3dd9 100644 --- a/examples/example_fit_multi_datasets.py +++ b/examples/example_fit_multi_datasets.py @@ -7,70 +7,76 @@ All minimizers require the residual array to be one-dimensional. Therefore, in the ``objective`` function we need to ``flatten`` the array before returning it. +TODO: this could/should be using the Model interface / built-in models! + """ import matplotlib.pyplot as plt import numpy as np from lmfit import Parameters, minimize, report_fit -from lmfit.models import GaussianModel -############################################################################### -# Create N simulated Gaussian data sets -N = 5 -np.random.seed(2021) -x = np.linspace(-1, 2, 151) -data = [] -for _ in np.arange(N): - params = Parameters() - params.add('amplitude', value=0.60 + 9.50*np.random.rand()) - params.add('center', value=-0.20 + 1.20*np.random.rand()) - params.add('sigma', value=0.25 + 0.03*np.random.rand()) - dat = GaussianModel().eval(x=x, params=params) + np.random.normal(size=x.size, scale=0.1) - data.append(dat) -data = np.array(data) -############################################################################### -# The objective function will extract and evaluate a Gaussian from the compound model -def objective(params, x, data, model): +def gauss(x, amp, cen, sigma): + """Gaussian lineshape.""" + return amp * np.exp(-(x-cen)**2 / (2.*sigma**2)) + + +def gauss_dataset(params, i, x): + """Calculate Gaussian lineshape from parameters for data set.""" + amp = params[f'amp_{i+1}'] + cen = params[f'cen_{i+1}'] + sig = params[f'sig_{i+1}'] + return gauss(x, amp, cen, sig) + + +def objective(params, x, data): """Calculate total residual for fits of Gaussians to several data sets.""" ndata, _ = data.shape resid = 0.0*data[:] # make residual per data set for i in range(ndata): - components = model.components[i].eval(params=params, x=x) - resid[i, :] = data[i, :] - components + resid[i, :] = data[i, :] - gauss_dataset(params, i, x) # now flatten this to a 1D array, as minimize() needs return resid.flatten() + ############################################################################### -# Create a composite model by adding Gaussians -model_arr = [GaussianModel(prefix=f'n{i+1}_') for i,_ in enumerate(data)] -model = sum(model_arr[1:], start=model_arr[0]) +# Create five simulated Gaussian data sets +np.random.seed(2021) +x = np.linspace(-1, 2, 151) +data = [] +for _ in np.arange(5): + amp = 0.60 + 9.50*np.random.rand() + cen = -0.20 + 1.20*np.random.rand() + sig = 0.25 + 0.03*np.random.rand() + dat = gauss(x, amp, cen, sig) + np.random.normal(size=x.size, scale=0.1) + data.append(dat) +data = np.array(data) ############################################################################### -# Prepare the fitting parameters and constrain n2_sigma, ..., nN_sigma to be equal to n1_sigma -fit_params = model.make_params() +# Create five sets of fitting parameters, one per data set +fit_params = Parameters() for iy, y in enumerate(data): - fit_params.add(f'n{iy+1}_amplitude', value=0.5, min=0.0, max=200) - fit_params.add(f'n{iy+1}_center', value=0.4, min=-2.0, max=2.0) - fit_params.add(f'n{iy+1}_sigma', value=0.3, min=0.01, max=3.0) + fit_params.add(f'amp_{iy+1}', value=0.5, min=0.0, max=200) + fit_params.add(f'cen_{iy+1}', value=0.4, min=-2.0, max=2.0) + fit_params.add(f'sig_{iy+1}', value=0.3, min=0.01, max=3.0) - if iy>0: - fit_params[f'n{iy+1}_sigma'].expr = 'n1_sigma' +############################################################################### +# Constrain the values of sigma to be the same for all peaks by assigning +# sig_2, ..., sig_5 to be equal to sig_1. +for iy in (2, 3, 4, 5): + fit_params[f'sig_{iy}'].expr = 'sig_1' ############################################################################### # Run the global fit and show the fitting result -out = minimize(objective, fit_params, args=(x, data, model)) +out = minimize(objective, fit_params, args=(x, data)) report_fit(out.params) - ############################################################################### # Plot the data sets and fits plt.figure() -for i, y in enumerate(data): - components = model.eval_components(params=out.params, x=x) - plt.plot(x, y, 'o', x, components[f'n{i+1}_'], '-') - -plt.show() \ No newline at end of file +for i in range(5): + y_fit = gauss_dataset(out.params, i, x) + plt.plot(x, data[i, :], 'o', x, y_fit, '-')