python - 具有容差的两个数据集的最长公共(public)子序列

Question

我有四个 1D np.array s: x1, y1, x2, y2 , 其中 x1和 y2具有相同的长度，也是 x2和 y2具有相同的长度，因为它们是数据集对应的 x 和 y 值。 len(x1)和 len(x2)总是不同的。假设 len(x1) > len(x2)目前。这两个数组总是有共同的值，但以一种特殊的方式:这些值并不完全相同，只是在一个公差范围内(因为数值错误等)。公差 = 0.01 的示例:

x1 = np.array([0, 1.01, 1.09, 1.53, -9.001, 1.2, -52, 1.011])
x2 = np.array([1, 1.1, 1.2, 1.5, -9, 82])

我只想保留共同的值(以宽容的方式)。使用较短的数组作为引用，即 x2在这种情况下。 x2 中的第一个值是 1 , 并且在 x1 中有对应的值，即 1.01 .下一个: 1.2在 x2 中也有对应的值, 1.2 .值 1.5没有对应的值，因为 1.53超出公差，因此将其过滤掉，等等。
完整的结果应该是:

x1 = np.array([1.01, 1.09, -9.001, 1.2])
x2 = np.array([1, 1.1, -9, 1.2])

为了更进一步，基于以这种方式过滤 x 值，我想过滤两个数据集的相同索引的 y 值，换句话说，我想找到两个数据集的最长公共(public)子序列。请注意，这里的排序很重要，因为它与 y 值有联系(如果我们先对 x 进行 argsort 排序并重新索引 x 和 y 并不重要)。
我根据 this 所做的尝试回答:

def longest_common_subseq(x1, x2, y1, y2, tol=0.02):
    # sort them first to keep x and y connected
    idx1 = np.argsort(x1)
    x1, y1 = x1[idx1], y1[idx1]
    idx2 = np.argsort(x2)
    x2, y2 = x2[idx2], y2[idx2]
    
    # here I assumed that len(x2) < len(x1)
    idx = (np.abs(x1[:,None] - x2) <= tol).any(axis=1)
    
    return x1[idx], x2[idx], y1[idx], y2[idx]

在这种情况下，y 值可以是任意的，只有形状必须与 x1 匹配和 x2 .例如:

y1 = np.array([0, 1, 2, 3, 4, 5, 6, 7])
y2 = np.array([-1, 0, 3, 7, 11, -2])

尝试运行上面的函数会引发 IndexError: boolean index did not match indexed array along dimension 0 .
我理解:索引数组的长度是错误的，因为 x1和 x2有不同的长度，到目前为止我做不到。有没有很好的方法来实现这一点？
编辑:
如果多个值在公差范围内，则应选择最接近的值。

Answer 1

一种简单的方法是找到所有元素之间的距离:

dist = np.abs(x1 - x2[:, None])

既然你说通常你不会在任何其他元素的公差范围内有多个元素，你可以这样做

i2, i1 = np.nonzero(dist < tol)

如果您有多个匹配项，您可以先修剪匹配项:

i1 = np.argmin(dist, axis=1)
i2 = np.flatnonzero(dist[np.arange(x2.size), i1] < tol)
i1 = i1[i2]

如果原始数据已排序，则索引也会排序(它们将是对角线)。这意味着您可以通过检查索引之间的间距来检查子序列长度。匹配序列的两个索引都将递增 1。

mask = (np.diff(i1) == 1) & (np.diff(i2) == 1)
# smear the mask to include both endpoints
mask = np.r_[False, mask] | np.r_[mask, False]
# pad the mask to ensure proper indexing and find the changeover points
locs = np.diff(np.r_[False, mask, False])
inds = np.flatnonzero(locs)
lengths = inds[1::2] - inds[::2]

您可以从上面的数量中找到最长运行的指数:

k = np.argmax(lengths)
start = inds[2 * k]
stop = inds[2 * k + 1]
longest_x1 = x1[i1[start:stop]]
longest_y1 = y1[i1[start:stop]]
longest_x2 = x2[i2[start:stop]]
longest_y2 = y2[i2[start:stop]]