Multiprocessing in Julia

The easiest way to start using multiprocessing in Julia is to use the Distributed.jl package. This is the standard package to facilitate the multiprocessing paradigm. To add the package, type the following into the REPL:

using Pkg; Pkg.add("Distributed");

As with multithreading, we should check how many processes are available to do work. This is done with the command:

using Distributed;
nprocs()

By default, Julia launches only a single process when running. There are a few ways in which we can use more, the simplest is to use the addprocs function:

n_processes = 4;
addprocs(n_processes);

This will spawn n_processes locally to act as workers. The total number of processes reported by nprocs() will be one higher than the number of workers, as this includes the master control process. These notes were built with one worker started per available thread rather than the four above, so the numbers below will differ from what you see on your own machine:

@show nprocs()
@show nworkers()

nprocs() = 9
nworkers() = 8

Note: One can start Julia with multiple processes using julia -p 3, which will start Julia with 4 total processes. Be aware that one should also run the following, to make sure that each worker has access to all the packages:

julia Copy

@everywhere begin
    using Pkg
    Pkg.activate(".")
end

The @everywhere macro executes the command on all workers. Note the begin ... end block: @everywhere only applies to the single expression that follows it, so @everywhere using Pkg; Pkg.activate(".") would activate the project on the master process only.

The bread and butter of Distributed.jl are the following macros and functions:

• addprocs(n) - Spawns n available processes. This command launches processes locally by default, but can also be used to launch processes on other machines via SSH.
• @everywhere - A macro that runs the subsequent expression on all connected processes. This is most useful for using statements and include calls to load function definitions on the workers.
• pmap - A multiprocessing, parallel implementation of the standard map function in Julia. This enables easy use of multiprocessing for an embarrassingly parallel problem.
• @distributed - A macro which is useful for parallelising a for loop, with an optional, built-in, mechanism for performing reductions.
• @sync - A macro that halts progression until all enclosed uses of @async, @spawn, @spawnat and @distributed are complete. Note that a @distributed loop which specifies a reduction operator already blocks until it can return the reduced value, so @sync is only strictly needed for the non-reducing form.

Chunked Implementation

Let’s start off with implementing a chunked multiprocessing version of the Monte Carlo pi estimation. We can start off by having an array of numbers representing the number of darts to throw on each worker. Once this array is calculated we can map those numbers onto an estimate of $\pi$.

function est_pi_mc_serial(n)
    n_c = zero(typeof(n))
    for _ in 1:n
        x = rand() * 2 - 1
        y = rand() * 2 - 1
        r2 = x*x+y*y
        if r2 <= 1
            n_c += 1
        end
    end
    return 4 * n_c / n
end

@everywhere function est_pi_mc_serial_worker(n)
    n_c = zero(typeof(n))
    for _ in 1:n
        x = rand() * 2 - 1
        y = rand() * 2 - 1
        r2 = x*x+y*y
        if r2 <= 1
            n_c += 1
        end
    end
    return 4 * n_c / n
end

function est_pi_mc_multiprocessing_chunked(n)
    num_workers = nworkers()

    chunk_size = div(n, num_workers, RoundUp)
    iter = Iterators.partition(1:n, chunk_size)
    # Create a mapping function from a range to an estimate
    _f(range) = est_pi_mc_serial_worker(length(range)) * length(range) / 4
    num_inside = pmap(_f, iter)

    n_c = sum(num_inside)
    return 4 * n_c / n
end
nothing

Note that when using multiprocessing, the function est_pi_mc_serial_worker does not exist on all processes by default. This is a common error when people use the multiprocessing paradigm, as we have to load the function definitions on each worker before they are used. We used @everywhere to define the function on all workers. In larger projects, you can abstract away all the functions in a separate file/module and load it with:

@everywhere include("support_code.jl");

Now that we have the code loaded, we can benchmark the multiprocessing approach:

n = 100_000
mc_pi_serial_time = @belapsed est_pi_mc_serial($n)

function est_pi_mc_threaded_chunked(n)
    n_threads = Threads.nthreads()
    num_inside = zeros(Float64, n_threads)
    chunk_size = div(n, n_threads, RoundUp)
    iter = collect(enumerate(Iterators.partition(1:n, chunk_size)))
    Threads.@threads for info in iter
        i, idx_range = info
        n_block = length(idx_range)
        pi_est = est_pi_mc_serial(n_block)
        num_inside[i] = pi_est*n_block/4
    end
    n_c = sum(num_inside)
    return 4 * n_c / n
end

mc_pi_threaded_chunked_time = @belapsed est_pi_mc_threaded_chunked($n)
mc_pi_mp_chunked_time = @belapsed est_pi_mc_multiprocessing_chunked($n)
println("MP Speedup vs Serial: ", mc_pi_serial_time/mc_pi_mp_chunked_time)
println("Threaded/MP Relative Speed: ", mc_pi_threaded_chunked_time/mc_pi_mp_chunked_time)

MP Speedup vs Serial: 0.44430268718602983
Threaded/MP Relative Speed: 0.07945937325572124

All three implementations throw exactly the same number of darts, so any difference between them is pure overhead. At this value of n that overhead dominates: the multiprocessing version is actually a little slower than the plain serial version, while the multithreaded version — splitting the same work into the same number of chunks — is close to an order of magnitude faster than it. The work handed to each worker here is only a fraction of a millisecond, which is not enough to hide the cost of shipping a task out to another process and its result back again. One should remember that latency costs are much higher in multiprocessing than in multithreading, and a multiprocessing implementation only starts to pay for itself once each chunk carries substantially more work than this.

Using @distributed

