查找一组索引，将一个 NumPy ndarray 的行映射到另一个
Finding a set of indices that maps the rows of one NumPy ndarray to another
我有两个结构化的二维numpy数组，原则上相等，意思是
A = numpy.array([[a1,b1,c1],
[a2,b2,c2],
[a3,b3,c3],
[a4,b4,c4]])
B = numpy.array([[a2,b2,c2],
[a4,b4,c4],
[a3,b3,c3],
[a1,b1,c1]])
不是
意义上的
numpy.array_equal(A,B) # False
numpy.array_equiv(A,B) # False
numpy.equal(A,B) # ndarray of True and False
但从某种意义上说，一个数组 (A) 是 原始 而在另一个数组 (B) 中，数据沿一个轴打乱（可能是沿着行或列）。
什么是 sort/shuffle B 匹配或等于 A 或者排序 A 等于 B 的有效方法？相等性检查确实不重要，只要将两个数组打乱以相互匹配即可。 A 因此 B 具有唯一的行。
我尝试了 view 方法来对两个数组进行排序
def sort2d(A):
A_view = np.ascontiguousarray(A).view(np.dtype((np.void,
A.dtype.itemsize * A.shape[1])))
A_view.sort()
return A_view.view(A.dtype).reshape(-1,A.shape[1])
但这显然在这里不起作用。需要对非常大的阵列执行此操作，因此性能和可伸缩性至关重要。
根据您的示例，您似乎同时打乱了所有列，因此存在映射 A→B 的行索引向量。这是一个玩具示例：
A = np.random.permutation(12).reshape(4, 3)
idx = np.random.permutation(4)
B = A[idx]
print(repr(A))
# array([[ 7, 11, 6],
# [ 4, 10, 8],
# [ 9, 2, 0],
# [ 1, 3, 5]])
print(repr(B))
# array([[ 1, 3, 5],
# [ 4, 10, 8],
# [ 7, 11, 6],
# [ 9, 2, 0]])
我们想要恢复一组索引，idx，使得 A[idx] == B。当且仅当 A 和 B 不包含重复行时，这将是一个唯一映射。
一种有效的*方法是找到可以对 A 中的行进行词法排序的索引，然后找到 B[=60= 中每一行的位置] 将属于 A 的排序版本。 A useful trick 是将 A 和 B 视为一维数组，使用 np.void dtype 将每一行视为单个元素：
rowtype = np.dtype((np.void, A.dtype.itemsize * A.size / A.shape[0]))
# A and B must be C-contiguous, might need to force a copy here
a = np.ascontiguousarray(A).view(rowtype).ravel()
b = np.ascontiguousarray(B).view(rowtype).ravel()
a_to_as = np.argsort(a) # indices that sort the rows of A in lexical order
现在我们可以使用 np.searchsorted 对 B 中的每一行都属于 A[= 的排序版本的位置执行二进制搜索60=]:
# using the `sorter=` argument rather than `a[a_to_as]` avoids making a copy of `a`
as_to_b = a.searchsorted(b, sorter=a_to_as)
从A→B的映射可以表示为A→As→B[的复合=60=]
a_to_b = a_to_as.take(as_to_b)
print(np.all(A[a_to_b] == B))
# True
如果A和B不包含重复行，从B→A[=60=的逆映射]也可以使用
获得
b_to_a = np.argsort(a_to_b)
print(np.all(B[b_to_a] == A))
# True
作为单一函数：
def find_row_mapping(A, B):
"""
Given A and B, where B is a copy of A permuted over the first dimension, find
a set of indices idx such that A[idx] == B.
This is a unique mapping if and only if there are no repeated rows in A and B.
Arguments:
A, B: n-dimensional arrays with same shape and dtype
Returns:
idx: vector of indices into the rows of A
"""
if not (A.shape == B.shape):
raise ValueError('A and B must have the same shape')
if not (A.dtype == B.dtype):
raise TypeError('A and B must have the same dtype')
rowtype = np.dtype((np.void, A.dtype.itemsize * A.size / A.shape[0]))
a = np.ascontiguousarray(A).view(rowtype).ravel()
b = np.ascontiguousarray(B).view(rowtype).ravel()
a_to_as = np.argsort(a)
as_to_b = a.searchsorted(b, sorter=a_to_as)
return a_to_as.take(as_to_b)
基准：
In [1]: gen = np.random.RandomState(0)
In [2]: %%timeit A = gen.rand(1000000, 100); B = A.copy(); gen.shuffle(B)
....: find_row_mapping(A, B)
1 loop, best of 3: 2.76 s per loop
*成本最高的步骤是对行进行快速排序，平均为 O(n log n)。我不确定是否可以做得比这更好。
由于可以对其中一个数组进行打乱以匹配另一个数组，因此没有人阻止我们重新排列两个数组。使用 Jaime's Answer，我们可以 vstack 两个数组并找到唯一的行。然后 unique 返回的逆索引本质上是所需的映射（因为数组不包含重复行）。
为了方便，我们先定义一个unique2d函数：
def unique2d(arr,consider_sort=False,return_index=False,return_inverse=False):
"""Get unique values along an axis for 2D arrays.
input:
arr:
2D array
consider_sort:
Does permutation of the values within the axis matter?
Two rows can contain the same values but with
different arrangements. If consider_sort
is True then those rows would be considered equal
return_index:
Similar to numpy unique
return_inverse:
Similar to numpy unique
returns:
2D array of unique rows
If return_index is True also returns indices
If return_inverse is True also returns the inverse array
"""
if consider_sort is True:
a = np.sort(arr,axis=1)
else:
a = arr
b = np.ascontiguousarray(a).view(np.dtype((np.void,
a.dtype.itemsize * a.shape[1])))
if return_inverse is False:
_, idx = np.unique(b, return_index=True)
else:
_, idx, inv = np.unique(b, return_index=True, return_inverse=True)
if return_index == False and return_inverse == False:
return arr[idx]
elif return_index == True and return_inverse == False:
return arr[idx], idx
elif return_index == False and return_inverse == True:
return arr[idx], inv
else:
return arr[idx], idx, inv
我们现在可以定义我们的映射如下
def row_mapper(a,b,consider_sort=False):
"""Given two 2D numpy arrays returns mappers idx_a and idx_b
such that a[idx_a] = b[idx_b] """
assert a.dtype == b.dtype
assert a.shape == b.shape
c = np.concatenate((a,b))
_, inv = unique2d(c, consider_sort=consider_sort, return_inverse=True)
mapper_a = inv[:b.shape[0]]
mapper_b = inv[b.shape[0]:]
return np.argsort(mapper_a), np.argsort(mapper_b)
验证:
n = 100000
A = np.arange(n).reshape(n//4,4)
B = A[::-1,:]
idx_a, idx_b = row_mapper(A,B)
print np.all(A[idx_a]==B[idx_b])
# True
基准：
针对@ali_m 的解决方案
进行基准测试
%timeit find_row_mapping(A,B) # ali_m's solution
%timeit row_mapper(A,B) # current solution
# n = 100
100000 loops, best of 3: 12.2 µs per loop
10000 loops, best of 3: 47.3 µs per loop
# n = 1000
10000 loops, best of 3: 49.1 µs per loop
10000 loops, best of 3: 148 µs per loop
# n = 10000
1000 loops, best of 3: 548 µs per loop
1000 loops, best of 3: 1.6 ms per loop
# n = 100000
100 loops, best of 3: 6.96 ms per loop
100 loops, best of 3: 19.3 ms per loop
# n = 1000000
10 loops, best of 3: 160 ms per loop
1 loops, best of 3: 372 ms per loop
# n = 10000000
1 loops, best of 3: 2.54 s per loop
1 loops, best of 3: 5.92 s per loop
虽然可能还有改进的余地，但当前的解决方案比ali_m的解决方案慢2-3倍，而且可能有点混乱，而且两个数组都需要映射。只是认为这可能是一个替代解决方案。
我有两个结构化的二维numpy数组，原则上相等，意思是
A = numpy.array([[a1,b1,c1],
[a2,b2,c2],
[a3,b3,c3],
[a4,b4,c4]])
B = numpy.array([[a2,b2,c2],
[a4,b4,c4],
[a3,b3,c3],
[a1,b1,c1]])
不是
意义上的numpy.array_equal(A,B) # False
numpy.array_equiv(A,B) # False
numpy.equal(A,B) # ndarray of True and False
但从某种意义上说，一个数组 (A) 是 原始 而在另一个数组 (B) 中，数据沿一个轴打乱（可能是沿着行或列）。
什么是 sort/shuffle B 匹配或等于 A 或者排序 A 等于 B 的有效方法？相等性检查确实不重要，只要将两个数组打乱以相互匹配即可。 A 因此 B 具有唯一的行。
我尝试了 view 方法来对两个数组进行排序
def sort2d(A):
A_view = np.ascontiguousarray(A).view(np.dtype((np.void,
A.dtype.itemsize * A.shape[1])))
A_view.sort()
return A_view.view(A.dtype).reshape(-1,A.shape[1])
但这显然在这里不起作用。需要对非常大的阵列执行此操作，因此性能和可伸缩性至关重要。
根据您的示例，您似乎同时打乱了所有列，因此存在映射 A→B 的行索引向量。这是一个玩具示例：
A = np.random.permutation(12).reshape(4, 3)
idx = np.random.permutation(4)
B = A[idx]
print(repr(A))
# array([[ 7, 11, 6],
# [ 4, 10, 8],
# [ 9, 2, 0],
# [ 1, 3, 5]])
print(repr(B))
# array([[ 1, 3, 5],
# [ 4, 10, 8],
# [ 7, 11, 6],
# [ 9, 2, 0]])
我们想要恢复一组索引，idx，使得 A[idx] == B。当且仅当 A 和 B 不包含重复行时，这将是一个唯一映射。
一种有效的*方法是找到可以对 A 中的行进行词法排序的索引，然后找到 B[=60= 中每一行的位置] 将属于 A 的排序版本。 A useful trick 是将 A 和 B 视为一维数组，使用 np.void dtype 将每一行视为单个元素：
rowtype = np.dtype((np.void, A.dtype.itemsize * A.size / A.shape[0]))
# A and B must be C-contiguous, might need to force a copy here
a = np.ascontiguousarray(A).view(rowtype).ravel()
b = np.ascontiguousarray(B).view(rowtype).ravel()
a_to_as = np.argsort(a) # indices that sort the rows of A in lexical order
现在我们可以使用 np.searchsorted 对 B 中的每一行都属于 A[= 的排序版本的位置执行二进制搜索60=]:
# using the `sorter=` argument rather than `a[a_to_as]` avoids making a copy of `a`
as_to_b = a.searchsorted(b, sorter=a_to_as)
从A→B的映射可以表示为A→As→B[的复合=60=]
a_to_b = a_to_as.take(as_to_b)
print(np.all(A[a_to_b] == B))
# True
如果A和B不包含重复行，从B→A[=60=的逆映射]也可以使用
获得b_to_a = np.argsort(a_to_b)
print(np.all(B[b_to_a] == A))
# True
作为单一函数：
def find_row_mapping(A, B):
"""
Given A and B, where B is a copy of A permuted over the first dimension, find
a set of indices idx such that A[idx] == B.
This is a unique mapping if and only if there are no repeated rows in A and B.
Arguments:
A, B: n-dimensional arrays with same shape and dtype
Returns:
idx: vector of indices into the rows of A
"""
if not (A.shape == B.shape):
raise ValueError('A and B must have the same shape')
if not (A.dtype == B.dtype):
raise TypeError('A and B must have the same dtype')
rowtype = np.dtype((np.void, A.dtype.itemsize * A.size / A.shape[0]))
a = np.ascontiguousarray(A).view(rowtype).ravel()
b = np.ascontiguousarray(B).view(rowtype).ravel()
a_to_as = np.argsort(a)
as_to_b = a.searchsorted(b, sorter=a_to_as)
return a_to_as.take(as_to_b)
基准：
In [1]: gen = np.random.RandomState(0)
In [2]: %%timeit A = gen.rand(1000000, 100); B = A.copy(); gen.shuffle(B)
....: find_row_mapping(A, B)
1 loop, best of 3: 2.76 s per loop
*成本最高的步骤是对行进行快速排序，平均为 O(n log n)。我不确定是否可以做得比这更好。
由于可以对其中一个数组进行打乱以匹配另一个数组，因此没有人阻止我们重新排列两个数组。使用 Jaime's Answer，我们可以 vstack 两个数组并找到唯一的行。然后 unique 返回的逆索引本质上是所需的映射（因为数组不包含重复行）。
为了方便，我们先定义一个unique2d函数：
def unique2d(arr,consider_sort=False,return_index=False,return_inverse=False):
"""Get unique values along an axis for 2D arrays.
input:
arr:
2D array
consider_sort:
Does permutation of the values within the axis matter?
Two rows can contain the same values but with
different arrangements. If consider_sort
is True then those rows would be considered equal
return_index:
Similar to numpy unique
return_inverse:
Similar to numpy unique
returns:
2D array of unique rows
If return_index is True also returns indices
If return_inverse is True also returns the inverse array
"""
if consider_sort is True:
a = np.sort(arr,axis=1)
else:
a = arr
b = np.ascontiguousarray(a).view(np.dtype((np.void,
a.dtype.itemsize * a.shape[1])))
if return_inverse is False:
_, idx = np.unique(b, return_index=True)
else:
_, idx, inv = np.unique(b, return_index=True, return_inverse=True)
if return_index == False and return_inverse == False:
return arr[idx]
elif return_index == True and return_inverse == False:
return arr[idx], idx
elif return_index == False and return_inverse == True:
return arr[idx], inv
else:
return arr[idx], idx, inv
我们现在可以定义我们的映射如下
def row_mapper(a,b,consider_sort=False):
"""Given two 2D numpy arrays returns mappers idx_a and idx_b
such that a[idx_a] = b[idx_b] """
assert a.dtype == b.dtype
assert a.shape == b.shape
c = np.concatenate((a,b))
_, inv = unique2d(c, consider_sort=consider_sort, return_inverse=True)
mapper_a = inv[:b.shape[0]]
mapper_b = inv[b.shape[0]:]
return np.argsort(mapper_a), np.argsort(mapper_b)
验证:
n = 100000
A = np.arange(n).reshape(n//4,4)
B = A[::-1,:]
idx_a, idx_b = row_mapper(A,B)
print np.all(A[idx_a]==B[idx_b])
# True
基准： 针对@ali_m 的解决方案
进行基准测试%timeit find_row_mapping(A,B) # ali_m's solution
%timeit row_mapper(A,B) # current solution
# n = 100
100000 loops, best of 3: 12.2 µs per loop
10000 loops, best of 3: 47.3 µs per loop
# n = 1000
10000 loops, best of 3: 49.1 µs per loop
10000 loops, best of 3: 148 µs per loop
# n = 10000
1000 loops, best of 3: 548 µs per loop
1000 loops, best of 3: 1.6 ms per loop
# n = 100000
100 loops, best of 3: 6.96 ms per loop
100 loops, best of 3: 19.3 ms per loop
# n = 1000000
10 loops, best of 3: 160 ms per loop
1 loops, best of 3: 372 ms per loop
# n = 10000000
1 loops, best of 3: 2.54 s per loop
1 loops, best of 3: 5.92 s per loop
虽然可能还有改进的余地，但当前的解决方案比ali_m的解决方案慢2-3倍，而且可能有点混乱，而且两个数组都需要映射。只是认为这可能是一个替代解决方案。