How to increase the efficiency of DLA algorithm?

Question

The algorithm:

1. Place seed at the center of the canvas.
2. Release a random walker from the edge.
3. Random walker sticks to the neighboring sites of the seed/previous points.
4. Repeat N(particles) times.

For a 500x500 matrix with 50k iterations, it would easily take one full day to compute. How to reduce the computation time drastically?

Why 50k particles, cause it's a part of the assignment!

Logically program is correct.

I tried profiling, it didn't tell me a whole lot. The while loop is causing the issue? Any suggestions to improve the efficiency of my code?

import time
import numpy as np
import random
import matplotlib.pyplot as plt 
from numba import jit, cuda, prange

#from numba import roc
starttime = time.time()
arr = [h,w] = [500,500] #track of height and width of canvas
particles = 50000
sticking_coeff = 0.5
canvas = np.zeros((h, w)).astype(int) #the blank n * n matrix
canvas[h//2,w//2] = 1 #make center element is 1
#// to get the floor value
stick = []
stick.append([h//2 + 1, w//2])#below 1 
stick.append([[h//2 - 1, w//2]])#above 1 
stick.append([h//2, w//2 + 1])#right of 1 
stick.append([h//2, w//2 - 1])#left of 1 

#@jit(nopython=True, parallel=True)
#@roc.jit(device=True)
#@numba.jit()
@jit(cache=True)
def walk(A, B, canvas):# A: row B: Column
    while True:
        #print('while')
        x = np.random.randint(4)#0: row 1: column
        if x == 0:#forward option
            A += 1
            #print('A+')
        elif x == 1:
            B += 1
            #print('B+')
        elif x == 2:
            A -= 1 #reverse option
            #print('A-')
        else:
            B -= 1
            #print('B-')

        #take care of overflow 
        if A < 0:
            A = 0
            #print('A=0')
        elif B < 0:
            B = 0
            #print('B=0')
        if A >= h:
            A -= 1
            #print('A-1')
        elif B >= w:
            B -= 1
            #print('B-1')

        positi = [A , B]    
        if positi in stick:
            if np.random.rand() < sticking_coeff:
                for site in [[positi[0] + 1, positi[1]],
                            [positi[0] - 1, positi[1]],
                            [positi[0], positi[1] + 1],
                            [positi[0], positi[1] - 1]]:
                    if site not in stick:
                        stick.append(site)
                canvas[positi[0] , positi[1]] = 1   
                break  # <-- Also, this is require, I think
            else:
                continue

    return canvas

for i in prange(particles):

        print('particle ',i+1)
        selec = random.sample(set(['A','B','C','D']),1)
        #pos1 = randrange(0, len(arr))
        pos = np.random.randint(0, len(arr))
        print(selec)

        if selec == ['A']:
            #arr[:,0]#first column
            #arr[pos,0] = 1
            walk(pos,0,canvas)

        elif selec == ['B']:
            #arr[:,-1]#last column
            #arr[pos,-1] = 1
            walk(pos,-1,canvas)

        elif selec == ['C']:   
            #arr[0,:]#first row
            #arr[0,pos] =1
            walk(0,pos,canvas)

        else:    
            #arr[-1,:]#last row
            #arr[-1,pos] = 1
            walk(-1,pos,canvas)
plt.matshow(canvas)
plt.savefig("rand_walk.png",dpi=2000)
print('That took {} seconds'.format(time.time() - starttime))

Answer 1

You might consider to use a tuple of directions, something like this:

directions = ((1,0), (0,1), (-1,0), (0,-1))

In your while loop, you can do:

direction = directions[np.random.randint(4)]
A += direction[0]
B += direction[1]

If you switch from tuple to numpy arrays (direction and your current position), you could also use numpy.add which might be faster. If this really improve performence you have to measure. To do this, you might use a seed for your random generator, to get reproducible code.

The border check might be faster with NumPy too, using abslolut - if its ok, to change the behavior that way (your praticles wont "clue" at the border, but bounce back).

Precalculation

You could go a step further, by precalculating all posible paths for several steps.

# each direction has same propability - making things easier
# you could simply create a list with an entry for each path of a tree diagram,
# as each path has the same probability
def add_one_step(steps: list):
    random_step = [(1,0), (0,1), (-1,0), (0,-1)]
    if steps == []:
        return random_step
    result = []
    for i, steps_entry in enumerate(steps):
        for step in random_step:
            a = steps_entry[0]+step[0]
            b = steps_entry[1]+step[1]
            result.append((a, b))
    return result


def get_multiple_steps(n=5):
    final_directions = []
    while n > 0:
        final_directions = add_one_step(final_directions)
        n -= 1
    return final_directions


# be careful about how many steps to precalculate! The list lengths go with 4**n
precalculated_steps = []
for i in range(12):
    precalculated_steps.append(get_multiple_steps(i))

You could use such precalculated values, to do several steps in one go.

n = 10
direction = precalculated_steps[n][np.random.randint(4**n)]
A += direction[0]
B += direction[1]

Thats the most simple aproach. If you want to go to higher n values, you have to think about how to reduce precalculation time (this is just a simple brute force calculation of all paths), and how to safe each result only once (and how many times it occurs).

The tricky part is your border and the stick. You have to chose a fitting matrix size depending on your current distance to border and stick.

Outlook - binomial distribution and minimal distance matrix

If you understood the above and really want to speed up your code, I would suggest to look at:

• Binomial distribution
  • you can split stepping in two seperate movements (x and y) - than your radom walk over n steps is simply a binominal distribution
  • with this, you could precalc a quadratic matrix for each n = 1 .. 249 with the probability to reach a position
• precalc minimal distance from sticking and border for each position; you have to update this only when a particle finally sticks - and you only need to check for each position, if this stick-position is closer than the previous value; This matrix is used to get the max number of steps you can do in one go.

By using the minimal distance matrix to determine the maximal number of steps allowed to do, and then using the fitting precalculated stepping matrix the solution should be ready in minutes or seconds. If you can be more lax with the border condition, it would speed up the algorithm even more.

Further disscussion and code example can be found in another question.

Errors

After understanding what your code is intended to do, I think your for-loop over all paricles is broken.

First:

pos = np.random.randint(0, len(arr))  # will set pos to 0 or 1
# should be:
pos = np.random.randint(0, arr[0])
# or better
random_pos_x = np.random.randint(0, arr[0])
random_pos_y = np.random.randint(0, arr[1])

Also, you only start particles at the upper and left border?

walk(pos,-1,canvas)  # why -1?
...
walk(-1,pos,canvas)  # why -1?

# should be:
walk(random_pos_x ,arr[1]-1, canvas)
walk(arr[0]-1, random_pos_y, canvas)

The name arr is an example for a bad name, because it hints to "array" (thus it was hard to spot the mistake above). Better would be samething like area. Or use h and w (personally I would prefer to name thus as height and width; thats easier to read in most cases)

Answer 2

There are a couple of obvious (but minor) contributors to inefficiency here:

• You search an increasing number of stick positions for each particle. If I understand the code correctly, once a particle has "stuck" to a position, no other particle will ever "stick" there, and you can remove the original position from the set.

• You could use a set instead of a list. I don't know whether this would be an optimization or a pessimization. My guess is the latter.

• You are using a lot of lists in places you could use tuples. For example,

---

positi = [A , B]

could be just

positi = (A, B)

And then instead of doing list accesses back into the thing you just built, you can just write what you mean:

stick.remove((A, B))
for site in ((A+1, B), (A-1, B), (A, B+1), (A, B-1)):
    if site not in stick:
        stick.append(site)
canvas[A, B] = 1

If you made stick a set, then you would eliminate that condition because sets are uniqued automatically:

stick.remove((A, B))
for site in ((A+1, B), (A-1, B), (A, B+1), (A, B-1)):
    stick.add(site)
canvas[A, B] = 1

---

Finally, your "take care of overflow" part is doing a lot of redundant tests. Consider rewriting it as

    x = np.random.randint(4)
    if x == 0:
        if (A < h-1): A += 1
    elif x == 1:
        if (B < h-1): B += 1
    elif x == 2:
        if (A > 0) : A -= 1
    else:
        if (B > 0): B -= 1

However, fundamentally, if you must use this brute-force algorithm and it must be fast, you should probably switch to a compiled language and/or a language that supports simultaneous multithreading.