如何对两个数组进行排序以获得平滑度答案

【问题标题】：How to sort two arrays for smoothness如何对两个数组进行排序以获得平滑度
【发布时间】：2018-04-03 19:07:17
【问题描述】：

假设我有两个大小相同的数组A 和B。为了确定性，我们假设它们是二维的。两个数组中的值代表一些数据，这些数据应该平滑地依赖于数组中的位置。但是，A 中的一些值已与B 中的相应值交换。我想扭转这种交换，但我很难找到一个标准来告诉我什么时候应该交换两个值。

一个例子应该会有所帮助。这是一些创建两个这样的数组并随机交换它们的一些元素的python代码

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import random

### create meshgrid ###
x = np.linspace(-10,10,15);
y = np.linspace(-10,10,11);

[X,Y] = np.meshgrid(x,y);

### two sufficiently smooth functions on the meshgrid ###
A = -X**2-Y**2;
B = X**2-Y**2-100;

### plot ###
ax=plt.subplot(2,2,1)
im1=ax.imshow(A,extent=[-10, 10, -10, 10])
ax.set_title('A')
ax2=plt.subplot(2,2,2)
im2=ax2.imshow(B,extent=[-10, 10, -10, 10])
ax2.set_title('B')

### randomly exchange a few of the elements of A and B ###
for i in np.arange(0,15):
    for j in np.arange(0,11):
        randNumb = random.random();
        if randNumb>0.8:
            mem=A[j,i];
            A[j,i] = B[j,i];
            B[j,i] = mem;

### plot for comparison ###
ax=plt.subplot(2,2,3)
im1=ax.imshow(A,extent=[-10, 10, -10, 10])
ax2=plt.subplot(2,2,4)
im2=ax2.imshow(B,extent=[-10, 10, -10, 10])

plt.show()

这将导致以下情节：

上面两个图是原始数组 A 和 B，下面两个是改组后的版本。现在的练习是反转这个过程，即从打乱后的版本创建原始数组 A 和 B。

关于我所说的“平滑”的说明。当然，该算法只有在原始数据实际上足够平滑时才有效，这意味着一个数组中的相邻数据点对于所有点的值都足够接近。一个解决方案可能会假设是这种情况。

还要注意，在上面的示例中，这个练习非常容易通过肉眼完成。我在实现这一点时遇到的问题是找到一个好的标准来告诉我何时切换到单元格。

注意，反向转换当然可以重新标记 A 和 B。

【问题讨论】：

标准应该可能是图像中差异（的绝对值）的总和。只要交换减少了差异的总和，直到您不能再进行交换，您就可以使用一些效率较低的算法交换值对，但我不确定这是否一定会收敛。

标签： python arrays algorithm sorting

【解决方案1】：

一种可靠的方法是将 4 个相邻像素值的平均值与每个像素中实际存在的值进行比较。也就是说，对于每个像素，计算 A 和 B 中 4 个相邻像素的平均值，并将其与 both A 和 @987654324 中该像素的实际值进行比较@。以下条件效果很好，实际上是一种最小二乘法：

if (  (A[i, j] - A_mean)**2 + (B[i, j] - B_mean)**2
    > (A[i, j] - B_mean)**2 + (B[i, j] - A_mean)**2
    ):
    # Do swap

这里，A_mean 和 B_mean 是 4 个相邻像素的平均值。

要考虑的另一件重要事情是扫描所有像素不一定足够：可能会发生在一次扫描之后，校正交换使得上述条件可以识别更多应该交换的像素.因此，我们必须遍历数组，直到找到“稳定状态”。

完整代码

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import random

### create meshgrid ###
x = np.linspace(-10,10,15);
y = np.linspace(-10,10,11);

[X,Y] = np.meshgrid(x,y);

### two sufficiently smooth functions on the meshgrid ###
A = -X**2-Y**2;
B = X**2-Y**2-100;

### plot ###
ax=plt.subplot(3,2,1)
im1=ax.imshow(A,extent=[-10, 10, -10, 10])
ax.set_title('A')
ax2=plt.subplot(3,2,2)
im2=ax2.imshow(B,extent=[-10, 10, -10, 10])
ax2.set_title('B')

### randomly exchange a few of the elements of A and B ###
for i in np.arange(0,15):
    for j in np.arange(0,11):
        randNumb = random.random();
        if randNumb>0.8:
            mem=A[j,i];
            A[j,i] = B[j,i];
            B[j,i] = mem;

