lmfit分步函数和步长

Question

我想要适合在图像中的2D形状。在过去，我已经在Python语言中使用lmfit并将2D函数/数据包装到1D中，成功地做到了这一点。在这种情况下，2D模型是一个平滑函数(具有高斯轮廓的环)。现在我正在尝试做同样的事情，但是使用了一个“非平滑函数”，它并没有像预期的那样工作。

这就是我想要做的(猜测和拟合是一样的)：

​

​

我故意改变了猜测的参数，以便轻松地查看它是否按预期移动，但什么也没有发生。

我注意到，如果我使用2D高斯，而不是瑞士国旗，这是一个平滑的函数，它可以很好地工作(参见下面的MWE )：

​

​

所以我猜这个问题与瑞士国旗功能不平滑有关。我试图通过添加高斯滤波器(模糊)来使其平滑，但它仍然不起作用，即使瑞士国旗的地块变得非常模糊。

过了一段时间后，我开始思考，也许使用lmfit的步长(不管是谁在后台)太小了，不能在瑞士国旗上产生任何变化。我想试着将步长增加到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, y, center, 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, y, *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], min = 0, 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',
    aspect = 'equal',
)
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'],fit_results.params['yc']), side = fit_results.params['side']),
)

mpl.manager.show()

而this是data.txt的内容。

Answer 1

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

默认情况下，lmfit使用leas square方法。让我们将其改为方法'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需要所有可变参数的有限界限“的消息。

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

​

​

Answer 2

lmfit中的所有拟合算法( scipy.optimize也是如此)，包括“全局优化器”在内，实际上都是针对连续变量(双精度)的。当试图找到最优参数值时，大多数算法将在该值中进行非常小的一步(在~1.e-7级别)以确定导数，然后将导数用于下一次猜测最优值。

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

所谓的“全局求解器”，如差分进化，盆地跳跃，shgo，认为导数方法可能导致“假最小值”，因此将“喷雾参数空间”与大量的候选值，然后使用不同的策略来改进这些结果中的最佳，以找到最佳值。一般来说，这些程序运行起来要慢得多(OTOH运行时很便宜！)对于可能存在多个“最小值”并且你真的想找到其中最好的一个的问题，或者很难对起始值进行合理的猜测的问题，非常好用。

对于您的问题，很明显您可以猜测起始值(中心像素必须在图像上，所以可以猜测“中间”)，并且从图像中似乎没有太多可能找到的错误最小值。这意味着可能不需要全局求解器的开销。

另一种方法是允许您的形状对象在图像中的任何连续中心居中，而不仅仅是整数像素。当然，您必须将其映射到离散图像，但它不需要完全打开/关闭。使用像scipy.special.erf()和erfc()这样的S型函数将允许你仍然有一个从“开”到“关”的过渡，但宽度很小但有限，会渗入相邻的像素。这将足以让拟合找到一个连续的(因此，子像素！)中心位置的值。在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值，width为0.5 (即半个像素)几乎肯定会允许拟合查找edge1和edge2的子整数值。(旁注:强制width参数固定或强制其为正，无论是在代码中还是在参数级别)。

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