Since we are effectively performing a map of empty arguments into the sum of darts in and out, we can look at an alternative implementation using the @distributed macro:

function est_pi_mc_multiprocessing(n)
    n_c = @sync @distributed (+) for _ = 1:n
        # Choose random numbers between -1 and +1 for x and y
        x = rand() * 2 - 1
        y = rand() * 2 - 1
        # Work out the distance from origin using Pythagoras
        r2 = x*x+y*y
        # Count point if it is inside the circle (r^2=1)
        r2 <= 1 # Last line indicates term to reduce
    end
    return 4 * n_c / n
end

mc_pi_mp_time = @belapsed est_pi_mc_multiprocessing($n)
println("MP Simple/Chunked Relative Speed: ", mc_pi_mp_chunked_time/mc_pi_mp_time)

MP Simple/Chunked Relative Speed: 1.3175741339006646

We see that this implementation is not only much simpler, but is also somewhat faster than the chunked implementation. This is because @distributed splits the range across the workers for us, so each worker still runs one long, cache-friendly local loop, exactly as in the chunked version — the inner loops are doing the same work. What differs is how that work is dispatched: pmap serialises a task and a result for each chunk through its dynamic scheduler, whereas @distributed performs a single static split and folds the partial results in with the reduction operator as each worker finishes. Reach for the manual chunking only when you need control over how the work is divided; for a straightforward reduction over a range, @distributed is both shorter and quicker.

Suggested Multiprocessing Pattern

During development, one may switch between using serial, multithreading and multiprocessing patterns. For this reason, it is often best to write the inner loop of an algorithm in a separate function. This avoids having to maintain several similar copies of implementations of the bulk of your code. For this reason, we suggest developing your code in the following way:

1. Write all inner loop functions, e.g. a single run of a Monte-Carlo simulation, inside a single (or small number of) file(s). Make sure you have a single file (which could simply include various other files), which loads all function definitions of any functions that may be run in parallel.
2. Have a separate file for executing your code, which runs all the @everywhere macros. This file should run @everywhere include("allfunctions.jl"), where “allfunctions.jl” is the relative path to the file with all function definitions that are needed by the parallel code.
3. Make use of packages like Transducers.jl or Folds.jl.

General Parallel Pattern

As a quick example, let’s write a function which emulates pmap with a custom flag to indicate whether to use multithreading, multiprocessing or serial implementations:

# Declare an enum for the different map types
abstract type AbstractExecutionMethod end
struct MultiprocessingEx <: AbstractExecutionMethod end
struct MultithreadingEx <: AbstractExecutionMethod end
struct SerialEx <: AbstractExecutionMethod end

custom_map(::MultiprocessingEx, mapping_fn, c...) = pmap(mapping_fn, c...)
function custom_map(::MultithreadingEx, mapping_fn, c...)
    # Do some work to infer return type of function
    return_types = Base.return_types(
        mapping_fn,
        Tuple{eltype.(c)...}
        )

    if length(return_types) == 1
        return_type = return_types[begin]
    else
        return_type = Union{return_types...}
    end
    # Create a container to store results
    container = Array{return_type}(undef, size(first(c)))
    Threads.@threads for i in eachindex(container)
        container[i] = mapping_fn(map(x -> x[i], c)...)
    end
    return container
end
custom_map(::SerialEx, mapping_fn, c...) = map(mapping_fn, c...)

# Change the default implementation in your code
custom_map(mapping_fn, c...) = custom_map(SerialEx(), mapping_fn, c...)
nothing

We can then implement our naive Monte-Carlo algorithm with a single implementation:

function throw_dart()
    x = rand() * 2 - 1
    y = rand() * 2 - 1
    r2 = x*x+y*y
    return r2 <= 1
end

function throw_darts(n)
    total = zero(typeof(n))
    for _ in 1:n
        total += throw_dart()
    end
    return total
end

@everywhere function throw_dart_w()
    x = rand() * 2 - 1
    y = rand() * 2 - 1
    r2 = x*x+y*y
    return r2 <= 1
end

@everywhere function throw_darts_w(n)
    total = zero(typeof(n))
    for _ in 1:n
        total += throw_dart_w()
    end
    return total
end

get_blocks(::SerialEx) = 1
get_blocks(::MultithreadingEx) = Threads.nthreads()
get_blocks(::MultiprocessingEx) = nworkers()
function estimate_pi(parallel_type, n)
    n_blocks = min(n, get_blocks(parallel_type))
    chunk_size = div(n, n_blocks, RoundUp)
    iter = Iterators.partition(1:n, chunk_size)
    num_darts = [length(r) for r in iter]
    if parallel_type isa MultiprocessingEx
        dart_results = pmap(throw_darts_w, num_darts)
    elseif parallel_type isa MultithreadingEx
        dart_results = Vector{Int}(undef, length(num_darts))
        Threads.@threads for i in eachindex(num_darts)
            dart_results[i] = throw_darts(num_darts[i])
        end
    else
        dart_results = map(throw_darts, num_darts)
    end
    return 4 * sum(dart_results) / n
end

n = 1_000_000
println("Serial:")
@btime estimate_pi($(SerialEx()), $n)
println("Multithreaded:")
@btime estimate_pi($(MultithreadingEx()), $n)
println("Multiprocessing:")
@btime estimate_pi($(MultiprocessingEx()), $n)

Serial:
  4.198 ms (8 allocations: 192 bytes)
Multithreaded:
  626.960 μs (59 allocations: 3.45 KiB)
Multiprocessing:
  1.044 ms (873 allocations: 32.88 KiB)

These implementations are not always optimal, and one can specialise an implementation that can be improved using a different method. For example, if the input is smaller than a certain size, one will always use the serial version of the code.