优化洪水相似算法（渗流理论）

Question

这里是我的python代码：

from numpy import *
from copy import *

def Grid(s, p):
    return random.binomial(1, p, (s,s))

def InitialSpill(G, i, j):
    G[i, j] = 2

def Fillable(G, i, j):
    if i > 0 and G[i - 1, j] == 2:
            return True
    if j > 0 and G[i, j - 1] == 2:
            return True
    if i < len(G) - 1 and G[i + 1, j] == 2:
            return True
    if j < len(G) - 1 and G[i, j + 1] == 2:
            return True
    return False

def Fill(G):
    F = copy(G)
    for i in range(len(G)):
        for j in range(len(G)):
            if F[i, j] == 2:
                G[i, j] = 3 # 3 denote a "dry" cell
            elif F[i, j] == 1 and Fillable(F, i, j):
                G[i, j] = 2 # 2 denote a "filled" cell

def EndReached(G): # Check if all filled cells are dry and if no cells are fillable
    for i in range(len(G)):
        for j in range(len(G)):
            if (G[i, j] == 1 and Fillable(G, i, j)) or G[i, j] == 2:
                    return False
    return True

def Prop(G): # yield the ratio between dry and total fillable cells
    (dry, unrch) = (0, 0)
    for e in G:
        for r in e:
            if r == 1:
                unrch += 1
            if r == 3:
                dry += 1
    if unrch == 0 and dry < 2:
        return 0
    return dry / (unrch + dry)

def Percolate(s, p, i, j): #Percolate a generated matrix of size n, probability p
    G = Grid(s, p)
    InitialSpill(G, i, j)
    while not EndReached(G):
        Fill(G)
    return Prop(G)

def PercList(s, i, j, n, l):
    list_p = linspace(0, 1, n)
    list_perc = []
    for p in list_p:
        sum = 0
        for i in range(l):
            sum += Percolate(s, p, i, j)
        list_perc += [sum/l]
    return (list_p, list_perc)

其思想是用一个矩阵来表示可渗透域，其中：

• 0是一个满的、无法填充的单元格
• 1是一个空的、可填充的单元格
• 2是填充的单元格（将变干=>；3）
• 3是干电池（已填充，因此无法填充）

我想把干电池/可填充电池的比率表示为p的函数（矩阵中电池充满的概率）。在

但是，我的代码效率极低（即使使用很小的值，也要花很多时间才能完成）。在

如何优化它？在

Answer 1

这段代码的效率不高，因为当您使用numpy时，大多数计算都是按元素进行的（2D数组元素上的双循环，等等），这在Python中很慢。

有两种可能的方法，加快速度

1. 您可以将代码矢量化，以便它对数组使用优化的numpy运算符。例如，不是使用返回布尔值的函数Fillable(G, i, j)，而是定义一个函数Fillable(G)，它一次为所有i，j索引返回一个布尔numpy数组（不使用循环）

   def Fillable(G):
        res = np.zeros(G.shape, dtype='bool')
        # if i > 0 and G[i - 1, j] == 2:
        # find indices where this is true
        imask, jmask = np.where(G == 2)
        imask, jmask = imask.copy(), jmask.copy() # make them writable
        imask -= 1   # shift i indices
        imask, jmask = imask[imask>=0], jmask[jmask>=0]
        res[imask, jmask] = True
        # [..] do the other conditions in a similar way
        return res
   

   这样我们就可以从Fill(G)函数中删除双循环，例如

   ^{pr2}$

   几乎大部分代码都可以用类似的方式重写。

2. 另一种选择是使用Cython。在这种情况下，代码结构可以保持不变，简单地添加类型将大大加快速度。例如，Cython的优化函数Fill将是

   cpdef int Fill(int [:,::1] G):
       cdef int [:,::1] F = G.base.copy()
       cdef int i, j, N
       N = G.base.shape[0]
       for i in range(N):
           for j in range(N):
               if F[i, j] == 2:
                   G[i, j] = 3 
               elif F[i, j] == 1 and Fillable(F, i, j):
                   G[i, j] = 2 
       return 0
   

在任何情况下，您应该首先profile your code，以查看哪些函数调用花费的时间最多。仅仅优化一个或两个临界函数就有可能使这个过程加快一个数量级。