使用 'scipy.optimize.minimize' 而不是单个参数来最小化参数列表
Minimizing a list of parameters using 'scipy.optimize.minimize' instead of individual parameters
我正在寻找一种方法来将伪光合成轮廓与观察到的光谱轮廓相匹配。一般的方法是在给定波长处存在对轮廓有贡献的理论线。在多个重叠配置文件的情况下，所有配置文件的总和就是结果。每个配置文件都由高斯宽度、洛伦兹因子、波长偏移和振幅描述。重叠的配置文件除了振幅外所有参数都相同。
到目前为止，我为每种情况下的贡献行数 (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 个值的列表传递（如果是多行，可能会更多），但需要传递一个值，它本身就是一个列表（我希望我是有道理的）。话虽这么说，我不知道如何解决它。欢迎任何建议。
我同时解决了。结果是传递给最小化的初始值自动变成了一个 numpy 数组（至少我不认为我自己把它变成了一个），所以我需要做的就是用 pseudovoigt_sum 函数启动 sig_gam_sft_amp = list(sig_gam_sft_amp) 一切正常。如果有人有类似的问题，我会在这里留下答案。
我正在寻找一种方法来将伪光合成轮廓与观察到的光谱轮廓相匹配。一般的方法是在给定波长处存在对轮廓有贡献的理论线。在多个重叠配置文件的情况下，所有配置文件的总和就是结果。每个配置文件都由高斯宽度、洛伦兹因子、波长偏移和振幅描述。重叠的配置文件除了振幅外所有参数都相同。
到目前为止，我为每种情况下的贡献行数 (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 个值的列表传递（如果是多行，可能会更多），但需要传递一个值，它本身就是一个列表（我希望我是有道理的）。话虽这么说，我不知道如何解决它。欢迎任何建议。
我同时解决了。结果是传递给最小化的初始值自动变成了一个 numpy 数组（至少我不认为我自己把它变成了一个），所以我需要做的就是用 pseudovoigt_sum 函数启动 sig_gam_sft_amp = list(sig_gam_sft_amp) 一切正常。如果有人有类似的问题，我会在这里留下答案。