Fixing a common distance between two lines of a doublet for a double-component fit.

[ZIHABAKE Daniel]

Hello lmfit,

i want to constrain the distance between two doublet to be the same for two Gaussian-components.  In my code, i specify the difference between centers  of line 4 and line 5 to be the same with the one between line 3 and 6. But when  i print the value and i calculate , i find the differences are not the same, something that may be caused by error values for centers.  

How can i constrain my parameters D1 and D2 to rely on values of centers of the lines without taking into account errors of centers of lines ?

this is what i think that could help .

there is attached data. 

Thank you.

Daniel

from numpy import loadtxt
import numpy as np
from numpy import linspace
import matplotlib.pyplot as plt
from lmfit import Parameters
import statistics
from lmfit.models import GaussianModel, LinearModel
from matplotlib.axis import Axis
from scipy.interpolate import interp1d
from scipy.optimize import fixed_point

params={'legend.fontsize':'18','axes.labelsize':'18',
'axes.titlesize':'18','xtick.labelsize':'18',
'ytick.labelsize':'18','lines.linewidth':2,'axes.linewidth':2,'animation.html': 'html5'}
plt.rcParams.update(params)
plt.rcParams.update({'figure.max_open_warning': 0})
plt.rcParams['figure.figsize'] = (12,8)


data = loadtxt('/home/zihabake/Desktop/work/Fits/Data1.txt')
x = data[:, 0]
y = data[:, 1]

fun = interp1d(x=x, y=y, kind=25)
x2 = np.linspace(start=5000, stop=5950, num=5000)
y2 = fun(x2)

mask = (x2>= 5875) & (x2 <= 5910)
mask = (x2>= 5840) & (x2 <= 5925)

xarray, yarray = x2[mask], y2[mask]
x=xarray
y=yarray
     

#  Note: use builtin models
gmodel = (LinearModel(prefix='off_')+GaussianModel(prefix='g4_')+GaussianModel(prefix='g5_')+GaussianModel(prefix='g3_')+GaussianModel(prefix='g6_'))

#  Note: many fewer bounds:
params = Parameters()
params.add('off_intercept', value=1)
params.add('off_slope', value=0)

params.add('g5_amplitude', value=-3 )
params.add('g5_center', value=5890, min=5888, max=5892)
params.add('g5_sigma', value=3, min=0)

params.add('alpha1', value=1.5, vary=True, min=1, max=2.00)
params.add('g4_amplitude', expr='g5_amplitude/alpha1')
params.add('g4_center', value=5896, min=5893, max=5900)
params.add('g4_sigma', expr='g5_sigma')



params.add('g3_amplitude', value=-3)
params.add('g3_center', value=5887, min=5885, max=5893)
params.add('g3_sigma', value=3, min=0)
params.add('alpha2', value=1.5, vary=True, min=1, max=2.00)
params.add('g6_amplitude', expr='g3_amplitude/alpha2')
params.add('g6_center', value=5894, min=5890, max=5897)
params.add('g6_sigma', expr='g3_sigma')




params.add('D2', expr='g6_center-g3_center')
params.add('D1', expr='g4_center-g5_center')

params.add('D2', value=6.0000, vary=False)
params.add('D1', expr='D2')


result = gmodel.fit(y, params, x=x)
print(result.fit_report())

comps = result.eval_components(x=x)

fig, axes = plt.subplots(2, 1, gridspec_kw={'hspace': 0, 'height_ratios': [3,1]}, figsize=(12, 8))


axes[0].plot(x, y, 'k', label='observed')
axes[0].set_ylim(ymin=-0.5, ymax=1.25)
axes[0].plot(x, result.best_fit, 'r--', label='best fit')

axes[0].plot(x, comps['g3_']+comps['off_'], 'b-', label='g3')
axes[0].plot(x, comps['g4_']+comps['off_'], 'g-', label='g4')
axes[0].plot(x, comps['g5_']+comps['off_'], 'g-', label='g5')
axes[0].plot(x, comps['g6_']+comps['off_'], 'b-', label='g6')

axes[0].plot(x, comps['off_'], 'k-', label='offset')

residual=result.best_fit-y
axes[1].plot(x,residual)
axes[1].set_ylim(ymin=-0.1, ymax=0.1)
axes[1].set_xlabel(r'wavelength [$\rm \AA$]')
axes[0].set_ylabel(r'flux [$\rm erg / (cm^{2} \, s \, \AA)$]')
axes[1].axhline (y = 0, color='k')

plt.show()

[Matt Newville]

Hi Daniel, 

On Fri, Aug 19, 2022 at 7:05 AM ZIHABAKE Daniel <dziha...@gmail.com> wrote:

Hello lmfit,

i want to constrain the distance between two doublet to be the same for two Gaussian-components.  In my code, i specify the difference between centers  of line 4 and line 5 to be the same with the one between line 3 and 6. But when  i print the value and i calculate , i find the differences are not the same, something that may be caused by error values for centers.  

How can i constrain my parameters D1 and D2 to rely on values of centers of the lines without taking into account errors of centers of lines ?

this is what i think that could help .

Each variable parameter needs to be used in the fit.  You have something like:

    params.add('g5_center', value=5890, min=5888, max=5892)

    params.add('g4_center', value=5896, min=5893, max=5900)

    params.add('g3_center', value=5887, min=5885, max=5893)

    params.add('g6_center', value=5894, min=5890, max=5897)

    params.add('D2', expr='g6_center-g3_center')

    params.add('D1', expr='g4_center-g5_center')
    params.add('D2', value=6.0000, vary=False)

    params.add('D1', expr='D2')

which has 'g3_center',  'g4_center',  'g5_center', and  'g6_center' as variable parameters, and then define 'D2' as the difference between the two variables with  'g6_center-g3_center', and similarly define 'D1' to be  'g4_center-g5_center'.  But then you fix 'D2' to a value of 6.0 and redefine 'D1' to be 'D1'.

I think what you intend is 

    params.add('D1', value=3)

    params.add('g5_center', value=5890, min=5888, max=5892)    

    params.add('g4_center', expr='g5_center + D1')

    params.add('D2', expr='D1')

    params.add('g3_center', value=5887, min=5885, max=5893)

    params.add('g6_center', expr='g3_center + D2')

That is:

    1.  make a *variable* 'D1' for the difference of 'g5_center' and 'g4_center'.

    2.  make 'g5_center' a variable.

    3.  constrain 'g4_center' to be 'g5_center + D1'. 

There will still be 2 variables for 'g4_center' and 'g4_center'.

Using the same approach for peaks 'g3_center' and 'g6_center', but now also constrain 'D1' to be equal to 'D2', so that really only 3 variables ('g5_centerl', 'g3_center', and 'D1') control the centers of the 4 peaks.  

-Matt

[ZIHABAKE Daniel]

Thank you very much Matt, it is giving what I wanted.

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbr_7%2B491gb63pNx9FNa5%2Bs8i29-iujg4p66N9qtCdoQfw%40mail.gmail.com.