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

Question

我有四个一维 np.array:x1、y1、x2、y2，其中 x1 和 y2 code> 具有相同的长度，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])

我只想保留共同的值(value)观(以宽容的方式)。使用较短的数组作为引用，在本例中为 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: bool 索引与维度 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]

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

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]]