4 Parallel Simulation in Python

The goal of this chapter is much the same as the previous one. We have a parallel program, sim_template.py, for running a simulation experiment on a computer with only a few cores such as a personal machine. The development this time is in Python. In the next section we give a listing for the code and provide a cursory explanation of how to use it. Then, in the subsequent sections, we illustrate its use with a couple of examples and give a more detailed treatment of how the code works.

4.1 sim_template.py

The code listing below, sim_template.py, represents a Python version of the program sim_template.R developed in the previous chapter.

Listing 4.1. Python simulation template.

from sim_worker import worker
from numpy.random import default_rng, SeedSequence
from multiprocessing import Pool

def sim_template(seed, n, n_reps, n_cores, *args):
    # create a list of random number generators
    ss = SeedSequence(seed)
    child_seeds = ss.spawn(n_cores)
    rg_list = [default_rng(s) for s in child_seeds]
    
    # calculate the number of replications for each core  
    base_reps = n_reps // n_cores
    remainder = n_reps % n_cores
    reps = [base_reps] * n_cores

    # load balance
    for i in range(remainder):
        reps[i] += 1

    # information to be passed to workers
    set_up = [[0] * 4 for _ in range(n_cores)]
    for i in range(n_cores):
        set_up[i][0] = rg_list[i]
        set_up[i][1] = n
        set_up[i][2] = reps[i]
        set_up[i][3] = args

    with Pool(processes = n_cores) as pool:
        result = pool.starmap(worker, set_up)

    return result

Much like its R analog, understanding this Python spin on parallel simulation reduces to understanding the function signature:

sim_template(seed, n, n_reps, n_cores, *args)

The arguments are

seed: the user supplied seed for the random number generator used in the program,
n: the size for each sample,
n_reps: the number of replicate samples of size n to be used in the simulation experiment,
n_cores: the number of worker processes to apply to the computations and
*args: the ubiquitous Python catch-all for a variable number of arguments. The purpose of *args here is to pass in parameters that are needed to create data from a statistical model.

Beyond sim_template we need a Python module sim_worker.py that implements your particular estimation methodology. A sketch of what it should look like is given below.

def worker(rg, n, n_reps, param_tuple):
    U = rg.random((n, n_reps))
    YOUR CODE GOES HERE

The function takes four arguments: the random number generator, the sample size, the number of samples, and a tuple deriving from the *args sim_template argument. The first line of code uses the random method for the random number generator to create an n by n_reps array of random uniform deviates. These values will then be used along with possibly a quantile function transformation from Section 2.2 and parameter values in args to create data following some statistical model and compute values for a corresponding statistic. That’s what your code will do. The Python scripts sim_template.py and sim_worker.py should be in the same directory.

4.2 Estimating \(\pi\) revisited

In Chapter 2 we saw that we could estimate \(\pi\) by throwing darts at a square two units on a side centered at the origin with a circle drawn in its interior. We tossed the darts and then counted the fraction falling in the circle. Four times this fraction serves as an estimator of \(\pi\).

This section carries out a virtual version of dart throwing on a computer. Specifically, we randomly generate dart coordinates and check if they fall in the unit circle. The Python worker code below performs this estimation operation.

def worker(rg, n, n_reps, param_tuple):
    x = 2 * rg.random((n, n_reps)) - 1
    y = 2 * rg.random((n, n_reps)) - 1

    inside = x*x + y*y <= 1
    pi_hat = 4 * inside.mean(axis = 0)

    return pi_hat

Running sim_template from the IDLE console with this worker function produces results like

import numpy as np
out = sim_template(1234, 100, 1000000, 4)
out

## [array([3.2 , 3.4 , 3.08, ..., 2.92, 3.12, 3.  ], shape=(250000,)),
## array([3.24, 3.2 , 3.04, ..., 3.32, 3.12, 3.44], shape=(250000,)),
## array([3.24, 3.12, 3.16, ..., 3.08, 3.28, 3.28], shape=(250000,)),
## array([3.2 , 3.28, 3.36, ..., 3.08, 3.4 , 3.12], shape=(250000,))]

pi_hat = np.concatenate(out).mean()
pi_hat

