numpy 矢量化并行更新

Question

我想用以下简单的规则实现一个一维元胞自动机：

1. 如果单元格为 1，相邻单元格为 0，则向右移动
如果单元格为 1，
2. 相邻单元格为 1，则不要移动
3. 所有单元格根据其状态同时更新。
4. 我们有封闭的边界条件。这意味着最后一个像元的相邻像元是第一个像元。

例如：

0 1 1 0 1

更新后：

1 1 0 1 0

我的解决方案是

def update(cells):
    neighbors = np.roll(cells,-1)
    dim = len(cells)
    tmp_cells = np.zeros(dim)
    for i,j in  enumerate(cells):
        if j and not neighbors[i]:
            tmp_cells[i], tmp_cells[(i+1)%dim] = 0, 1
        elif j:
            tmp_cells[i] = 1
    return tmp_cells

这工作正常，但解决方案没有利用np.arrays的所有可能性，并简化为一个简单的list算法。

我以为我可以在cells和neighbors之间找到一个简洁的逻辑，但显然我现在必须睡觉。

一些想法？

Answer 1

要在不循环的情况下获取单元格的值，您需要知道其两侧的邻居。你需要左边，因为如果你是0你的新值取决于你的左邻居，而如果你是1你的新值取决于你的右邻居。

您可以详尽地编写所有 3 单元格组合，对吧？换句话说：

000 -> 0
001 -> 0
010 -> 0 # move to the right
011 -> 1 # stay put
100 -> 1 # left neighbor has moved
101 -> 1 # left neighbor has moved
110 -> 0 # move to the right
111 -> 1 # stay put

您可以轻松地将该表转换为布尔函数。我们可以简化它，但让我们一开始就愚蠢：-x & y & z | x & -y & -z | x & -y & z | x & y & z .

就是这样：

left = np.roll(cells, -1)
right = np.roll(cells, 1)
return (np.logical_not(left) & cells & right | # ...)

现在，您当然需要简化布尔方程*，但这应该可以帮助您入门。

_{*或者退后一步，重新考虑规则。如果您是0，则新值始终是从左邻居复制的;如果你是1，它总是从你的右邻居那里复制的。您可以使用布尔运算符的组合来编写它，但使用掩码赋值可能更简单：result[cells] = left[cells]; result[notcells] = right[notcells] .}

Answer 2

Abarnert 建议的滚动可以工作，但会创建两个副本，从性能角度来看并不理想。

为了获得最佳性能，最好使用视图。像这样：

import numpy as np
cells = [0, 1, 1, 0, 1]
cells = np.concatenate([[0],cells, [0]])    #add padding cells
def boundary(cells):
    """enforce boundary conditions"""
    cells[0] = cells[-2]
    cells[-1] = cells[1]
    return cells
cells = boundary(cells)
centre = cells[1:-1]
left   = cells[0:-2]
right  = cells[2:]
#add your logic here...
print (np.logical_not(left) & centre & right)

Answer 3

基于Eelco Hoogendoorn关于视图性能的评论和abarnet在这里对我有用的解决方案的伟大逻辑提示：

def update_parallel(cells):
    cells = boundary(cells)
    center = cells[1:-1]
    left = cells[0:-2]
    right = cells[2:]
    ones = (center == 1)
    zeros = (center == 0)
    result = np.copy(center)
    result[zeros] = left[zeros]
    result[ones] = right[ones]
    num_moves = sum(np.logical_xor(center, result)) / 2   # number of moves per step
    return result, num_moves