Monte Carlo estimation of pi using distarray

monte-carlo-pi.py

"""
Estimate pi using a Monte Carlo method.

If we imagine a unit circle inscribed within a unit square, the ratio of
the area of the circle to the area of the square is pi/4. So if a point
is chosen at random within the square, it has a pi/4 probability of
being inside the circle too.

So we choose one N points in the square, count how many are in the
circle and divide by N to give an estimation of pi/4. We then multiply
by 4 to get pi.

The covergence of this method is very slow O(n**-0.5).
Usage:
    $ python pi.py 
"""
from distarray.client import Context
from distarray import odin

# Number of points.
N = 1e6

# Get the context of the engines.
context = Context()

# Construct the function we want to run on each engine.
@odin.local
def pi_montecarlo(n):
    """Get an estimation of pi on each engine."""
    import numpy
    x = numpy.random.rand(n)
    y = numpy.random.rand(n)
    r = numpy.hypot(x, y)
    return 4*numpy.sum(r < 1)/float(n)

# Divide our number of points across the engines.
N_on_each_engine = N/len(context.view)
# Run our function on each engine. We get back a list of estimates.
pi_estimates = pi_montecarlo(N_on_each_engine)
# Average our pi estimates print.
print(sum(pi_estimates)/len(pi_estimates))

@cowlicks what would this look like in just plain 'ol numpy? What would it be using just IPython.parallel without any ODIN / distarray stuff? What benefit does using distarray bring?

Also, what about making the local function just return a random array, and doing the Pi computation with that random array at the Client level? This would make it closer to the plain NumPy solution, which would be good to highlight as well.

We're just getting a feel for things here, so this is mostly exploratory. I'd like to get a strong sense of what value distarray brings over serial NumPy and basic IPython.parallel with these demos. Oh, and performance comparisons between serial NumPy and distarray...