MPI for parallelization in R

workers <- 10
NCPU <- workers+1
Njob <- 100
work.id.vec <- 1:Njob

job_sleep <- function(work.id){
  library(data.table)
  requireNamespace("pbdMPI") #load to get accurate timing for first job.
  start <- Sys.time()
  my_rank <- tryCatch({
    pbdMPI::comm.rank()
  }, error=function(e)NA_integer_)
  Sys.sleep(work.id/100)
  data.table(my_rank, work.id, start, end=Sys.time(), proc=factor(Sys.getpid()))
}

save(work.id.vec, job_sleep, file="job_sleep.RData")
timing_tables <- list()
diff_seconds <- function(...)as.numeric(difftime(..., units="s"))
polygon.mine <- function
### Make a Positioning Method that places non-overlapping speech
### polygons at the first or last points.
(top.bottom.left.right,
### Character string indicating what side of the plot to label.
  offset.cm=0.1,
### Offset from the polygon to the most extreme data point.
  padding.cm=0.05,
### Padding inside the polygon.
  custom.colors=NULL
### Positioning method applied just before draw.polygons, can set
### box.color and text.color for custom colors.
){
  if(is.null(custom.colors)){
    custom.colors <- directlabels::gapply.fun({
      rgb.mat <- col2rgb(d[["colour"]])
      d$text.color <- with(data.frame(t(rgb.mat)), {
        gray <- 0.3*red + 0.59*green + 0.11*blue
        ifelse(gray/255 < 0.5, "white", "black")
      })
      d
    })
  }
  opposite.side <- c(
    left="right",
    right="left",
    top="bottom",
    bottom="top")[[top.bottom.left.right]]
  direction <- if(
    top.bottom.left.right %in% c("bottom", "left")
  ) -1 else 1
  min.or.max <- if(
    top.bottom.left.right %in% c("top", "right")
  ) max else min
  if(top.bottom.left.right %in% c("left", "right")){
    min.or.max.xy <- "x"
    qp.target <- "y"
    qp.max <- "top"
    qp.min <- "bottom"
    padding.h.factor <- 2
    padding.w.factor <- 1
    limits.fun <- ylimits
    reduce.method <- "reduce.cex.lr"
  }else{
    min.or.max.xy <- "y"
    qp.target <- "x"
    qp.max <- "right"
    qp.min <- "left"
    padding.h.factor <- 1
    padding.w.factor <- 2
    limits.fun <- directlabels::xlimits
    reduce.method <- "reduce.cex.tb"
  }
  list(
    hjust=0.5, 
    function(d,...){
      ## set the end of the speech polygon to the original data point.
      for(xy in c("x", "y")){
        extra.coord <- sprintf(# e.g. left.x
          "%s.%s", opposite.side, xy)
        d[[extra.coord]] <- d[[xy]]
      }
      ## offset positions but do NOT set the speech polygon position
      ## to the min or max.
      d[[min.or.max.xy]] <- d[[min.or.max.xy]] + offset.cm*direction
      d
    },
    "calc.boxes",
    reduce.method,
    function(d, ...){
      d$h <- d$h + padding.cm * padding.h.factor
      d$w <- d$w + padding.cm * padding.w.factor
      d
    },
    "calc.borders",
    function(d,...){
      do.call(rbind, lapply(split(d, d$y), function(d){
        directlabels::apply.method(directlabels::qp.labels(
          qp.target,
          qp.min,
          qp.max,
          directlabels::make.tiebreaker(min.or.max.xy, qp.target),
          limits.fun), d)
      }))
    },
    "calc.borders",
    custom.colors,
    "draw.polygons")
}
my.polygons <- list(cex=0.7, vjust=0, polygon.mine(
  "top", offset.cm=0.2, custom.colors=list(box.color="white")))
