如何编写pythonic代码来修改数组以使其在对其行进行操作的操作下关闭答案

【问题标题】：How to write a pythonic code to modify an array in order to make it closed under an operation acting on its rows如何编写pythonic代码来修改数组以使其在对其行进行操作的操作下关闭
【发布时间】：2020-08-13 21:09:39
【问题描述】：

我有一个数组 A[:,:]，以及一个作用于 A 的行 x=A[:,J] 的操作 Op。因此，对于 A 的每一行 x，我得到 Op(x)。如果 Op(x) 不是 A 的一行，我将它附加到 A。我这样做直到 A 在 Op 下关闭（我假设 Op 不会产生永无止境的循环，即 Op 在一定数量下关闭迭代次数）。在这个过程结束时，给定扩展的 A 在 Op 下闭合，我还想要排列 Pindex 使得 Op(A[:,J])=A[:,Pindex(J)]。

我已经能够编写一个 Python 代码来做到这一点：

import numpy as np

A=np.array([[0,2,3],
            [0,-3,-1],
            [0,4,3]])

def Op(x):
    return [0,-x[2],x[1]-x[2]]



A=A.tolist()

last=len(A)
Pindex=[]

for i,x in enumerate(A):
    found=False 
    xOp=Op(x)
    for j,y in enumerate(A):
        if np.array_equal(y,xOp):
            Pindex.append(j)
            found=True
            break
    if not found:
        A.append(xOp)
        Pindex.append(last)
        last+=1

A=np.asarray(A)      


print A      
print Pindex
print A[Pindex]

但是，在我看来，它并不是很“pythonic”。我想它可以改进，让它更快。有什么建议吗？

附：这是更大代码的一部分，我需要在其中使用数组。我需要将数组转换为列表，因为我需要更新我正在迭代的对象的长度。也许只有使用数组才能做到这一点的更聪明的方法。

附言我不确定问题的标题。如果您有建议，我可以更改它。

【问题讨论】：

在这里使用np.array 有什么意义吗？为什么不使用列表呢？
由于这是一个工作代码，但唯一的问题是性能，我建议将其发布在codereview.stackexchange.com 上。如果您这样做，请删除该问题的副本。
@hpaulj 那是因为这是更大代码的一部分，我需要使用数组。我需要将数组转换为列表，因为我需要更新我正在迭代的对象的长度。也许只有使用数组才能做到这一点的更聪明的方法
@Georgy。到目前为止，我无法发布新问题。我去做。一旦我这样做，我将删除这个问题，如果它没有收到任何答案。
@Gregory 我认为这个问题值得保留。它与本网站上的所有其他“我如何矢量化这个”问题并没有根本不同。

标签： python numpy

【解决方案1】：

如果您矢量化op，您当然可以矢量化循环的内容（python 约定是使用snake_case）。例如，给定您的起始数组，

A = np.array([[0, 2, 3],
              [0,-3,-1],
              [0, 4, 3]])

首先定义op，它可以同时作用于整个事物：

def op(x):
    return np.stack((x[:, 0], -x[:, 2], x[:, 1] - x[:, 2]), axis=1)

您可以通过首先屏蔽A 中已经重复的部分来组合数组。使用@Untubu's asvoid approach 结合in1d。

dt = np.dtype((np.void, A.dtype.itemsize * A.shape[-1]))

def asvoid(arr):
    arr = np.ascontiguousarray(arr)
    if np.issubdtype(arr.dtype, np.floating):
        arr += 0.
    return arr.view(dt).squeeze()

要构建数组，你的循环可以是这样的：

while True:
    B = op(A)
    mask = np.in1d(asvoid(B), asvoid(A), invert=True)
    if mask.any():
        A = np.concatenate((A, B[mask]), axis=0)
    else:
        break

完成后，您将拥有两种排列的行：A 和 op(A)。您可以使用类似https://stackoverflow.com/a/42232761/2988730 中的技术计算从一个到另一个的排序索引：

oa = np.argsort(asvoid(A), axis=0)
ob = np.argsort(asvoid(B), axis=0)
ib = np.empty_like(oa)
ib[ob] = np.arange(B.shape[0])
J = oa[ib]

您也可以增量构建J，因为您知道在给定迭代中A 的每个元素已经有一个J，而上一次迭代（B）的增量可能或可能不是。事实上，如果保证关闭，您应该只需要在不断缩小的B 上调用op。唯一的问题是您必须按排序顺序维护A（argsorted 很好），因为您附加了B 的元素以方便搜索。

【讨论】：