Random Walk - Parallel processing

Question

I'm currently implementing a monte-carlo method to solve a diffusion equation. The solution can be expressed as a mathematical expectation of phi(W) where phi is a function (varies accordingly to the diffusion equation) and W is a symmetric random walk stopped at the boundary of the domain. To evaluate the function at a point x, I need to start every walk in the expectation from x.

I want to evaluate the function at a large number of points. So this is what I do:

1. I start simultaneously one walk from every point
2. I use the same steps for every point
3. Once each walk has reached the boundary, I stop
4. I repeat these operations N times, in order to approximate the mathematical expectation by the frequency over N simulations.

My code (Python) looks like this :

for k in range(N): #For each "walk"
    step = 0
    while not(every walk has reach the boundary):
        map(update_walk,points) #update the walk starting from each x in points
        incr(step)

The problem is : it is extremely long since N can be large and the number of points also. I am looking for ANY solution that will help me optimize this code.

I have thought of parallel processing (each walk is independent) with IPython but I did not succeed because it is inside a function (it returned an error like

"could not launch function 'f' because it was not found as 'file.f' " but 'f' is define within file.big_f)

Answer 1

This is one of the few cases I might not use parallel computing. I think it'd be faster to just use numpy.

>>> import numpy as np
>>> def walk(x, n=100, box=.5, delta=.2):
...   np.random.seed()
...   w = np.cumsum(x + np.random.uniform(-delta,delta,n))
...   w = np.where(abs(w) > box)[0]
...   return w[0] if len(w) else n
... 
>>> pwalk = np.vectorize(walk)
>>> pwalk(np.zeros(10))
array([10,  25, 4,   5, 100,  6,  28,   6,  25,  23])

I think you should know how to get the expectation value from there.

You could also have passed a tuple to the last argument of np.random.uniform, and then not needed to use np.vectorize.

Of course, if you wanted to use parallel computing, then you could pick a good map function, and call walk from the map instead of using vectorize as I did above.

>>> from pathos.multiprocessing import ProcessingPool as Pool
>>> p = Pool(4)
>>> p.map(walk, [0]*10)
[8, 7, 39, 7, 36, 7, 22, 27, 18, 31]

Using pathos, so the map can be called easily from the interpreter. https://github.com/uqfoundation/pathos