How to increase the efficiency of DLA algorithm?

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))

• Can you explain, when your while loop is supposed to finish. How often do you run the loop? – natter1 Jan 20 at 12:06
• It will run till all 50,000 particles are stuck to the center element, in a random sequence that is. – Caleb jones Jan 20 at 12:32
• You now, when I testet for the first time, the while loop run 312778 times when called the first time, 1417116 times when called the second ... – natter1 Jan 20 at 13:18
• Welcome to Code Review! Please edit your question so that the title describes the purpose of the code, rather than its mechanism. We really need to understand the motivational context to give good reviews. Thanks! – Toby Speight Jan 20 at 13:52
• What is the "DLA algorithm"? And do you need to simulate each particle individually, or could you perhaps use math to approximate the expected distribution? – Quuxplusone Jan 20 at 14:42

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
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:
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)

• You could go a step further, by precalculating a probablity matrix for several steps....COULD YOU ELABORATE WITH A CODE? – Caleb jones Jan 21 at 1:41
• @Caleb I updated my answer – natter1 Jan 21 at 11:39
• codeshare.io @natter please have a look, my function placement isn't right. – Caleb jones Jan 21 at 16:13
• @Caleb what do you mean function placement is not right? What you are missing is an implementation of your border (right now your particles could move far away outside the box) and a check, how far you are away from stick (you should not do more steps in one go than the distance to stick). Esp. the latter can be tricky. – natter1 Jan 21 at 17:03
• @natt1 Sweet, I guess it does what it's supposed to do. But again I can't figure out where I have messed up since its left-biased for every run. (codeshare.io/anV13j) <if you can execute it and see what I'm talking about. and compare it what its supposed to look like DLA – Caleb jones Jan 21 at 18:37

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)):
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.

• 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 – Caleb jones Jan 21 at 4:02
• This part can be made much better, any way I can improve it? – Caleb jones Jan 21 at 4:05