## np.float64(3.14141632)

Here we ran sim_template with 4 processes and the output is a 4 component list. Each element of the list is an array of 250000 estimates of \(\pi\) obtained from 100 virtual dart throws. The arrays are concatenated and the mean of all the million estimators is the resulting approximation to \(\pi\). Again, the square-root law from Chapter 2 means we only achieved three decimals of accuracy with a hundred million dart throws.

Some timing experiments that employed the perf_counter timer from the Python time module produced the results in Table 4.1. The table gives the median execution times across 10 runs for computing \(\pi\) estimates with sim_template when the sample size is 100 with a million replications.

Table 4.1: Table 4.2: Median elapsed execution times (seconds) for the parallel \(\pi\) simulation. Interquartile ranges are shown in parentheses. Speedup is computed relative to the one-core execution on the same machine.
Cores	iMac M4 Median (IQR)	Speedup	ASUS i7 Median (IQR)	Speedup
1	1.452 (0.011)	1.000	1.972 (0.011)	1.000
2	1.020 (0.046)	1.422	1.154 (0.005)	1.709
3	0.774 (0.016)	1.875	0.903 (0.004)	2.185
4	0.668 (0.013)	2.174	0.769 (0.008)	2.566
5	0.614 (0.005)	2.365	0.682 (0.018)	2.894
6	0.588 (0.004)	2.470	0.630 (0.007)	3.128
7	0.583 (0.010)	2.491	0.607 (0.005)	3.247
8	0.597 (0.006)	2.432	0.606 (0.009)	3.257

The worker count varies from 1 to 8. The ASUS machine reaches a speedup of about 3.25 with seven or eight processes, while the iMac reaches a maximum speedup of about 2.5 with seven processes.

4.3 Cauchy Data Revisited

Let us now return to the Cauchy scale parameter model of Section 3.3. In that instance data is obtained by sampling from the Cauchy distribution with the objective being estimation of the scale parameter \(\theta\). The estimator we consider is half of the sample interquartile range. The corresponding sim_worker.py code in this instance might look like

import numpy as np

def worker(rg, n, n_reps, param_tuple):
    theta = param_tuple[0]

    U = rg.random((n, n_reps))
    # quantile transformation
    X = theta * np.tan(np.pi * (U - 0.5))

    Q = np.quantile(X, [0.25, 0.75], axis=0)
    theta_hat = (Q[1] - Q[0]) / 2

    return theta_hat

Here a trial or postulated value for the true scale parameter comes in from the *args argument to sim_template and is passed on to worker in theta = param_tuple[0].

Then,

X = theta * np.tan(np.pi * (U - 0.5))

uses the Cauchy quantile function to transform the uniform random numbers in the U array into Cauchy random variates with scale parameter theta. That generates the data. To compute the parameter estimators we use the vectorized (sample) quantile function from the NumPy library to compute the two quartiles and take half their difference as the estimator. Note that axis = 0 applies the function to the columns of X.

This choice of sim_worker produces results like

sim_template(1234, 100, 20, 4, 5)

## [array([5.62082898, 5.9548292 , 5.25999085, 4.61180851, 5.32471512]),
## array([4.69824636, 5.99497178, 6.15112504, 5.02016374, 5.01611561]),
## array([5.35159213, 6.2273778 , 4.5618758 , 5.2565794 , 5.83274694]),
## array([4.80614598, 5.43726919, 4.91259206, 5.42029505, 5.47352615])]

The results correspond to Cauchy scale parameter estimators for 20 samples of size 100 using four processes with the true value of the scale parameter set to 5.

Some timing experiments produced results like those in Table 4.3 that shows the median elapsed execution time across ten runs for sim_template using the Cauchy interquartile range version of sim_worker with samples of size 100, 1,000,000 replications and anywhere from 1 to 8 cores on our 8-core iMac and 20-core Windows ASUS laptop machines. The best speed-up, ratio of serial time to parallel time, is 3.78 for the iMac using 7 or 8 processes and 4.74 for the ASUS using 8.

