Numpy:如何模拟晶格中的传播

Question

以下是我的低效Python代码：

import numpy as np 
import random  
matrix = np.empty( (100,100), dtype=bool )
matrix[:,:] = False
matrix[50,50] = True
def propagate(matrix, i, j):
    for (di,dj) in [ (1,0), (0,1), (-1,0), (0,-1) ]:
        (ni,nj) = (i+di, j+dj)
        if matrix[ni,nj] and flip_coin_is_face():
              matrix[i,j] = True 
def flip_coin_is_face():
    return random.uniform(0,1) < 0.5
for k in xrange(1000):
   for i in xrange(1,99):
      for j in xrange(1,99):
          propagate(matrix, i, j)

其基本上从矩阵的中心传播True状态。由于我用Python编写循环和传播规则，所以这当然非常缓慢。

我的问题是，如何使用Numpy索引使其尽可能快？

Answer 1

我能想到一种方法，但它与原始代码不同。也就是说，您可以在每个步骤数组（k循环）中过滤1，将每个值传播到其neiboughhours，即掷骰子4次1，并评估下一个步骤数组。每个操作都可以用numpy一行来完成（对于矩阵使用where、reshape、+和*），因此不会有内部循环。

不同之处在于，我们并没有考虑在一个步骤中传播的值，同时评估所有的更改。事实上，它会减慢传播速度，我想这是明显的，就完成所有矩阵所需的步骤而言。

如果这种方法可以，我可以想出一些代码。

Answer 2

我不擅长numpy，但要在矩阵中"传播"，可以使用广度优先搜索之类的方法。如果你以前没有使用过，它看起来是这样的：

import Queue
def neighbors(i, j, mat_shape):
    rows = mat_shape[0]
    cols = mat_shape[1]
    offsets = [(-1, 0), (1, 0), (0, 1), (0, -1)]
    neighbors = []
    for off in offsets:
        r = off[0]+i
        c = off[1]+j
        if 0 <= r and r <= rows and 0 <= c and c <= cols:
            neighbors.append((r,c))
    return neighbors
def propagate(matrix, i, j):
    # 'parents' is used in two ways. first, it tells us where we've already been
    #  second, it tells us w
    parents = np.empty(matrix.shape)
    parents[:,:] = None
    # first-in-first-out queue. initially it just has the start point
    Q = Queue.Queue()
    # do the first step manually; start propagation with neighbors
    matrix[i,j] = True
    for n in neighbors(i,j,matrix.shape):
        Q.put(n)
        parents[n[0],n[1]] = (i,j)
    # initialization done. on to the propagation
    while not Q.empty():
        current = Q.get() # get's front element and removes it
        parent = parents[current[0], current[1]]
        matrix[current[0], current[1]] = matrix[parent[0], parent[1]] and flip_coin_is_face()
        # propagate to neighbors, in order
        for next in neighbors(current[0], current[1], matrix.shape):
            # only propagate there if we haven't already
            if parents[next[0], next[1]] is None:
                parents[next[0], next[1]] = current
                Q.put(next)
    return matrix

你可能会更聪明，尽早切断传播（因为一旦传播到False，它就再也不会得到True）。但对于100x100，这应该足够快。