提高lmfit中多数据集拟合的速度？字符串调用limi

Question

我想用形状（161，14）拟合数据集，其中行是能量方向，cols是同一光谱在不同实验条件下的重复。在

数据集中应该有3个不同的峰值，所以我建立了一个由三个voigts组成的复合模型。目标是要有共享参数，这样voigts的中心和宽度是相同的。在

我找到了这个相关的问题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

^{pr2}$

Answer 1

好吧，分析通常很棘手，但是您没有包括分析的结果或fit报告，因此也很难做出充分的响应。在

14000个功能评估看起来确实很多。但是我不知道Voigt参数的初始值有多真实。定义三个Voigt函数，然后约束所有参数都相同，这似乎有点奇怪。混合创建一个组合Model，然后使用{}，这看起来也很奇怪。在

对于一个简化的（但应该是相关的？）基于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个函数评估和

^{pr2}$

这对我来说意味着fit需要花费一些时间来计算约束表达式（对于Voigt来说，这对于fwhm和height来说有点复杂），但是它并没有支配运行时。我认为你有更多的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)

{fwm>这两个函数的调用次数都不一样。在

我建议尝试一个类似的策略，也许删除height和{}的参数提示，也许使用

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

看看这是否能改善你的跑步时间。在

我怀疑理解fit为什么需要这么多迭代可能会更有用。如果您有多个Voigt峰值，您可能需要检查它们是交换位置还是重叠。。。在