lmfit 步进函数和步长

Question

我想在图像中拟合 2D 形状。过去，我使用 Python 中的 lmfit 并将 2D function/data 包装成 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, 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 使用最小二乘法。让我们将其更改为方法 '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，认为导数方法会导致“假最小值”，“喷雾参数space”也会如此许多候选值，然后使用不同的策略来优化这些结果中的最佳值以找到最佳值。一般来说，这些比运行慢得多（OTOH 运行时间很便宜！）并且非常适合可能存在多个“最小值”并且您真的想找到其中最好的问题，或者很难准确地猜测起始值。

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

另一种方法是让您的形状对象以图像中任何连续的中心为中心，而不仅仅是以整数像素为中心。当然，您确实必须将其映射到离散图像，但不需要完全 on/off。使用像 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 值，0.5 的 width（即半个像素）几乎肯定会允许拟合以找到 edge1 和 edge2 的子整数值. （另外：在代码或参数级别强制固定 width 参数或强制其为正）。

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

lmfit 步进函数和步长

lmfit Stepped functions and Step size

python

lmfit