Я хочу согласовать набор данных с формой (161,14), где строки - это направление энергии, а столбцы - повторения одного и того же спектра с различными экспериментальными условиями.
В наборе данных должно быть 3 разных пика, поэтому я настроил составную модель из трех голосов.Цель состоит в том, чтобы иметь общие параметры, такие, чтобы центр и ширина voigts были одинаковыми.
Я нашел этот связанный вопрос Python и lmfit: как согласовать несколько наборов данных с общими параметрами?
Однако здесь параметры жестко запрограммированы, поэтому я попытался, как показано ниже.
import h5py
import numpy as np
from lmfit import Parameters, minimize, report_fit
from lmfit.models import VoigtModel, LinearModel
from matplotlib import pyplot as plt
import cProfile
mods = None
c = [530., 531.5, 533.]
c_win = 1
sigma = 0.2
gamma = 0.2
gamma_min = 0.1
gamma_max = 1.
def objective(params, x, data):
""" calculate total residual for fits to several data sets held
in a 2-D array, and modeled by Gaussian functions"""
nx, ndata = data.shape
resid = 0.0 * data[:]
# nx = 1
# make residual per data set
for i in range(ndata):
resid[:, i] = data[:, i] - mods[i].eval(params,x=x)
# resid = data - mods[0].eval(params, x=x)
# now flatten this to a 1D array, as minimize() needs
# print(resid.sum())
return resid.flatten()
def make_param(v, params):
for i in range(3):
v[i].set_param_hint('amplitude', value=1e3)
v[i].set_param_hint('center', value=c[i], min=c[i] - c_win, max=c[i] + c_win)
v[i].set_param_hint('sigma', vary=False, value=sigma)
v[i].set_param_hint('gamma', vary=True, expr='', value=gamma, min=gamma_min, max=gamma_max)
params += v[i].make_params()
f = h5py.File("../../analysis.h5", "a")
raw = f["rawdata"]
proc = f["processed"]
spec_group = raw["Co0001_0042O1s_4600"]
specs = spec_group['sweeps'][()]
x = spec_group['x_b'][()]
specs2 = np.zeros((161, 14))
specs2[:, :] = specs[:, 0, :]
l0 = LinearModel(prefix="l0_")
v0 = VoigtModel(prefix="p0_")
v1 = VoigtModel(prefix="p1_")
v2 = VoigtModel(prefix="p2_")
v = [v0, v1, v2]
params = Parameters()
mod0 = l0 + v0 + v1 + v2
params += l0.make_params(intercept=3000, slope=0)
make_param(v, params)
specs2 = specs2[:, ::4]
mods = [mod0]
for i in range(1, specs2.shape[1]):
l0 = LinearModel(prefix="l0_%i" % i)
v0 = VoigtModel(prefix="p0_%i" % i)
v1 = VoigtModel(prefix="p1_%i" % i)
v2 = VoigtModel(prefix="p2_%i" % i)
params += l0.make_params(intercept=3000, slope=0)
v = [v0, v1, v2]
make_param(v, params)
params['p0_%icenter' % i].expr = 'p0_center'
params['p1_%icenter' % i].expr = 'p1_center'
params['p2_%icenter' % i].expr = 'p2_center'
params['p0_%igamma' % i].expr = 'p0_gamma'
params['p1_%igamma' % i].expr = 'p1_gamma'
params['p2_%igamma' % i].expr = 'p2_gamma'
params['p0_%isigma' % i].expr = 'p0_sigma'
params['p1_%isigma' % i].expr = 'p1_sigma'
params['p2_%isigma' % i].expr = 'p2_sigma'
mods += [l0 + v0 + v1 + v2]
cProfile.run('result = minimize(objective, params, args=(x, specs2))')
# result = minimize(objective, params, args=(x, specs2))#,method='ampgo')
report_fit(result)
plt.figure()
plt.plot(x, specs2[:, 0], x, mods[0].eval(result.params, x=x))
plt.plot(x, specs2[:, -1], x, mods[-1].eval(result.params, x=x))
high = np.max(x)
low = np.min(x)
plt.xlim(high, low)
plt.show()
Код выполняется и подходит удовлетворительно, однако это занимает очень много времени.
Итак, я сделал cprofile, и, похоже, большую часть времени это разбор строк.Предназначено ли это или есть способ сократить это время?
Также я заметил, что для этих 4 спектров нужно было провести 14125 оценок.Довольно много, верно?Я делаю фундаментальную ошибку в том, как я определяю параметры, или другая минимизация лучше для этой конкретной проблемы?
Отчет по профилированию и подгонке: https://pastebin.com/pveD6sRe
First lines of the profiling sorted by total time:
ncalls tottime percall cumtime percall filename:lineno(function)
10163844/1288568 21.010 0.000 42.870 0.000 asteval.py:279(run)
226048 9.725 0.000 32.178 0.000 model.py:775(make_funcargs)
18309888 8.781 0.000 13.574 0.000 model.py:769(_strip_prefix)
169536 6.870 0.000 6.870 0.000 lineshapes.py:63(voigt)
1695690 4.731 0.000 54.555 0.000 parameter.py:745(_getval)