python - 提高 lmfit 中多个数据集拟合的速度？字符串调用限制

Question

我想要拟合形状为 (161,14) 的数据集，其中行是能量方向，列是不同实验条件下相同光谱的重复。

数据集中应该有3个不同的峰值，所以我建立了三个voigts的复合模型。目标是具有共享参数，使得 voigt 的中心和宽度相同。

我发现了这个相关问题Python and lmfit: How to fit multiple datasets with shared parameters?

但是这里的参数是硬连线的，所以我尝试如下。

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)

Answer 1

嗯，分析通常很棘手，但您没有包含分析结果或拟合报告，因此也很难做出充分的响应。

14000 次函数评估确实看起来很多。但我不知道 Voigt 参数的初始值有多现实。定义三个 Voigt 函数然后将所有参数限制为相同似乎有点奇怪。混合创建复合模型然后使用 lmfit.minimize 看起来也非常奇怪。

对于基于 lmfit 示例的简化(但应该相关？)案例(数据来自 https://github.com/lmfit/lmfit-py/blob/master/examples/test_peak.dat ):

#!/usrbin/env python
from numpy import loadtxt
import cProfile
from lmfit.models import  VoigtModel

data = loadtxt('examples/test_peak.dat') # from lmfit/examples
x = data[:, 0]
y = data[:, 1]

mod = VoigtModel()
pars = mod.guess(y, x=x)
pars['gamma'].set(value = 2, vary=True, expr=None)

cProfile.run("out= mod.fit(y, pars, x=x)", sort=1)

print(out.fit_report(min_correl=0.25))

我得到了 54 个函数评估和配置文件输出

42228 function calls (37487 primitive calls) in 0.054 seconds
Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        3    0.021    0.007    0.021    0.007 {built-in method numpy.dot}
33225/150    0.006    0.000    0.013    0.000 asteval.py:279(run)
       59    0.004    0.000    0.004    0.000 lineshapes.py:63(voigt)
 1050/150    0.001    0.000    0.012    0.000 asteval.py:581(on_binop)
  300/225    0.001    0.000    0.008    0.000 asteval.py:744(on_call)
     8364    0.001    0.000    0.001    0.000 {built-in method builtins.isinstance}
      156    0.001    0.000    0.001    0.000 {built-in method builtins.compile}
      412    0.001    0.000    0.013    0.000 parameter.py:740(_getval)

这对我来说表明拟合需要花费一些时间来评估约束表达式(对于 Voigt 来说，fwhm 和 height 有点复杂)，但事实并非如此。 t 主导运行时。我认为你有更多的 Voigt 函数和更多的函数评估，所以它可能更重要。

如果我明确简化约束表达式是不正确的

mod = VoigtModel()
mod.param_hints['fwhm']['expr'] = 'sigma'
mod.param_hints['height']['expr'] = 'amplitude'

pars = mod.guess(y, x=x)
pars['gamma'].set(value = 2, vary=True, expr=None)
cProfile.run("out= mod.fit(y, pars, x=x)", sort=1)

然后我看到分析输出

16723 function calls (16361 primitive calls) in 0.045 seconds
Ordered by: internal time
 ncalls  tottime  percall  cumtime  percall filename:lineno(function)
      3    0.022    0.007    0.022    0.007 {built-in method numpy.dot}
     59    0.005    0.000    0.005    0.000 lineshapes.py:63(voigt)
450/150    0.001    0.000    0.002    0.000 asteval.py:279(run)
    412    0.001    0.000    0.004    0.000 parameter.py:740(_getval)

所以它没有对 asteval 进行那么多的调用，但它似乎也没有运行得更快(FWIW，函数 eval 的数量是相同的)。

我可能建议尝试类似的策略，也许删除 height 和 fwhm 的参数提示，也许使用

 mod.param_hints.pop('fwhm')
 mod.param_hints.pop('height')

看看这是否会改善您的运行时间。

我怀疑了解为什么你的拟合需要如此多次迭代可能会更有用。如果您有多个 Voigt 峰，您可能需要检查它们是否交换位置或重叠...