我正在寻找一种方法来适应伪voigt合成剖面到观测到的光谱剖面。一般的方法是，在给定的波长有理论线，有助于轮廓。在多个重叠轮廓的情况下，结果是所有轮廓的总和。每个剖面由高斯宽度、洛伦兹因子、波长位移和振幅来描述。重叠剖面除振幅外，其他参数均相同。
到目前为止，我为每一种情况的贡献行数(1-12)创建了单独的函数，这些函数没有问题，因为我知道参数的确切数量并逐一指定它们。然而，这导致了大量的代码行和缺乏通用性。我想做一个函数，它可以拟合任意数量的贡献行参数。我能想到的唯一方法就是把所有的参数放在一个列表中，然后把它们的初始值传递给最小化器。我认为这是我能提供的描述问题的最少的代码，尽管它仍然很长。请忽略一些奇怪的变量名(为了一致性，它们与我老板在他的代码中使用的相同)。
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()
我假设问题是我将初始参数作为4个值的列表传递(如果是多行则可能更多)，但需要传递一个值，这是一个列表本身(我希望我是有意义的)。话虽如此，我也不知道该如何修复它。欢迎提出任何建议。
我在此期间解决了它。事实证明，传递给最小化的init值会自动转换为numpy数组(至少我认为我自己没有将其转换为numpy数组)，因此我所需要做的就是使用sig_gam_sft_amp = list(sig_gam_sft_amp)启动pseudovoigt_sum函数，然后一切都工作了。如果有人有类似的问题，我将把答案留在这里。