lmfit 阶梯函数和步长

问题描述

我想在图像中放入 2D 形状。过去，我在 Python 中使用 lmfit 成功完成了此操作，并将 2D 函数/数据包装到 1D。在那种情况下，2D 模型是一个平滑函数（具有高斯轮廓的环）。现在我正在尝试做同样的事情，但使用“非平滑功能”并且它没有按预期工作。

这就是我想要做的（猜测和安装是一样的）：

我有意改变了猜测的参数，以便轻松查看它是否按预期移动，但没有任何反应。

我注意到，如果我使用 2D 高斯函数而不是瑞士国旗，这是一个平滑的函数，这可以正常工作（请参阅下面的 MWE）：

所以我猜测问题与瑞士国旗功能不流畅有关。我试图通过添加高斯滤波器（模糊）使其平滑，但它仍然不起作用，即使瑞士国旗图变得非常模糊。

一段时间后，我想到可能使用 lmfit（o 无论在后台的人）的步长太小，无法对瑞士国旗产生任何变化。我想尝试将步长增加到 1，但我不知道该怎么做。

这是我的 MWE（对不起，它仍然很长）：

import numpy as np
import myplotlib as mpl # https://github.com/SengerM/myplotlib
import lmfit

def draw_swiss_flag(fig,center,side,**kwargs):
    fig.plot(
        np.array(2*[side] + 2*[side/2] + 2*[-side/2] + 2*[-side] + 2*[-side/2] + 2*[side/2] + 2*[side]) + center[0],np.array([0] + 2*[side/2] + 2*[side] + 2*[side/2] + 2*[-side/2] + 2*[-side] + 2*[-side/2] + [0]) + center[1],**kwargs,)

def swiss_flag(x,y,center: tuple,side: float):
    # x,y numpy arrays.
    if x.shape != y.shape:
        raise ValueError(f'<x> and <y> must have the same shape!')
    flag = np.zeros(x.shape)
    flag[(center[0]-side/2<x)&(x<center[0]+side/2)&(center[1]-side<y)&(y<center[1]+side)] = 1
    flag[(center[1]-side/2<y)&(y<center[1]+side/2)&(center[0]-side<x)&(x<center[0]+side)] = 1
    return flag

def gaussian_2d(x,side):
    return np.exp(-(x-center[0])**2/side**2-(y-center[1])**2/side**2)

def wrapper_for_lmfit(x,x_pixels,y_pixels,function_2D_to_wrap,*params):
    pixel_number = x # This is the pixel number in the data array
    # x_pixels and y_pixels are the number of pixels that the image has. This is needed to make the mapping.
    if (pixel_number > x_pixels*y_pixels - 1).any():
        raise ValueError('pixel_number (x) > x_pixels*y_pixels - 1')
    x = np.array([int(p%x_pixels) for p in pixel_number])
    y = np.array([int(p/x_pixels) for p in pixel_number])
    return function_2D_to_wrap(x,*params)

data = np.genfromtxt('data.txt') # Read data
data -= data.min().min()
data = data/data.max().max()

guessed_center = (data.sum(axis=0).argmax()+11,data.sum(axis=1).argmax()+11) # I am adding 11 in purpose.
guessed_side = 19

model = lmfit.Model(lambda x,xc,yc,side: wrapper_for_lmfit(x,data.shape[1],data.shape[0],swiss_flag,(xc,yc),side))
params = model.make_params()
params['xc'].set(value = guessed_center[0],min = 0,max = data.shape[1])
params['yc'].set(value = guessed_center[1],max = data.shape[0])
params['side'].set(value = guessed_side,min = 0)
fit_results = model.fit(data.ravel(),params,x = [i for i in range(len(data.ravel()))])

mpl.manager.set_plotting_package('matplotlib')
fit_plot = mpl.manager.new(
    title = 'Data vs fit',aspect = 'equal',)
fit_plot.colormap(data)
draw_swiss_flag(fit_plot,guessed_center,guessed_side,label = 'Guessed')
draw_swiss_flag(fit_plot,(fit_results.params['xc'],fit_results.params['yc']),fit_results.params['side'],label = 'Fitted')

swiss_flag_plot = mpl.manager.new(
    title = 'Swiss flag plot',)
xx,yy = np.meshgrid(np.array([i for i in range(data.shape[1])]),np.array([i for i in range(data.shape[0])]))
swiss_flag_plot.colormap(
    z = swiss_flag(xx,yy,center = (fit_results.params['xc'],side = fit_results.params['side']),)

mpl.manager.show()

和this是data.txt的内容。

解决方法

看来你的代码没问题。问题是，正如您已经猜到的，lmfit 使用的算法不能很好地处理非平滑数据。

默认情况下 lmfit 使用 leas squares 方法。让我们将其更改为方法 'differential_evolution'。

params['side'].set(value=guessed_side,min=0,max=len(data))
fit_results = model.fit(data.ravel(),params,x=[i for i in range(len(data.ravel()))],method='differential_evolution'
                        )

请注意，我需要为最大值添加一些有限值，以防止出现“differential_evolution 需要对所有不同参数进行有限界”消息。

切换到进化算法后，拟合现在看起来像这样：

,

lmfit 中的所有拟合算法（以及 scipy.optimize 就此而言），包括“全局优化器”，真正适用于连续变量（双精度）。当试图找到最优参数值时，大多数算法会在值中做一个非常小的步骤（在 ~1.e-7 级别）来确定导数，然后将使用该导数对最优值进行下一次猜测价值。

您看到的问题是您的模型函数使用参数值作为离散值 - 作为使用 int() 的数组的索引。如果对参数值进行了微小的更改，则不会检测到结果的任何变化 - 算法将决定拟合结果不依赖于该值的微小变化。

所谓的“全局求解器”，如微分进化、盆地跳跃、shgo，认为导数方法会导致“假最小值”，因此会“喷洒参数空间”，其中包含大量候选值，然后使用不同的策略来优化这些结果中的最佳结果以找到最佳值。一般来说，这些运行速度要慢得多（OTOH 运行时很便宜！）并且非常适合可能有多个“最小值”并且您真的想找到其中最好的问题，或者对起始值进行合理猜测的问题很辛苦。

对于您的问题，很明显您可以猜测起始值（中心像素必须在图像上，例如，所以可能猜测“中间”），并且从图像中似乎很可能没有可能会发现很多错误的最小值。这意味着可能不需要全局求解器的费用。

另一种方法是让您的形状对象以图像中的任何连续中心为中心，而不仅仅是以整数像素为中心。当然，您必须将其映射到离散图像，但不需要完全打开/关闭。使用像 scipy.special.erf() 和 erfc() 这样的 sigmoidal 函数将允许您仍然有从“开”到“关”的过渡，但宽度小但有限，渗入相邻像素。这足以让拟合找到中心位置的连续（因此，亚像素！）值。在 1-d 中，这可能看起来像 ::

from scipy.special import erf
def smoothed_window(x,edge1,edge2,width):
    return (erf((x-edge1)/width) + erf((edge2-x)/width))/2.0

对于整数 x 值，0.5（即半个像素）的 width 几乎肯定会允许拟合找到 edge1 和 edge2 的子整数值。 （另外：在代码中或在参数级别强制固定 width 参数或强制其为正数）。

我没有尝试将其扩展到您更复杂的“瑞士国旗”功能，但它应该是可能的，并且也适用于拟合中心值。