36-651/751: Parallelism in Practice

Python and multiprocessing

import multiprocessing

NUM_TRIALS = 500

def evaluator_on_data(data):
    def fit_and_evaluate(tuning_params):
        for _ in range(NUM_TRIALS):
            ## do stuff...

        return test_error

    return fit_and_evaluate

all_tuning_params = range(...)

## by default, makes as many processes as there are CPUs
pool = multiprocessing.Pool()

test_errors = pool.map(evaluator_on_data(my_dataset), all_tuning_params)

pool.map can take arbitrary functions and arbitrary sequences, and will automatically divide the sequence of data up to be given to different processes in the pool. This does have caveats:

1. The function (the first argument) must be pickle-able, because multiprocessing pickles the function to pass it to the other processes. (Remember, processes don’t share memory, so the function has to be sent to them.) Unfortunately, while functions can be pickled, object methods cannot be.
2. The data must also be pickle-able.

3. Global variables are problematic. The NUM_TRIALS variable is available in all the pool processes. Why? fork! We can read its value without causing problems, but note that if we modify it inside fit_and_evaluate, each process will have its own value, not a shared value.

   But this doesn’t work for every kind of variable.

   For example, suppose you’re using PostgreSQL and write

   import psycopg2
   import multiprocessing
   
   conn = psycopg2.connect(database = "some_database", user = "username",
                           password = "password", host = "localhost")
   cur = conn.cursor()
   
   pool = multiprocessing.Pool()
   
   results = pool.map(some_function_using_cur, some_big_dataset)
   

   This will not work correctly, because PostgreSQL cursor objects are not meant to be shared. You need to create one cursor per process.

   Also be careful if you have open files or datasets and try to write to them from multiple processes.

4. Global variables are also nice. Suppose your parallel involves some 10GB dataset; you want to run your analysis on 8 cores, but you do not have 80GB of RAM. If you load the dataset before creating the Pool, that’s fine: as long as the worker processes do not modify it, only read from it, they will not create their own copies of it.

R and parallel

In R, the parallel package works in a similar way. Instead of creating a pool of processes, you make a “cluster”, then submit tasks to the cluster for computation. The cluster can just be a bunch of processes on one computer – created via fork, for example – or it can consist of R processes on many different machines, which parallel connects to via SSH.

Let’s start with the simple case. We can make a fork cluster and submit tasks to it:

library(parallel)

## argument is the number of processes to make. detectCores() documentation
## actually suggests you don't do this, because it can return NA if it doesn't
## know the number of cores, and can return smaller numbers if e.g. your laptop
## power-saving settings mean it wants to use fewer cores than it actually has.
## Also it may be impolite to use every core on a shared machine.
cl <- makeForkCluster(detectCores())


result <- clusterApply(cl, big_dataset, analyze_data_fn)

## or:

result <- parLapply(cl, big_dataset, analyze_data_fn)
result <- parSapply(cl, big_dataset, analyze_data_fn)
## ... and so on.

stopCluster(cl) # kill the cluster processes

Note that fork clusters do not work on Windows, because there is no fork on Windows.

On a Mac, forking causes problems with graphics: the Quartz graphics window can only be controlled from one process at a time, so it will be disabled in the children.

We can alternately make a PSOCK cluster, which creates new R processes not by forking but by creating Rscript processes from scratch and sending them code to evaluate over sockets.

library(parallel)

cl <- makePSOCKcluster(c("hydra1.stat.cmu.edu", "hydra2.stat.cmu.edu"),
                       user="your SSH username", homogeneous=TRUE,
                       rscript="/usr/bin/Rscript")

## then do stuff as usual

The rscript argument tells R where to find the Rscript program on the servers; otherwise it tries to find it in the same place it exists on your computer, which on Macs is a bizarre location that does not exist on Linux.

This will create one process per host; repeat the hostname in the first argument if you want to run multiple processes per server. Please do not use all eight cores on every hydra.

Note that if you have to send large datasets for computation, R will have to take the time to send those datasets to each server in the cluster. You might instead arrange to have the data in files on each server (e.g. in your home directory) and have the R processes load the data from there.

