Futureverse:
Worry-Free Parallelization in R

Henrik Bengtsson

University of California, San Francisco

R Foundation, R Consortium, Bioconductor

@HenrikBengtsson



RaukR 2025, Visby, Sweden (2025-06-12)

Why do we parallelize code?

Parallel & distributed processing can be used to:

  • speed up processing (wall time)
  • lower memory footprint (per machine)
  • avoid data transfers (compute where data lives)
  • Other reasons, e.g. asynchronous UI

How do we parallelize code?

We may choose to parallelize on:

  • our personal laptop or work desktop computer (single user)
  • a shared powerful computer (multiple users)
  • across many computers, e.g. in the office or in the cloud
  • high-performance compute (HPC) cluster (multiple users) with a job scheduler, e.g. Slurm, Son of Grid Engine (SGE)
  • peer-to-peer (P2P) network of computers

R: A Functional Programming Language

R is a functional programming language

Mathematics:

\[ f(x) = 2 \cdot x \\ y = f(1) \\ z = f(3) \]

R:

f <- function(x) 2*x
y <- f(1)
z <- f(3)

A function has code, arguments, and a value

f <- function(x) {
  2*x
}

A function has:

  • code, e.g. 2*x
  • arguments (“parameters”), e.g. x
  • a value

When we call the function;

y <- f(4)

R evaluates the code with the value of argument x (= 4), returns the value (= 8), and assigns the value to variable y.

A function can take multiple arguments

f <- function(x, b) {
  b*x
}

A function has:

  • code, e.g. b*x
  • arguments (“parameters”), e.g. x and b
  • a value

When we call the function;

y <- f(4, 2)

R evaluates the code with the values of arguments x (= 4) and b (= 2), returns the value, and assigns the value to variable y.

Arguments can have default values

f <- function(x, b = 2) {
  b*x
}

When we call the function;

y <- f(4)

R evaluates the code with the value of argument x (= 4), the default value of b (= 2), returns the value, and assigns the value to variable y.

Summation

\[ h(\mathbf{x}) = \sum_{i=1}^{n} x_i ; \mathbf{x} = (x_1, x_2, \ldots, x_n) \\ \]

h <- function(x) {
  sum <- 0
  n <- length(x)
  for (i in 1:n) sum <- sum + x[i]
  sum
}

Summation - emulate slowness

\[ h(\mathbf{x}) = \sum_{i=1}^{n} x_i ; \mathbf{x} = (x_1, x_2, \ldots, x_n) \\ \]

h <- function(x) {
  sum <- 0
  for (i in seq_along(x)) {
    Sys.sleep(1.0)   # one-second delay
    sum <- sum + x[i]
  }
  sum
}
> h(1:3)  # ~3 seconds
[1] 6

Summation can be done in parts

\[ h(\mathbf{x}) = \sum_{i=1}^{n} x_i \]

\[ = \sum_{i=1}^{k} x_i + \sum_{i=k+1}^{n} x_i \\ \]

\[ = h(\mathbf{x}_{1:k}) + h(\mathbf{x}_{(k+1):n}) \]

Summation can be done in parts

\[ h(\mathbf{x}) = \sum_{i=1}^{n} x_i \\ \mathbf{x} = (1, 2, \ldots, 10) \]

\[ y = h(\mathbf{x}) \]

\[ a = h(\mathbf{x}_{1:5}) \\ b = h(\mathbf{x}_{6:10}) \\ y = a + b \]

\[ y = 55 \]

Summation of parts can be done in any order

\[ ~\\ ~\\ y = h(\mathbf{x}) \]

y <- h(x)

\[ a = h(\mathbf{x}_{1:5}) \\ b = h(\mathbf{x}_{6:10}) \\ y = a + b \]

a <- h(x[1:5])
b <- h(x[6:10])
y <- a + b

\[ b = h(\mathbf{x}_{6:10}) \\ a = h(\mathbf{x}_{1:5}) \\ y = a + b \]

b <- h(x[6:10])
a <- h(x[1:5])
y <- a + b

