2

I'm trying to simulate a 2-D random walk in python with boundaries (the particle/object will can't cross the boundaries and must return back). However my version is not vectorised and is very slow. How do I implement it without(or minimising) the use of loops.

Here is my approach

def bound_walk():

    Origin = [0, 0] #Starting Point

    #All Possible directiom
    directions = ((-1, 1), (0, 1), (1, 1),
                 (-1, 0)        , (1, 0),
                 (-1, -1), (0, -1), (1, -1))


    #Directions allowed when x-coordinate reaches boundary
    refelectionsx = ((-1, 1), (0, 1),
                 (-1, 0),(-1, -1), (0, -1))

    #Directions allowed when y-coordinate reaches boundary
    refelectionsy = ((-1, 0)        , (1, 0),
                 (-1, -1), (0, -1), (1, -1))

    points = [(0, 0)]
    for i in range(20000):
        direction = choice(directions)
        reflection1 = choice(refelectionsx)
        reflection2 = choice(refelectionsy)
        if Origin[0]>50: #Boundary==50
            Origin[0] += reflection1[0]
        elif Origin[0]<-50:
            Origin[0] -= reflection1[0]
        else:
            Origin[0] += direction[0]

        if Origin[1]>50:
            Origin[1] += reflection2[1]
        elif Origin[1] < -50:
            Origin[1] -= reflection2[1]
        else:
            Origin[1] += direction[1]
        points.append(Origin[:])
    return points
TaOne
  • 21
  • 2

2 Answers2

1

Here is one approach that is fast, but not 100% equivalent to your implementation. The difference is that in my implementation, at the boundary the chance of going in one of the directions parallel to the boundary is half that of the directions retreating from the boundary. That is arguably the better model if you think of the directions as being the result of binning continuous directions, because the boundary cuts the relevant bins in half.

If you try it you'll find that it does 10 million steps more or less instantly.

The trick is that we simply "unroll" space, so we can simulate an uncostrained random walk which is cheap and then in the end we fold it back to the bounding rectangle.

# parameters
>>> directions = np.delete(np.indices((3, 3)).reshape(2, -1), 4, axis=1).T - 1
>>> boundaries = np.array([(-50, 50), (-50, 50)])
>>> start = np.array([0, 0])>>> steps = 10**7
>>>
# "simulation" 
>>> size = np.diff(boundaries, axis=1).ravel()
>>> 
>>> trajectory = np.cumsum(directions[np.random.randint(0, 8, (steps,))], axis=0)
>>> trajectory = np.abs((trajectory + start - boundaries[:, 0] + size) % (2 * size) - size) + boundaries[:, 0]
>>> 
# some sanity checks
# boundaries are respected
>>> print(trajectory.min(axis=0))
[-50 -50]
>>> print(trajectory.max(axis=0))
[50 50]
# step size looks ok
>>> print(np.diff(trajectory, axis=0).min(axis=0))
[-1 -1]
>>> print(np.diff(trajectory, axis=0).max(axis=0))
[1 1]
# histograms of time spent at coordinates looks flat
>>> print(np.bincount(trajectory[:, 0] - boundaries[0, 0]))
[ 50276 100134 100395 100969 101218 101388 101708 100688 101460 102667
 103613 103652 103540 103296 102676 102105 102766 102855 101786 101246
 101442 101152 101020 100498 100637 100588 100100  99745 100034  99878
  99120  98076  98193  98126  97715  98317  98343  97693  97391  96854
  96576  96906  96423  96445  96779  96672  96376  95747  95732  95881
  96833  97149  98490  99692  99519  98800  99497 100070 100065  99816
  99838 100470 100466 100887 100461 100033  99405  99425 100537 100227
 100796 101668 101218 101413 101559 101258 101416 101292 100567 100022
 100266 100770 100882 100519 100326 100795 101066 101293 101667 101666
 101040 101221 101019 100868 101681 100778 100121  98500  98174  98308
  49254]
>>> print(np.bincount(trajectory[:, 1] - boundaries[1, 0]))
[ 52316 104725 104235 103801 102936 102269 102604 102557 102514 103063
 102130 101805 101699 102285 102456 102464 102590 104010 103502 103105
 102784 102927 103430 104750 104671 104836 104547 103280 102131 101548
 101173 101806 101345 101959 101525 101061 101260 100774 100126  98806
  99209 100105  99686 100418 101056 101434 101078 101680 103042 103732
 103003 102047 100832 100489 100809 100429 101325 102420 102282 102205
 101341 100644  99827  99482  98931  98588  97911  97981  97053  96794
  96818  97364  97025  97093  97807  98594  98280  98406  98474  98516
  98555  98713  98381  98296  97600  97374  97423  97092  96238  95771
  95547  95325  94710  94115  93332  92219  91309  91780  92399  92345
  45461]
Paul Panzer
  • 51,835
  • 3
  • 54
  • 99
  • @Paul_Panzer could you please take a look at my question https://stackoverflow.com/questions/74009046/how-to-generate-random-placement-and-trajectory-inside-of-polygon-in-vectorized – illuminato Oct 10 '22 at 15:32
-1

What about using a function (which maps any point on the plane to a point in your desired region) on the output of the random walk. eg.

def map_function(x, lower_bound, upper_bound):
   return (x - lower_bound) % (upper_bound - lower_bound)

This would allow the full flexibility on how you intend to have the boundary conditions (eg. "hard reflecting" or more "soft" -> more likely to slow down /reverse when approaching boundary...)

=> This is equivalent to the answer from @Paul Panzer when you choose the function so it basically "folds" the plane on the desired region

allamiro
  • 58
  • 1
  • 10
Richard Hueck
  • 1
  • 1
  • 1
    Please provide additional details in your answer. As it's currently written, it's hard to understand your solution. – Community Aug 30 '21 at 13:03