2

I have the following situation. There are M independent random walkers on the discrete domain 0, 1, ..., L. We do this for N identical domains. This results in a matrix X where X[i, j] is the position of walker i on domain j. To make a random step, I add an identically shaped matrix with random +1 and -1's to matrix X. Then I deal with the edges. This works well.

However, I want to extend this model to have solid particles that can't pass through each other. This is shown in 2 cases.

  1. One particle is at position i, the second is at position i+1. The first particle moves to the right, while the second moves to the left.
  2. One particle is at position i, the second is at position i+2. The first particle moves to the right, while the second particle moves to the left.

If I do all steps independently, I can check each step manually to see if it's a legal step. However, this is bad O(M^2N) performance. Is there a more efficient way to detect which matrix element pairs X[i,j], X[k, j] result in two particles passing through each other, preferably in a vectorized way? In this way, I can make the simulation skip these steps.

R van Genderen
  • 127
  • 1
  • 6
  • What is the expected behaviour for case 1 and 2? Do you want to sample till there is no collision, or do you want to resolve it manually? ie the first walker goes and then the second ones moves in the other direction? – scleronomic Sep 25 '20 at 16:41
  • I feel like the easiest way to do this is just to not update these two walkers for that time step. – R van Genderen Oct 01 '20 at 11:00

1 Answers1

0

I guess you have to use some form of loops to achieve this, but maybe this helps:

import numpy as np
import matplotlib.pyplot as plt

L = 50
N = 30
W = 20
n_steps = 1000

# Initialize all walkers on the left side
wn0 = np.arange(W)[:, np.newaxis].repeat(N, axis=1)

# Set up the plot
fig, ax = plt.subplots()
worlds = np.zeros((N, L))
worlds[np.arange(N)[np.newaxis, :], wn0] = np.arange(W)[:, np.newaxis]
h = ax.imshow(worlds, cmap='gray_r')  # cmap='tab20')
ax.set_xlabel('Distance in 1D World')
ax.set_ylabel('Ensemble of Worlds')

for _ in range(n_steps):

    r = np.where(np.random.random(wn0.shape) < 0.5, 1, -1)

    wn1 = wn0 + r
    wn1 = np.clip(wn1, 0, L-1)

    # Case 1
    rest_mat = np.zeros_like(wn0, dtype=bool)
    for i in range(W):
        for j in range(i+1, W):
            rest_mat[[[i], [j]], np.logical_and(wn0[i] == wn1[j], wn1[i] == wn0[j])] = True
    wn1[rest_mat] = wn0[rest_mat]

    # Case 2, go from 0->W and than from W->0 to make sure all duplicates are gone
    for i in np.hstack((range(W), range(W)[-2::-1])):
        temp = (wn1[i] == wn1).sum(axis=0) > 1
        wn1[i, temp] = wn0[i, temp]

    # for wn1_j in wn1.T:  Check if there are collisions
    #     assert len(np.unique(wn1_j)) == W

    wn0 = wn1

    worlds = np.zeros((N, L))
    worlds[np.arange(N)[np.newaxis, :], wn1] = np.arange(W)[:, np.newaxis]
    h.set_data(worlds)
    plt.pause(0.1)

Random Walk GIF

scleronomic
  • 4,392
  • 1
  • 13
  • 43