Based on my answer to another code review question I wrote my own python script for DLA (Diffusion Limited Aggregation).
The DLA algorithm task:
- Place seed at the center of the canvas.
- Release a random walker from the edge.
- Random walker sticks to the neighboring sites of the seed/previous points.
- Repeat N(particles) times.
Also, it was stated to use a 500x500 grid and simulate 50000 particles. And particles where not allowed to leave the area.
A brute force approach is very timeconsuming especially for increasing area size. Therefore I choose a different approach.
- Precalculate possible solutions when doing n steps in one go using multinomial distribution
- calculate a distance_matrix for all positions in the area, with the max. allowed steps that make sure, we do not leave the area (by more than one step; thats needed, or border would become sticking position) and also do not collide with another particle. This also gives automatically information about all stick positions (value == 0; sticking particles get value == -1).
- the distance matrix only has to be recalculated when a particle sticks (and also only a limited part of the whole area). While this is time consuming, it only happens N times (N=50000 particles)
- when a particle does a step, I first check the max distance and then get the corresponding precalculated possible movements and use np.random to choose from that list
This approach could solve the task in 45 minutes (the brute force approach was stated to take more than a day):
Note that 50000 particles seems to be too high a number for a 500x500 area, because in the end, particles are stopped at the border (but I wanted to reproduce the given task ^^).
For the precalculating, I used the module multiprocessing, for larger areas (200x200 and above) this improved the precalculating speed by a factor 6 with 6 cores on my test system. More than 95% of calculation time is spent in
position_list[pos_y, pos_x] += rv.pmf([n_left, n_right, n_down, n_up])
In the end, the precalculation only took 3% of the total simulation time, so I guess I can not improve efficency here (beside maybe saving precalculated values in a file, which could be loaded again for later simulations).
More than 1/3 of the total computation time is spent with DLASimulator.get_random_step()
I have a feeling, that there is still a possibility for significant improvement.
Feedback to readability, efficency and general remarks are welcome.
import bisect
from collections import defaultdict
from multiprocessing.pool import Pool
from multiprocessing.spawn import freeze_support
from typing import Optional, List
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import multinomial
def main():
np.random.seed(0)
dla_simulator = DLASimulator(area_size=(500, 500), max_steps=200) # max_steps is 167 for 500x500
dla_simulator.simulate(particles=50000)
dla_simulator.plot()
print("done")
class DLASimulator:
def __init__(self, area_size: tuple, max_steps: int, initial_particles_pos: Optional[List[tuple]] = None, max_factor=1.0):
self.area_size = area_size
if initial_particles_pos is None:
initial_particles_pos = [(area_size[0] // 2, area_size[1] // 2)]
self.initial_particles_pos = initial_particles_pos
self.distance_matrix = None # only to help code completion
self.init_distance_matrix()
self.is_first_particle = True
self.set_stick_particles(self.initial_particles_pos)
self.max_steps = min(max_steps, self.distance_matrix.max()) # no need to precalculate more than max distance
self.max_steps = max(int(self.max_steps*max_factor), 1) # no need to precalculate more than max distance
self.list_of_position_probability_dicts = generate_list_of_position_probability_dicts(self.max_steps)
def reset(self):
self.init_distance_matrix()
self.set_stick_particles(self.initial_particles_pos)
def simulate(self, particles: int = 1000, reset: bool = False):
print("start simulation")
if reset:
self.reset()
for _ in range(particles):
particle = self.get_new_particle()
self.simulate_particle(particle)
print(f"particle {_} sticked")
def simulate_particle(self, particle):
while True:
distance = self.distance_matrix[particle[1], particle[0]]
if distance == 0:
self.set_stick_particle(particle)
break
elif distance < 0:
# raise Exception("There is already a particle on This position")
break
else:
pos_x, pos_y = self.get_random_step(distance)
particle[0] = min(max(particle[0] + pos_x, 0), self.area_size[0]-1)
particle[1] = min(max(particle[1] + pos_y, 0), self.area_size[1]-1)
# calc distance to border for all positions (allowed number of steps in one go for each position)
def init_distance_matrix(self):
size_x, size_y = self.area_size
self.distance_matrix = np.zeros(self.area_size, np.int16)
for pos_x in range(size_x):
for pos_y in range(size_y):
self.distance_matrix[pos_y, pos_x] = min(pos_x + 1, pos_y + 1, size_x-pos_x, size_y - pos_y)
def set_stick_particles(self, particles):
for particle in particles:
self.set_stick_particle(particle)
def set_stick_particle(self, particle):
pos_x, pos_y = particle
self.distance_matrix[pos_y, pos_x] = -1 # no other particle allowed on this position
self.update_distance_positions(particle)
def update_distance_positions_in_rect(self, x_min, x_max, y_min, y_max, particle):
particle_x, particle_y = particle
for pos_x in range(x_min, x_max+1):
for pos_y in range(y_min, y_max + 1):
distance_to_stick = max(abs(pos_x - particle_x) + abs(pos_y - particle_y) - 1, 0)
self.distance_matrix[pos_y, pos_x] = min(self.distance_matrix[pos_y, pos_x], distance_to_stick)
def update_distance_positions(self, particle):
particle_x, particle_y = particle
x_min, y_min = 0, 0
x_max, y_max = self.area_size[0]-1, self.area_size[1]-1
# go in 4 directions; as soon as distance-value is smaller than distance from particle -> stop
for i in range(particle_x-1, -1, -1):
if self.distance_matrix[particle_y, i] <= particle_x - 1 - i:
x_min = i+1
break
for i in range(particle_y-1, -1, -1):
if self.distance_matrix[i, particle_x] <= particle_y - 1 - i:
y_min = i+1
break
for i in range(particle_x+1, self.area_size[0]):
if self.distance_matrix[particle_y, i] <= i - (particle_x + 1):
x_max = i-1
break
for i in range(particle_y+1, self.area_size[1]):
if self.distance_matrix[i, particle_x] <= i - (particle_y + 1):
y_max = i-1
break
self.update_distance_positions_in_rect(x_min, x_max, y_min, y_max, particle)
def get_new_particle(self) -> list:
random = np.random.randint(4)
if random == 0:
pos_x = 0
pos_y = np.random.randint(self.area_size[1])
elif random == 1:
pos_x = self.area_size[0] - 1
pos_y = np.random.randint(self.area_size[1])
elif random == 2:
pos_x = np.random.randint(self.area_size[0])
pos_y = 0
elif random == 3:
pos_x = np.random.randint(self.area_size[0])
pos_y = self.area_size[1] - 1
else:
raise Exception("Something went wrong in get_new_particle()")
return [pos_x, pos_y]
def get_random_step(self, distance):
n_step_dict = self.list_of_position_probability_dicts[min(distance, self.max_steps)]
random = np.random.random()
# search on numpy array seems to be slower then bisect on list
# index = np.searchsorted(n_step_dict["cumulative probability"], random)
index = bisect.bisect_left(n_step_dict["cumulative probability"], random)
pos_x = n_step_dict["pos_x"][index]
pos_y = n_step_dict["pos_y"][index]
return pos_x, pos_y
def plot(self):
canvas = np.zeros(self.area_size)
for x in range(self.area_size[0]):
for y in range(self.area_size[1]):
if self.distance_matrix[y, x] == -1:
canvas[y, x] = 1
plt.matshow(canvas)
plt.show()
plt.savefig("rand_walk_500particles.png", dpi=150)
def calc_symmetric_positions(position_list: dict) -> dict:
result = {}
for key, value in position_list.items():
pos_x = key[0]
pos_y = key[1]
result[-pos_x, pos_y] = position_list[(pos_x, pos_y)]
result[-pos_x, -pos_y] = position_list[(pos_x, pos_y)]
result[pos_x, -pos_y] = position_list[(pos_x, pos_y)]
result[pos_y, pos_x] = position_list[(pos_x, pos_y)]
result[-pos_y, pos_x] = position_list[(pos_x, pos_y)]
result[-pos_y, -pos_x] = position_list[(pos_x, pos_y)]
result[pos_y, -pos_x] = position_list[(pos_x, pos_y)]
return result
# no need for complete matrix, only list with actually reachable positions;
# creation is 1.055 times faster (compared to full matrix)
def calc_position_probability_list_symmetric(number_of_steps) -> dict:
result = {"cumulative probability": [], "pos_x": [], "pos_y": []}
# position_list = defaultdict(lambda:0.0) # multiprocessing needs to pickle -> can't have a lambda
position_list = defaultdict(float)
if number_of_steps == 0:
result["cumulative probability"].append(1)
result["pos_x"].append(0)
result["pos_y"].append(0)
return result
p_list = [0.25, 0.25, 0.25, 0.25] # all directions have same probability
rv = multinomial(number_of_steps, p_list)
# this loop finds all reachable positions
for n_right in range(number_of_steps + 1):
for n_left in range(number_of_steps - n_right + 1):
pos_x = n_right - n_left
if pos_x > 0: # due to symmetry
continue
for n_down in range(number_of_steps - n_left - n_right + 1):
n_up = number_of_steps - n_left - n_right - n_down
pos_y = n_up - n_down
# if pos_y > 0: # unnecessary; pos_x is already <=0
# continue
if pos_y > pos_x: # due to symmetry
continue
# this takes up >95% of computation time in calc_position_probability_list_symmetric!
# thus optimizing loop with impact on readability not helpful
# todo: check with pypy
position_list[pos_y, pos_x] += rv.pmf([n_left, n_right, n_down, n_up])
position_list.update(calc_symmetric_positions(position_list))
# create cummulative distribution values:
# todo: consider sorting before creating cumulative, to allow faster algorithm SLASimulator.in get_random_step()
sum=0
for key, value in position_list.items():
sum += value
result["cumulative probability"].append(sum)
result["pos_x"].append(key[0])
result["pos_y"].append(key[1])
# no speedup for bisect (np.array is even slower then list)
# result["cumulative probability"] = np.array(result["cumulative probability"])
# result["cumulative probability"] = tuple(result["cumulative probability"])
assert np.isclose(sum, 1.0)
return result
def generate_list_of_position_probability_dicts_single_process(max_number_of_steps: int) -> List[dict]:
result = []
for i in range(max_number_of_steps+1):
result.append(calc_position_probability_list_symmetric(i))
return result
def generate_list_of_position_probability_dicts(max_number_of_steps: int, multiprocessing=True) -> List[dict]:
if not multiprocessing:
return generate_list_of_position_probability_dicts_single_process(max_number_of_steps)
with Pool() as pool: # processes=None (or no argument given) -> number returned by os.cpu_count() is used.
result = pool.map(calc_position_probability_list_symmetric, range(max_number_of_steps+1))
return result
if __name__ == '__main__':
freeze_support() # needed for windows
main()
(the newest version can be found here)