make_gg <- function(package, expr){
  before <- Sys.time()
  result <- eval(expr)
  if(!is.data.table(result))result <- rbindlist(result)
  after <- Sys.time()
  DT <- data.table(before, after, result)
  timing_tables[[package]] <<- DT
  library(ggplot2)
  rect_dt <- unique(DT[, .(before, after)])
  seconds_possible <- workers*rect_dt[, diff_seconds(after, before)]
  seconds_used <- sum(DT[, diff_seconds(end, start)])
  ggplot()+
    ggtitle(sprintf("%s efficiency = %.1f%%", package, 100*seconds_used/seconds_possible))+
    theme_bw()+
    geom_rect(aes(
      xmin=before, xmax=after,
      ymin = -Inf, ymax=Inf),
      data=rect_dt,
      alpha=0.5,
      fill="grey")+
    geom_segment(aes(
      start, proc,
      xend=end, yend=proc),
      data=DT)+
    geom_point(aes(
      start, proc),
      shape=1,
      data=DT)+
    directlabels::geom_dl(aes(
      start, proc,
      label.group=paste(work.id, proc),
      label=work.id),
      data=DT,
      method=my.polygons)+
    scale_x_continuous("Time (seconds)")
}
library(data.table)

mirai::daemons(NULL)

mirai::daemons(workers)

make_gg("mirai", mirai::mirai_map(work.id.vec, job_sleep)[])

future::plan("multisession", workers=workers)
make_gg("future", future.apply::future_lapply(work.id.vec, job_sleep))

script_sleep <- function(prefix){
  code <- gsub("SCRIPT", prefix, 'library(data.table)
load("SCRIPT_sleep.RData")
ret_dt_list <- pbdMPI::task.pull(work.id.vec, SCRIPT_sleep)
ret_dt <- rbindlist(ret_dt_list)
if(pbdMPI::comm.rank() == 0)fwrite(ret_dt, "SCRIPT_sleep.csv")
pbdMPI::finalize()
')
  sleep.R <- paste0(prefix, "_sleep.R")
  cat(code, file=sleep.R)
  cmd <- sprintf("mpiexec -np %d Rscript --vanilla %s", NCPU, sleep.R)
  sleep.csv <- paste0(prefix, "_sleep.csv")
  unlink(sleep.csv)
  system(cmd)
  fread(sleep.csv, colClasses=list(factor="proc"))
}

make_gg("pbdMPI", script_sleep("job"))

add_seconds <- function(DT, sec.cols = c("start","end","before","after"))DT[
, paste0(sec.cols, ".seconds") :=
    lapply(sec.cols, function(x)diff_seconds(get(x), before))
]
(compare_dt <- data.table(package=names(timing_tables))[
, add_seconds(timing_tables[[package]])
, by=package])

rect_dt <- unique(compare_dt[, .(package, before.seconds, after.seconds)])
ggplot()+
  theme_bw()+
  geom_blank(aes(
    -Inf, ""),
    data=unique(compare_dt[, .(package)]))+
  geom_rect(aes(
    xmin=before.seconds, xmax=after.seconds,
    ymin = -Inf, ymax=Inf),
    data=rect_dt,
    alpha=0.5,
    fill="grey")+
  geom_segment(aes(
    start.seconds, proc,
    xend=end.seconds, yend=proc),
    data=compare_dt)+
  geom_point(aes(
    start.seconds, proc),
    shape=1,
    data=compare_dt)+
  directlabels::geom_dl(aes(
    start.seconds, proc,
    label.group=paste(work.id, proc),
    label=work.id),
    data=compare_dt,
    method=my.polygons)+
  facet_grid(package ~ ., labeller=label_both, scales="free")+
  scale_y_discrete("process")+
  scale_x_continuous(
    "Seconds since line before apply",
    limits=c(0, NA))

rank_sleep <- function(work.id){
  library(data.table)
  requireNamespace("pbdMPI") #load to get accurate timing for first job.
  start <- Sys.time()
  my_rank <- pbdMPI::comm.rank()
  Sys.sleep(my_rank/10)
  data.table(rank=my_rank, work.id, start, end=Sys.time(), proc=factor(Sys.getpid()))
}
save(work.id.vec, rank_sleep, file="rank_sleep.RData")
rank_dt <- script_sleep("rank")
ggplot()+
  theme_bw()+
  geom_segment(aes(
    start, rank,
    xend=end, yend=rank),
    data=rank_dt)+
  geom_point(aes(
    start, rank),
    shape=1,
    data=rank_dt)+
  directlabels::geom_dl(aes(
    start, rank,
    label.group=paste(work.id, rank),
    label=work.id),
    data=rank_dt,
    method=my.polygons)+
  scale_y_continuous(breaks=1:workers, limits=c(1, NCPU))+
  scale_x_continuous("Time (seconds)")

sessionInfo()