### plot for comparison ###
ax=plt.subplot(3,2,3)
im1=ax.imshow(A,extent=[-10, 10, -10, 10])
ax2=plt.subplot(3,2,4)
im2=ax2.imshow(B,extent=[-10, 10, -10, 10])

### swap back ###
N, M = A.shape
swapped = True
while swapped:
    swapped = False
    for i in range(N):
        for j in range(M):
            A_mean = np.mean([A[i - 1    , j - 1    ],
                              A[i - 1    , (j + 1)%M],
                              A[(i + 1)%N, j - 1    ],
                              A[(i + 1)%N, (j + 1)%M],
                              ])
            B_mean = np.mean([B[i - 1    , j - 1    ],
                              B[i - 1    , (j + 1)%M],
                              B[(i + 1)%N, j - 1    ],
                              B[(i + 1)%N, (j + 1)%M],
                              ])
            if (  (A[i, j] - A_mean)**2 + (B[i, j] - B_mean)**2
                > (A[i, j] - B_mean)**2 + (B[i, j] - A_mean)**2
                ):
                # Do swap
                A[i, j], B[i, j] = B[i, j], A[i, j]
                swapped = True

### plot result ###
ax=plt.subplot(3,2,5)
im1=ax.imshow(A,extent=[-10, 10, -10, 10])
ax2=plt.subplot(3,2,6)
im2=ax2.imshow(B,extent=[-10, 10, -10, 10])

plt.show()

请注意，上述代码认为数组是周期性的，因为边界处的像素的相邻像素被选择为相对边界上的像素（您在示例中提供的数组就是这种情况）。当索引变为负数时，这种索引换行会自动发生，但不会在索引大于或等于数组的给定维度时发生，这就是使用模运算符 % 的原因。

作为奖励技巧，请注意我如何在不需要临时 mem 变量的情况下交换 A[i, j] 和 B[i, j]。此外，我的外环在第一个维度（长度为 11 的那个）之上，而我的内环在第二个维度（长度为 15 的那个）之上。这比执行反向循环顺序要快，因为现在每个元素都按连续顺序（它们在内存中实际存在的顺序）被访问。

最后，请注意，此方法不能保证始终有效。可能会发生如此多的附近像素被交换而无法找到“正确”的解决方案。然而，无论您选择什么标准来确定是否应交换两个像素都是如此。

编辑（非周期性数组）

对于非周期性数组，边界像素的相邻像素少于 4 个（边缘像素为 3 个，角像素为 2 个）。大致如下：

A_neighbors = []
B_neighbors = []
if i > 0 and j > 0:
    A_neighbors.append(A[i - 1, j - 1])
    B_neighbors.append(B[i - 1, j - 1])
if i > 0 and j < M - 1:
    A_neighbors.append(A[i - 1, j + 1])
    B_neighbors.append(B[i - 1, j + 1])
if i < N - 1 and j > 0:
    A_neighbors.append(A[i + 1, j - 1])
    B_neighbors.append(B[i + 1, j - 1])
if i < N - 1 and j < M - 1:
    A_neighbors.append(A[i + 1, j + 1])
    B_neighbors.append(B[i + 1, j + 1])
A_mean = np.mean(A_neighbors)
B_mean = np.mean(B_neighbors)

请注意，如果邻居越少，该方法的鲁棒性就会降低。您还可以尝试使用最近的 8 个像素（即包括对角线邻居），而不仅仅是 4 个。

【讨论】：

+1。对于我的特定示例，这似乎非常有效！在我接受这个答案之前，我将使用不同的输入对其进行更多测试，以了解它何时崩溃。
经过一些测试，这似乎工作得很好，而且我确实很欣赏该解决方案的优雅。只有一个问题：周期性边界条件。我的示例具有此功能不是故意的，不幸的是我的真实数据没有。有没有一种优雅的方法可以在越过数组末尾时不计算像素？我可能会在你的解决方案中添加一些可怕的 if 语句，但也许你知道一个不错的技巧？

【解决方案2】：

注意：您可能想在 MathOverflow 上重申您的问题。

首先，遍历所有水平和垂直的邻居。（是否也考虑对角线取决于您。）

然后，计算所有“邻居”值的差。

最后，有两种流行的选择：

A) 对所有差异的abs()（绝对值）求和
B) 对所有差的平方求和

您的优化目标是最小化该总和。

变体 A) 通常更直观，而变体 B) 通常更易于使用优化工具进行跟踪。

（A）的问题在于函数是光滑的，但不可微，而 B）是光滑且可微的，即 B）在数学分析下“表现”更好。）

【讨论】：