Optimizing flood-fill alike algorithm (percolation theory)

Question

Here my python code :

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)

The idea is to represent a percolatable field with a matrix, where :

• 0 is a full, unfillable cell
• 1 is an empty, fillable cell
• 2 is a filled cell (will become dry => 3)
• 3 is a dry cell (already filled, hence unfillable)

I want to represent the ratio of dry/total fillable cells as a function of p (the probability that a cell is full in the matrix).

However, my code is extremely unefficient (take a lot of time to complete, even with smalls values).

How can I optimize it ?

Answer 1

This code is not very efficient because while you use numpy, most of the calculations are done element wise (double loops over 2D array elements, etc.) which is slow in Python.

There are two possible approaches, to speed this up,

1. You vectorize your code so it uses optimized numpy operators for arrays. For instance, instead of having a function Fillable(G, i, j) that returns a boolean, you define a function Fillable(G) that returns a boolean numpy array for all the i, j indices at a time (without using loops),

   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
   

   This then allows us to remove the double loop from the Fill(G) function, for instance with,

   def Fill(G):
       G[G==2] = 3
       G[(G==1) & Fillable(G)] = 2
   

   Pretty much most of this code could be rewritten in a similar fashion.

2. Another alternative is to use Cython. In this case the code structure could be kept unmodified, and simply adding types would significantly speed things up. For instance, an optimized function Fill for Cython would be,

   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
   

In any case, you should first profile your code, to see which function calls take the most time. It's possible that just optimizing one or two critical function would speed this up by an order of magnitude.