I am looking for a way to fit a pseudovoigt synthetic profile to an observed spectral profile. The general approach is that there are theoretical lines at given wavelengths that contribute to the profile. In the case of multiple overlapping profiles, the sum of all the profiles is the result. Each profile is described by a Gaussian width, Lorentz factor, wavelength shift and amplitude. Overlapping profiles have all the parameters the same apart from the amplitude.
So far, I made separate functions for each case of the number of contributing lines (1-12), and those work without issues as I know the exact number of parameters and specify them one by one. However, this led to a lot of lines of code and a lack of generality. I want to make one function that can fit the parameters for an arbitrary number of contributing lines. The only approach I can think of using is putting all the parameters in one list and passing their initial values to the minimizer.
I think this is the least amount of code I can provide that describes the problem, even though it's still quite long. Please disregard some weird variable names (they are the same as the one my boss uses in his code, for consistency).
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
def pseudovoigt_point(fv, s, l, x0, xx):
# Calculate the profile y value for a given wavelength
# xx = wavelength
# x0 = central wavelength
# s = central intensity
# fv = full width at half maximum
# l = Lorentz factor between 0 and 1
def gauss(x_fun, sig_fun):
res = (pow(np.e, -1 * (pow(x_fun, 2) / (2 * pow(sig_fun, 2)))) / sig_fun) / np.sqrt(2 * np.pi)
return res
def lorenz(x_fun, gam_fun):
res = (gam_fun / np.pi) / (x_fun * x_fun + gam_fun * gam_fun)
return res
x = xx - x0
fl = fv * l
if l > 0.9999:
fg = 0
else:
fg = np.sqrt(pow(fv - 0.5346 * fl, 2) - pow(0.2166 * fl, 2))
f_base = pow(fg, 5) + 2.69269 * pow(fg, 4) * fl + \
2.42843 * pow(fg, 3) * pow(fl, 2) + \
4.47163 * pow(fg, 2) * pow(fl, 3) + \
0.07842 * fg * pow(fl, 4) + pow(fl, 5)
f = pow(f_base, 0.2)
gam = f / 2
sig = (f / np.sqrt(2 * np.log(2))) / 2
ffl = fl / f
eta = 1.36603 * ffl - 0.47719 * pow(ffl, 2) + 0.11116 * pow(ffl, 3)
if eta < 0.0001:
vp = gauss(x, sig)
vp0 = gauss(0, sig)
elif eta > 0.9999:
vp = lorenz(x, gam)
vp0 = lorenz(0, gam)
else:
vp = eta * lorenz(x, gam) + (1 - eta) * gauss(x, sig)
vp0 = eta * lorenz(0, gam) + (1 - eta) * gauss(0, sig)
y_new = s * vp / vp0
return y_new
def pseudovoigt_1_profile(x_in_range, sigma, amplitude, gamma_L, x0):
# Calculate a profile with given parameters in a given x range
full_y = []
for point in x_in_range:
y_res = pseudovoigt_point(sigma, amplitude, gamma_L, x0, point)
full_y.append(y_res)
return full_y
def pseudovoigt_sum(x_in_range, sig_gam_sft_amp, x_0):
# Sum up all the contributing profiles
# "sig_gam_sft_amp" is a list with float values of the Gaussian width, Lorentz factor, shift and amplitudes of the contributing lines
sigma = sig_gam_sft_amp[0]
sig_gam_sft_amp.pop(0)
gam = sig_gam_sft_amp[0]
sig_gam_sft_amp.pop(0)
sft = sig_gam_sft_amp[0]
sig_gam_sft_amp.pop(0)
amp = []
for i in range(0, len(sig_gam_sft_amp)):
amp.append(sig_gam_sft_amp[i])
# Each x0 is the resulting wavelength after the shift
x0 = []
for i in range(0, len(x_0)):
x0.append(x_0[i] - sft)
y_list = []
for i in range(0, len(amp)):
y_list.append(pseudovoigt_1_profile(x_in_range, sigma, amp[i], gam, x0[i]))
zipped_list = zip(*y_list)
y_res = [sum(item) for item in zipped_list]
return y_res
def fit_minimize(init_sigma, init_gamma, init_shift, init_amps, x_in_range, y_in_range, the_lines, fit_method):
init_values = [init_sigma, init_gamma, init_shift]
for i in range(0, len(init_amps)):
init_values.append(init_amps[i])
def fit_min(init_vals, xa, ya, lin):
err = []
y_res = pseudovoigt_sum(xa, init_vals, lin)
for j in range(0, len(ya)):
err.append(np.abs(ya[j] - y_res[j]))
error = np.sum(err)
return error
# Setting the bounds (not working without them either)
bounds_low = []
bounds_high = []
bounds_low.append(init_values[0] * 0.5)
bounds_high.append(init_values[0] * 2)
bounds_low.append(0)
bounds_high.append(1)
bounds_low.append(init_values[2] - 10)
bounds_high.append(init_values[2] + 10)
for i in range(0, len(init_amps)):
bounds_low.append(init_amps[i] * 0.5)
bounds_high.append(init_amps[i] * 2)
the_bounds = tuple(zip(bounds_low, bounds_high))
# Checking that the bounds are formated correctly
print("init_values: ", init_values)
print("the_bounds: ", the_bounds)
try:
result = minimize(fit_min, init_values, args=(x_in_range, y_in_range, the_lines), method=fit_method, bounds=the_bounds)
fit_results = []
for i in range(0, len(result.x)):
fit_results.append(result.x[i])
print("fit_results: ", fit_results)
if np.min(fit_results) < 0:
raise Exception("Sorry, no numbers below zero")
except Exception:
print("Exception!")
pass
# Initial values of the parameters for one of the profiles.
init_sigma_1 = 3.75
init_gamma_1 = 0.1
init_shift_1 = -1.58
init_amps_1 = [309.17298]
sig_gam_sft_amp_1 = [3.75, 0.1, -1.58, 309.17298]
x_in_range_1 = [3933.25, 3933.5, 3933.75, 3934.0, 3934.25, 3934.5, 3934.75, 3935.0, 3935.25, 3935.5, 3935.75, 3936.0, 3936.25, 3936.5, 3936.75, 3937.0, 3937.25, 3937.5, 3937.75, 3938.0, 3938.25, 3938.5, 3938.75, 3939.0, 3939.25, 3939.5, 3939.75, 3940.0, 3940.25, 3940.5, 3940.75, 3941.0, 3941.25]
y_in_range_1 = [182.8853759765625, 185.4575653076172, 192.58180236816406, 206.53501892089844, 234.7921905517578, 278.3829650878906, 321.2417297363281, 347.1834716796875, 349.72564697265625, 335.0989990234375, 315.3836975097656, 295.64593505859375, 274.190673828125, 246.56423950195312, 210.0041961669922, 165.1270294189453, 121.44742584228516, 84.3924560546875, 56.95466613769531, 38.16130828857422, 26.14503288269043, 18.620834350585938, 13.414580345153809, 10.1522798538208, 8.189888954162598, 6.834545135498047, 5.831571578979492, 4.917010307312012, 3.957012891769409, 3.1664135456085205, 2.325767993927002, 1.5979336500167847, 1.1800310611724854]
the_lines_1 = [3933.67]
fit_method_1 = "SLSQP"
fit_minimize(init_sigma_1, init_gamma_1, init_shift_1, init_amps_1, x_in_range_1, y_in_range_1, the_lines_1, fit_method_1)
y_in_range_2 = pseudovoigt_sum(x_in_range_1, sig_gam_sft_amp_1, the_lines_1)
# Plotting the original profile and the synthetic fit with the initial parameters to make sure the function itself works as intended
plt.plot(x_in_range_1, y_in_range_1, label="original profile")
plt.plot(x_in_range_1, y_in_range_2, label="synthetic profile")
plt.legend()
plt.show()
I assume that the problem is that I'm passing the initial parameters as a list of 4 values (could be more if it were multiple lines), but would need to pass one value which is a list itself (I hope I'm making sense). That being said, I don't have an idea of how to fix it. Any suggestion would be wellcome.
