Parallelizing a Python application

In order to recognize the advantages of parallelism we need an algorithm that is easy to parallelize, complex enough to take a few seconds of CPU time, understandable, and also interesting not to scare away the interested learner. Estimating the value of number \(\pi\) is a classical problem to demonstrate parallel programming.

The algorithm we present is a classical demonstration of the power of Monte Carlo methods. This is a category of algorithms using random numbers to approximate exact results. This approach is simple and has a straightforward geometrical interpretation.

We can compute the value of \(\pi\) using a random number generator. We count the points falling inside the blue circle M compared to the green square N. The ratio 4M/N then approximates \(\pi\).

import random
def calc_pi(N):
    """Computes the value of pi using N random samples."""
    ...
    for i in range(N):
        # take a sample
        ...
    return ...

import random

def calc_pi(N):
    M = 0
    for i in range(N):
        # Simulate impact coordinates
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)

        # True if impact happens inside the circle
        if x**2 + y**2 < 1.0:
            M += 1
    return 4 * M / N

%timeit calc_pi(10**6)

676 ms ± 6.39 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Before we parallelize this program, the inner function must be as efficient as we can make it. We show two techniques for doing this: vectorization using numpy, and native code generation using numba.

We first demonstrate a Numpy version of this algorithm:

import numpy as np

def calc_pi_numpy(N):
    # Simulate impact coordinates
    pts = np.random.uniform(-1, 1, (2, N))
    # Count number of impacts inside the circle
    M = np.count_nonzero((pts**2).sum(axis=0) < 1)
    return 4 * M / N

This is a vectorized version of the original algorithm. A problem written in a vectorized form becomes amenable to data parallelization, where each single operation is replicated over a large collection of data. Data parallelism contrasts with task parallelism, where different independent procedures are performed in parallel. An example of task parallelism is the pea-soup recipe in the introduction.

This implementation is much faster than the ‘naive’ implementation above:

%timeit calc_pi_numpy(10**6)

25.2 ms ± 1.54 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

import dask.array as da

def calc_pi_dask(N):
    # Simulate impact coordinates
    pts = da.random.uniform(-1, 1, (2, N))
    # Count number of impacts inside the circle
    M = da.count_nonzero((pts**2).sum(axis=0) < 1)
    return 4 * M / N

%timeit calc_pi_dask(10**6).compute()

4.68 ms ± 135 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Using Numba to accelerate Python code

Numba makes it easier to create accelerated functions. You can activate it with the decorator numba.jit.

import numba

@numba.jit
def sum_range_numba(a):
    """Compute the sum of the numbers in the range [0, a)."""
    x = 0
    for i in range(a):
        x += i
    return x

Let’s time three versions of the same test. First, native Python iterators:

%timeit sum(range(10**7))

190 ms ± 3.26 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Second, with Numpy:

%timeit np.arange(10**7).sum()

17.5 ms ± 138 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Third, with Numba:

%timeit sum_range_numba(10**7)

162 ns ± 0.885 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)

Numba is hundredfold faster in this case! It gets this speedup with “just-in-time” compilation (JIT) — that is, compiling the Python function into machine code just before it is called, as the @numba.jit decorator indicates. Numba does not support every Python and Numpy feature, but functions written with a for-loop with a large number of iterates, like in our sum_range_numba(), are good candidates.

%time sum_range_numba(10**7)
%time sum_range_numba(10.**7)

CPU times: user 58.3 ms, sys: 3.27 ms, total: 61.6 ms
Wall time: 60.9 ms
CPU times: user 5 µs, sys: 0 ns, total: 5 µs
Wall time: 7.87 µs

@numba.jit
def calc_pi_numba(N):
    M = 0
    for i in range(N):
        # Simulate impact coordinates
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)

        # True if impact happens inside the circle
        if x**2 + y**2 < 1.0:
            M += 1
    return 4 * M / N

%timeit calc_pi_numba(10**6)

13.5 ms ± 634 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)