How to improve performance for evaluating 'Pi`in Python

Question

I wrote the following code in Python to evaluate the value of Pi . It is called the Monte Carlo method . Obviously, increasing the number of samples, the code becomes slower, and I assume that the slowest part of the code is in the sample part. How can I do it faster?

 from __future__ import division import numpy as np a = 1 n = 1000000 s1 = np.random.uniform(0,a,n) s2 = np.random.uniform(0,a,n) ii=0 jj=0 for item in range(n): if ((s1[item])**2 + (s2[item])**2) < 1: ii = ii + 1 print float(ii*4/(n)) 

Do you offer other (supposedly faster) codes?

Answer 1

The bottleneck here is your for loop. Python for loops are relatively slow, so if you need to iterate over a million items, you can get more speed by avoiding them at all. In this case, it is quite easy. Instead of this:

 for item in range(n): if ((s1[item])**2 + (s2[item])**2) < 1: ii = ii + 1 

do the following:

 ii = ((s1 ** 2 + s2 ** 2) < 1).sum() 

This works because numpy has built-in support for optimized array operations. The loop happens in c instead of python, so it is much faster. I did a quick test so you can see the difference:

 >>> def estimate_pi_loop(x, y): ... total = 0 ... for i in xrange(len(x)): ... if x[i] ** 2 + y[i] ** 2 < 1: ... total += 1 ... return total * 4.0 / len(x) ... >>> def estimate_pi_numpy(x, y): ... return ((x ** 2 + y ** 2) < 1).sum() ... >>> %timeit estimate_pi_loop(x, y) 1 loops, best of 3: 3.33 s per loop >>> %timeit estimate_pi_numpy(x, y) 100 loops, best of 3: 10.4 ms per loop 

Here are some examples of possible operations, so you have a sense of how this works.

Array squaring:

 >>> a = numpy.arange(5) >>> a ** 2 array([ 0, 1, 4, 9, 16]) 

Adding Arrays:

 >>> a + a array([0, 2, 4, 6, 8]) 

Comparison of arrays:

 >>> a > 2 array([False, False, False, True, True], dtype=bool) 

Summing Boolean Values:

 >>> (a > 2).sum() 2 

As you probably understand, there are faster ways to evaluate Pi, but I admit that I have always admired the simplicity and effectiveness of this method.

Answer 2

You have assigned numpy arrays, so you should use them to your advantage.

 for item in range(n): if ((s1[item])**2 + (s2[item])**2) < 1: ii = ii + 1 

becomes

 s1sqr = s1*s1 s2sqr = s2*s2 s_sum = s1sqr + s2sqr s_sum_bool = s_sum < 1 ii = s_sum_bool.sum() print float(ii*4/(n)) 

Where do you square arrays by summing them, checking if the sum is less than 1, and then summing the logical values ​​(false = 0, true = 1) to get the total number that matches the criteria.

Answer 3

I reviewed senderle answer, but if you don't want to change your code too much:

numba is a library designed for this purpose.

Just define your algorithm as a function and add the @jit decorator:

  from __future__ import division import numpy as np from numba import jit a = 1 n = 1000000 s1 = np.random.uniform(0,a,n) s2 = np.random.uniform(0,a,n) @jit def estimate_pi(s1, s2): ii = 0 for item in range(n): if ((s1[item])**2 + (s2[item])**2) < 1: ii = ii + 1 return float(ii*4/(n)) print estimate_pi(s1, s2) 

On my laptop, it gets about 20 times acceleration for n = 100000000 and 3 times acceleration for n = 1000000 .