\[ y = 55 \]

All takes ~10 seconds

Concurrent summation of parts

\[ \text{concurrently} \left\{\begin{matrix} a = h(\mathbf{x}_{1:5}) \\ b = h(\mathbf{x}_{6:10}) \\ \end{matrix}\right. \\ \hspace{10ex} y = a + b \]

How can we calculate \(a\) and \(b\) concurrently in R?

TL;DR

a %<-% h(x[1:5])
b %<-% h(x[6:10])
y <- a + b

a and b takes 5 seconds each => everything takes ~5 seconds

Package ‘future’ adds support for concurrent programming in R


library(future)
plan(multisession)

a %<-% h(x[1:5])
b %<-% h(x[6:10])

y <- a + b


  • %<-% is called the future-assignment operator

Your turn: Setup

Step 1: Install packages

install.packages("futureverse", dependencies = TRUE)

Step 2: Define tic() and toc()

Create functions tic() and toc() to measure time by copy’n’pasting:

tic <- function() {
  tic_start <<- base::Sys.time()
}

toc <- function() {
  dt <- base::difftime(base::Sys.time(), tic_start)
  dt <- round(dt, digits = 1L)
  message(paste(format(dt), "since tic()"))
}

Your turn: Setup

These functions can be used as a timer, e.g.

tic()
Sys.sleep(1.5)
toc()
#> 1.5 secs since tic()
Sys.sleep(4.0)
toc()
#> 5.5 secs since tic()

Your turn: Setup

Step 3: Implement \(h()\) as above, but call it slow_sum() 

slow_sum <- function(x) {
  sum <- 0
  for (value in x) {
    Sys.sleep(1.0)    # one-second slowdown per value
    sum <- sum + value
  }
  sum
}

Your turn: Verify

Step 4: Test slow_sum()

  1. Create \(x = (1, 2, ..., 10)\) in R

  2. Call y <- slow_sum(x) and confirm it takes 10 seconds

  3. Call a <- slow_sum(x[1:5]) and confirm it takes 5 seconds

  4. Same for b <- slow_sum(x[6:10])

Your turn: Our very first concurrent evaluation of R code

Step 5: Sum parts concurrently

library(future)
plan(multisession)

a %<-% slow_sum(x[1:5])
b %<-% slow_sum(x[6:10])

y <- a + b

  1. Confirm correct result

  2. Confirm that it finishes in 5 seconds - use tic() and toc()

  3. How long does each line take? - use tic() and toc()

Concurrent Programming in R

Future … what?

  • A future is an abstraction for a value that will be available later
  • The state of a future is either unresolved or resolved
  • The value is the result of an evaluated expression

An R assignment:

v <- expr

A future-value assignment:

f <- future(expr)
v <- value(f)


y <- slow_sum(x)
f <- future({ slow_sum(x) })
y <- value(f)

References

Friedman & Wise (1976, 1977), Hibbard (1976), Baker & Hewitt (1977)

Your turn: future() and value()

library(future)
plan(multisession)

a %<-% slow_sum(x[1:5])
b %<-% slow_sum(x[6:10])

y <- a + b
library(future)
plan(multisession)

fa <- future( slow_sum(x[1:5])  )
fb <- future( slow_sum(x[6:10]) )

y <- value(fa) + value(fb)

  • Verify result

  • Verify total run time

  • Verify run time per line

Functions

\[ a = \text{slow_sum}(\mathbf{x}_{1:5}) \\ b = \text{slow_sum}(\mathbf{x}_{6:10}) \\ y = a + b \]

a <- slow_sum(x[1:5])
b <- slow_sum(x[6:10])
y <- a + b

Functions and Futures

\[ \text{concurrently} \left\{\begin{matrix} a = \text{slow_sum}(\mathbf{x}_{1:5}) \\ b = \text{slow_sum}(\mathbf{x}_{6:10}) \\ \end{matrix}\right. \\ \hspace{10ex} y = a + b \]

fa <- future( slow_sum(x[1:5])  )
fb <- future( slow_sum(x[6:10]) )
y <- value(fa) + value(fb)

Futureverse - A Friendly, Unifying Parallelization Framework in R

  • “Write once, run anywhere” - a simple unified API
  • 100% cross-platform - easy to install
  • works with any type of parallel backends
  • very well tested, lots of CPU mileage
    10 years, used by 1300 packages, top-0.6% downloaded

Low friction:

  • Automatically exports global variables
  • Automatically relays output, messages, and warnings
  • Proper parallel random number generation (RNG)
  • … and much more

Globals automatically found (99.5% worry free)

Static-code inspection by walking the abstract syntax tree (AST):

x <- rnorm(n = 100)          lobstr::ast( { slow_sum(x) } )
f <- future({ slow_sum(x) })     |  █─`{`
                   └─────────────|  └─█─ slow_sum
                                 |    └─ x

=> globals identified and exported to the worker:

  • slow_sum() - a function (also searched recursively)
  • x - a numeric vector of length 100

Globals can also be specified manually;

f <- future({ slow_sum(x) }, globals = c("slow_sum", "x"))

Apply the same function over and over

For-loop:

xs <- list(1:25, 26:50, 51:75, 76:100)
ys <- list()
for (i in seq_along(xs)) {
  ys[[i]] <- slow_sum(xs[[i]])
}


for (x in xs) {
  ys <- c(ys, slow_sum(x))
}

Apply the same function over and over

Base R:

xs <- list(1:25, 26:50, 51:75, 76:100)
ys <- lapply(xs, slow_sum)

Tidyverse:

library(purrr)
xs <- list(1:25, 26:50, 51:75, 76:100)
ys <- map(xs, slow_sum)

Apply the same function over and over

Plyr:

library(plyr)
xs <- list(1:25, 26:50, 51:75, 76:100)
ys <- llply(xs, slow_sum)

Foreach:

library(foreach)
xs <- list(1:25, 26:50, 51:75, 76:100)
ys <- foreach(x = xs) %do% slow_sum(x)

Pipes: Syntactic sugar to spice it up

Base R:

ys <- lapply(xs, slow_sum)

is identical to

ys <- xs |> lapply(slow_sum)

R can’t tell the difference! (e.g. try with lobstr::ast(...))

Tidyverse:

ys <- xs |> map(slow_sum)

Plyr:

ys <- xs |> llply(xs, slow_sum)

Your turn: Draw an Owl

ys <- xs |> map(slow_sum)
ys <- xs |> future_map(slow_sum)

Your turn: Start with a for-loop

Combine the idea of:

f <- future( slow_sum(x) )
v <- value(f)

with:

ys <- list()
for (i in seq_along(xs)) {
  ys[[i]] <- slow_sum(xs[[i]])
}

to make a “concurrent” for-loop.

Hint: You need two for-loops - one for future() and one for value()

Your turn: Concurrent for-loop to concurrent map

# Create futures
fs <- list()
for (i in seq_along(xs)) {
  fs[[i]] <- future( slow_sum(xs[[i]]) )
}

# Collect values
ys <- list()
for (i in seq_along(xs)) {
  ys[[i]] <- value( fs[[i]] ) 
}

Your turn: A concurrent map()

# Create futures
fs <- purrr::map(xs, function(x) { future(slow_sum(x)) })

# Collect values
ys <- purrr::map(fs, value)

value() is vectorized, i.e. it can operate on lists too

# Create futures
fs <- purrr::map(xs, function(x) { future(slow_sum(x)) })

# Collect values
ys <- value(fs)

ys <- future_map(xs, slow_sum)

User-friendly, higher-level functions

  • The concept of “futures” was invented in 1975 (sic!)
  • future(), value(), and resolved() are easy to understand, powerful constructs

These building blocks lay the foundation for higher-level functions:

  • future.apply, e.g. future_lapply() and future_replicate()
  • furrr, e.g. future_map() and future_map_dbl()
  • doFuture, e.g. foreach() %dofuture% { ... }
  • …
  • - Maybe your package will be next?

Futureverse allows you to stick with your favorite coding style

Parallel alternatives to traditional, sequential functions:

ys <- lapply(xs, slow_sum)                    # base R
ys <- future_lapply(xs, slow_sum)             # {future.apply}

ys <- map(xs, slow_sum)                       # {purrr}
ys <- future_map(xs, slow_sum)                # {furrr}

ys <- foreach(x = xs) %do% slow_sum(x)        # {foreach}
ys <- foreach(x = xs) %dofuture% slow_sum(x)  # {foreach}+{doFuture}

Yes, we can of course use base-R or magrittr pipes where we want to, e.g.

ys <- xs |> future_map(slow_sum)
ys <- xs %>% future_map(slow_sum)

Futures can be resolved anywhere!

Our ys <- future_map(xs, slow_sum) can be resolved on lots of different backends - it doesn’t matter where!

On the current machine

plan(sequential)

sequentially (default)

plan(multisession)
plan(multisession, workers=4)

parallelly via background R workers

plan(mirai_multisession)
plan(mirai_multisession, workers=4)

parallelly via future.mirai; low latency

plan(callr)
plan(callr, workers=4)

parallelly via future.callr; all memory is returned when as each future is resolved

Futures can be resolved anywhere!

On other machines

plan(cluster, workers=c("n1", "m2.uni.edu", "vm.cloud.org"))

distributed on other local or remote computer

plan(mirai)

distributed via future.mirai daemons queue

On high-performance compute (HPC) cluster

plan(batchtools_lsf)
plan(batchtools_openlava)
plan(batchtools_sge)
plan(batchtools_slurm)
plan(batchtools_torque)

parallelly on HPC job schedulers (LSF, OpenLava, TORQUE/PBS, Grid Engine SGE/AGE, Slurm); for long-running tasks; high latency

Your turn: Try different backends

Sequential:

library(future)
plan(sequential)  # default

xs <- list(1:5, 6:10, 11:15, 16:20)
ys <- future_map(xs, slow_sum)

Built-in multisession:

plan(multisession)
ys <- future_map(xs, slow_sum)

Add-on mirai-based multisession:

plan(future.mirai::mirai_multisession)
ys <- future_map(xs, slow_sum)

Messages, Warnings, and Errors




     

Some functions have side effects

g <- function(x, b = 2) {
  if (b <= 0) stop("Argument 'b' must be a positive scalar: ", b)
  sqrt(b * x)
}

What type of results may this function produce?

  • a value
> g(1)
[1] 2
  • an error
> g(1, b = 0)
Error in g(1, b = 0) :
Argument 'b' must be a positive scalar: 0
  • a value and a warning
> g(-1)
[1] NaN
Warning message:
In sqrt(b * x) : NaNs produced

Your turn: Messages, warnings, and errors with futures

g <- function(x, b = 2) {
  message("x = ", x)
  if (b <= 0) stop("Argument 'b' must be a positive scalar: ", b)
  Sys.sleep(3.0)
  sqrt(b * x)
}

f <- future({ g(1) })
v <- value(f)

f <- future({ g(1, b = 0) })
v <- value(f)

f <- future({ g(-1) })
v <- value(f)

Your turn: Errors terminate calls

xs <- c(1, 2, 3, 4)
bs <- c(4, 3, 2, 1)
ys <- purrr::map2(xs, bs, g)
g <- function(x, b = 2) {
  message("x = ", x)
  if (b <= 0)
    stop("Argument 'b' must be positive: ", b)
  Sys.sleep(3.0)
  sqrt(b * x)
}

xs <- c(1, 2, -3, 4)
bs <- c(4, 3,  2, 1)
ys <- purrr::map2(xs, bs, g)

xs <- c(1, 2,  3, 4)
bs <- c(4, 3, -2, 1)
ys <- purrr::map2(xs, bs, g)

Your turn: Errors cancel remaining futures

library(future)
plan(multisession, workers = 4)

xs <- c(1, 2, 3, 4)
bs <- c(4, 3, 2, 1)
ys <- furrr::future_map2(xs, bs, g)

xs <- c(1, 2, -3, 4)
bs <- c(4, 3,  2, 1)
ys <- furrr::future_map2(xs, bs, g)

xs <- c(1, 2,  3, 4)
bs <- c(4, 3, -2, 1)
ys <- furrr::future_map2(xs, bs, g)

Business as usual but in parallel

library(future)
plan(multisession, workers = 4)
xs <- c(1, 2, -3, 4)
bs <- c(4, 3, 2, 1)

> ys <- furrr::future_map2(xs, bs, g)
x = 1
x = 2
x = -3
x = 4
Warning message:
In sqrt(b * x) : NaNs produced
>

Objective: Minimal mental load

library(future)
plan(multisession, workers = 4)
xs <- c(1, 2, -3, 4)
bs <- c(4, 3, 2, 1)

> ys <- purrr::map2(xs, bs, g) |> hey_r_this_is_concurrent()
x = 1
x = 2
x = -3
x = 4
Warning message:
In sqrt(b * x) : NaNs produced
>

Things we haven’t talked about

foreach() is not a for loop!

The name foreach() tricks us …

The name foreach() tricks us to believe it’s a for-loop

fq <- fs::dir_ls(glob = "*.fq")
bam <- list()
for (ii in seq_along(fq)) {
  bam[[ii]] <- align(fq[[ii]])
}

We must never think of foreach() as a for-loop

The name foreach() tricks us to believe it’s a for-loop

library(foreach)
fq <- fs::dir_ls(glob = "*.fq")
bam <- list()
foreach(ii = seq_along(fq)) %dopar% {
  bam[[ii]] <- align(fq[[ii]])
}

This does not work because … foreach() is not a for-loop!

Repeat after me: foreach() is not a for-loop!


for (ii in 1:1000) {
  message("foreach() is not a for-loop!")
  Sys.sleep(2)
  beepr::beep("wilhelm")
}

foreach() is just like lapply() …


fq <- fs::dir_ls(glob = "*.fq")
bam <- list()
lapply(seq_along(fq), function(ii) {
  bam[[ii]] <- align(fq[[ii]])
})

This does not work, because the <- assignment is done inside a function

… and just like map() …


fq <- fs::dir_ls(glob = "*.fq")
bam <- list()
purrr::map(seq_along(fq), function(ii) {
  bam[[ii]] <- align(fq[[ii]])
})

This does not work, because the <- assignment is done inside a function

If foreach() had looked like …


fq <- fs::dir_ls(glob = "*.fq")
bam <- list()
foreach(seq_along(fq), function(ii) {
  bam[[ii]] <- align(fq[[ii]])
})

… it would be clear that the <- assignment is done inside a function.

Takehome: lapply(), map(), foreach() return values


fq <- fs::dir_ls(glob = "*.fq")

bam <- lapply(fq, align)

bam <- purrr::map(fq, align)

bam <- foreach(x = fq) %dopar% align(x)

Nested parallelization with some care

  • “Structured concurrency” -
    Futureverse protects against over-parallelization

Sequential
processing

Parallelization
with 4 workers

Nested parallelization
with 5 workers each
running 3 parallel tasks

HPC processing: future.batchtools


plan(future.batchtools::batchtools_slurm)                          ‎

fq <- fs::dir_ls(glob = "*.fq")        # 80 files; 200 GB each
bam <- furrr::future_map(fq, align)    # 1 hour each


{henrik: ~}$ squeue
Job ID   Name            User         Time Use S
-------- --------------- ------------ -------- -
606411   xray            alice        46:22:22 R
606638   future_map-5    henrik       00:52:05 R
606641   python          bob          37:18:30 R
606643   future_map-6    henrik       00:51:55 R
...

Making progress . . .




slides

Together we can make a bright future!

Future turns 10 years next week!

Series of exciting feature announcements at:

https://www.futureverse.org/blog.html