How to do 1D discrete collision detection as efficiently as possible?

Question

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.

Answer 1

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)