Table 4.3: Median elapsed execution times (seconds) for the parallel Cauchy simulation. Interquartile ranges are shown in parentheses. Speedup is computed relative to the one-core execution on the same machine.
Cores	iMac M4 Median (IQR)	Speedup	ASUS i7 Median (IQR)	Speedup
1	3.488 (0.021)	1.000	3.236 (0.016)	1.000
2	2.101 (0.038)	1.661	1.747 (0.030)	1.852
3	1.535 (0.016)	2.273	1.299 (0.003)	2.491
4	1.222 (0.008)	2.854	1.037 (0.012)	3.119
5	1.043 (0.005)	3.345	0.889 (0.014)	3.639
6	0.950 (0.013)	3.674	0.781 (0.010)	4.143
7	0.924 (0.014)	3.776	0.714 (0.024)	4.532
8	0.924 (0.014)	3.776	0.683 (0.007)	4.738

4.4 Program Details

Let us now describe how sim_template works. The program begins by bringing in the code for constructing data and computing estimators that resides in the file sim_worker.py. From the numpy.random module the default random number generator is then imported. NumPy’s default generator uses a PCG64 bit generator: c.f. Appendix C.

We want different processors to have their own random number generators that produce number streams that behave as independent sequences with no detectable correlation or overlap. The recommended approach for accomplishing this is to use the SeedSequence method, from the numpy.random module. SeedSequence plays the same role in the Python version of sim_template as that of generateSeedVectors in its R implementation. It uses hashing techniques to transform user-provided seeds into high-quality initial states for BitGenerators, the low level objects NumPy uses to produce raw random bits.

The n_reps samples must be allocated to the workers appropriately. In general, n_reps is not a multiple of n_cores, and there will be additional samples that must be treated by some of the processes. The next step is to load balance and portion out the additional samples to the first few workers. A list set_up of n_cores lists is then populated with a random number generator for each process, the sample size, the number of replications to be treated by each worker and the elements of the tuple args argument.

Finally, before going parallel, a pool of n_cores worker processes is created and starmap from the multiprocessing module is used to apply the worker function from sim_worker.py to set_up. This has the effect that each process in the pool will receive one of the sublists from set_up and that list’s content will be used as the arguments for sim_worker.

Because Python starts worker processes differently across operating systems, executable code that calls sim_template should be protected by an if __name__ == "__main__": guard: e.g.,

from sim_template import sim_template
if __name__ == "__main__":
    result = sim_template(1234, 100, 1000000, 4)

On systems that use the spawn method, each worker imports the main module when it starts. Without the guard, the imported module may call sim_template again, recursively creating new process pools. The guard ensures that the simulation is started only by the original Python process.

4.5 The GIL and Multiprocessing

One aspect of Python that is relevant for parallel computation is the presence of the Global Interpreter Lock (GIL). In a conventional GIL-enabled CPython build, only one thread at a time can execute Python bytecode. Compiled extensions such as NumPy may release the GIL while performing computations outside the interpreter. In other words, on most machines, multiple threads do not run Python code simultaneously.

The primary reason for the GIL is to simplify memory management. It accomplishes that by only allowing one thread to execute at any one moment. Multiple threads simply take turns executing, rather than running concurrently on different cores.

For this reason, we have avoided threads and instead relied on the multiprocessing module. This approach creates separate processes (rather than threads) each of which has its own Python interpreter and GIL thereby allowing computations to proceed in parallel. The Pool object used in sim_template.py manages these worker processes, distributes the simulation tasks, and controls the pool’s lifetime. When execution leaves the with block, Python shuts down the worker processes and releases their associated resources.

4.6 Chapter Summary

In this chapter we developed a Python program sim_template.py for running parallel simulation experiments. We illustrated its use by estimating \(\pi\) and computing interquartile range estimators of a Cauchy scale parameter. The speed-up from using multiple processes was notable.

The basic Python template can be used directly for a wide variety of simulation problems requiring only the choice of a suitable worker function. In Chapter 6 we show how it can be adapted for problems that require node-to-node communication.

The next chapter returns us to R. The point and interval estimation of a regression change point is investigated which provides a case study of how to use the R sim_template for a problem of research substance.

3 Parallel Simulation in R

5 Case Study: Segmented Regression with Bootstrapping in R