Be careful not to have each worker output to the same file and cause conflicts.

The threaded way, with OpenMP

OpenMP is a framework for multithreaded programming that’s largely intended to parallelize different kinds of for loops. It works by being a part of the C/C++/Fortran compiler and reading a specially formatted comment above the loop that tells OpenMP how it should be parallelized.

Because it works for C and C++, it’s available in Rcpp (see the example gallery here) and in Cython (see the documentation here). It requires passing an extra flag to the compiler, usually -fopenmp. If you forget to pass the flag, the compiler ignores the OpenMP directives and compiles everything to run single-threaded.

Let’s use an example from this article. Suppose you’re doing a calculation that involves heavy use of the multivariate normal density; you use the dmvnorm function from the mvtnorm package, but it’s slow. Profiling reveals that it’s a major part of the execution time of your code.

You investigate your code and see that you calculate the density for many points, and you could parallelize these calculations.

Here’s the implementation in Rcpp, using the RcppArmadillo package for its linear algebra:

// [[Rcpp::depends("RcppArmadillo")]]
// [[Rcpp::export]]
arma::vec dmvnrm_arma_mc(arma::mat x,
                         arma::rowvec mean,
                         arma::mat sigma,
                         bool logd = false,
                         int cores = 1) {
    omp_set_num_threads(cores);

    int n = x.n_rows;
    int xdim = x.n_cols;

    arma::vec out(n);
    arma::mat rooti = arma::trans(arma::inv(trimatu(arma::chol(sigma))));

    double rootisum = arma::sum(log(rooti.diag()));
    double constants = -(xdim/2) * log2pi;

    #pragma omp parallel for schedule(static)
    for (int i = 0; i < n; i++) {
        arma::vec z = rooti * arma::trans( x.row(i) - mean) ;
        out(i)      = constants - 0.5 * arma::sum(z%z) + rootisum;
    }

    if (!logd) {
        out = exp(out);
    }
    return(out);
}

Notice the #pragma. We’re telling OpenMP to run the loop in parallel on a “static” schedule, meaning it decides in advance how to divide up the data between threads. (There are other schedules, like “dynamic”, which gives data to each thread as it finishes its work, so if some loop iterations are much faster than others, the thread with those chunks of the data gets an opportunity to get more data to work on.)

There are other options that can be included in the #pragma line. There are a bunch – review the very detailed guide here – so I’ll just cover a few. Some can be used to declare whether variables inside the loop should be shared or private:

int a, b = 0;

#pragma omp parallel for private(a) shared(b)
for(a = 0; a < 50; ++a) {
  #pragma omp atomic
  b += a;
}

Here a is private – each thread gets its own independent copy – and b is shared between cores, with a note to OpenMP requesting it to use an atomic instruction to do the addition. There are other options for determining what the value of a should be after the loop; e.g. lastprivate(a) says that after the loop, a should have the value it had in the last loop iteration.

OpenMP also understands certain kinds of reductions, including +, i, and *:

int factorial(int number) {
  int fac = 1;

  #pragma omp parallel for reduction(*:fac)
  for(int n = 2; n <= number; ++n) {
    fac *= n;
  }

  return fac;
}

This calculates the factorial in parallel. (Obviously this isn’t terribly useful, because when number is large enough for the parallelization to help, fac will overflow the maximum size of an integer anyway.)

You can have loops with multiple reduction variables. You can also manually construct parallel loops with atomic operations, SIMD, and parallel chunks that aren’t loops (e.g. you can ask for the same chunk of code to be run on every thread).

Cython works by spitting out C or C++ code from the Cython code you write, and it supports OpenMP as well. For example:

def sum(double[:] big_array):
    cdef double s = 0.0

    for ii in prange(big_array.shape[0], nogil=True, schedule='static'):
        s += big_array[ii]

    return s

prange tells Cython that this is a parallel range, and so it generates a for loop with the #pragma attached. Cython automatically recognizes that s is a reduction variable and generates the appropriate #pragma.