Parallelization - making it work on a HPC

[JerBub]

Current version uses forking for paralellization. That could be a problem on HPC that have memory restrictions. Better way would be using PSOCK or something simmilar that uses less resources.

[gregorgorjanc]

@JerBub can you link to websites that have led you PSOCK and the error message on the HPC that led you to that website?

[janaobsteter]

So there are two way to parallelise - with or without forking. While forking is quicker to set up and works fine on a laptop, clusters would usually have restrictions regarding forking. If using forking, the entire environment is cloned which also spikes the memory requirement.

Apparently, we only need a minor tweak in the code. We can still use doParallel and foreach packages in R, but just need to add some code

cl <- makeCluster(N, type="PSOCK") # We need to add this
registerDoParallel(cl) # We already do this

foreach(i = 1:100) %dopar% {
    # work
}

stopCluster(cl) # We need to add this

I will make these adjustments and try on a cluster.

[gregorgorjanc]

Hmm, does this mean that we can accommodate both options at the same time in the code? Please test both on HPC and laptop and let's see what you find. THANKS!

[janaobsteter]

I've now implemented this kind of a solution

  if (simParamBee$nThreads > 1) {
     N <- as.numeric(simParamBee$nThreads)
     cl <- makeCluster(N, type="PSOCK")
     registerDoParallel(cl)

     clusterExport(cl, c("SP"))
   } else {
     registerDoParallel(cores = 1)
   }

Why two options? Because setting up this PSOCK (makeCluster()) takes time and is hence significantly slower for a small number of colonies (and unnecessary when using 1 thread). PSOCK also doesn't copy the whole environment so these global variables (SP) need to be exported to each "daughter". I'll test it and see how it works

[janaobsteter]

This was the associated error with the forking:

Error in mcfork() :
  unable to fork, possible reason: Cannot allocate memory
Calls: buildUp ... %dopar% -> <Anonymous> -> mclapply -> lapply -> FUN -> mcfork
Execution halted

[janaobsteter]

To add some explanation here. R uses two types of parallelisation, as described now in the H_Parallelisation vignette. The two types are FORK or PSOCK.

• Forking works by creating subprocesses that share memory and environment with the daughter process, which makes it really quick and efficient - but this is not possible on Windows machine and not allowed on HPCs.
• PSOCK creates each subprocess as a separate R session, and the subprocesses communicate though sockets. This requires copying of the entire environment to each subprocess, making it less efficient and slower. BUT, it's allowed on all machines and OS.

Let's set up an example, let's pretend this is all sequential:

founderGenomes <- quickHaplo(nInd = 20, nChr = 1, segSites = 50)
SP <- SimParamBee$new(founderGenomes)
basePop <- createVirginQueens(founderGenomes)
drones = createDrones(basePop, nInd = 100)
baseColonies <- createMultiColony(basePop)
baseColonies <- cross(baseColonies, drones = drones, crossPlan = "create")
queens <- selectInd(createVirginQueens(baseColonies[[1]], nInd = 100), n = 20, use = "rand")

So, possible options for parallelising this next step are:

1. Don't define any parallel back-end --> this throws a warning saying no parallel backend detected, running sequentially (at least for me on Ubuntu)

SP$nThreads = ncores
reQueen(x = baseColonies, queen = queens)

2. Specify the parallel backend via registerDoParallel - on Linux, I know this does forking --> because this was the default option before but was failing on HPC. On Windows, I don't know what it does

SP$nThreads = ncores
registerDoParallel(cores = SP$nThreads)
reQueen(x = baseColonies, queen = queens)
stopImplicitCluster()

3. Create a FORK cluster - would only work on Linux

SP$nThreads = ncores
cl <- makeCluster(SP$nThreads, type="FORK")
registerDoParallel(cl)
reQueen(x = baseColonies, queen = queens)
stopImplicitCluster()

4. Create a PSOCK cluster - would work everywhere

SP$nThreads = ncores
cl <- makeCluster(SP$nThreads, type="PSOCK")
registerDoParallel(cl)
reQueen(x = baseColonies, queen = queens)
stopImplicitCluster()

When profiling the options above on Linux, we saw that the fastest one was registerDoParallel().

This is confusing, since I thought that on Linux, this does the same as option 3. It's a bit more complex (ChatGPT):

• registerDoParallel() uses an implicit, OS-dependent backend. This calls makePSOCKcluster() internally on Windows, or uses multicore (fork) parallelism on Unix-like systems. It does not return a cluster object — it's handled implicitly.

So given that this acts differently on different OS, this could be a good default! However, when using only one core, this would slow things down on Windows -- ??? @JerBub , could you help with testing this?

However, if this function does forking on Linux, why is it then slower than option 3, if they both do the same thing?
Explanation: On Linux, the registerDoParallel() option is faster because it does not actually create a cluster in the same sense as makeCluster().
Instead, it uses pure multicore forking via parallel::mclapply, which is lighter and faster than any cluster object. It sets the backend to an implicit multicore backend, which uses:

👉 parallel::mcfork()
👉 parallel::mcexit()
👉 parallel::mclapply() under the hood.

This is not equivalent to makeCluster(type="FORK"). It is even lighter.

I suggest:

• test registerDoParallel() on Windows without setting any core or using detectCores(). If it does not cause significant "delays", use that as default. Probably best to initialise this in the SP
• this would not work on HPC though, but we have a vignette for this (we could name it HPC parallelisation) and put something in the manual/notes as well

[gregorgorjanc]

Great job @janaobsteter !!!! Looking forward to hear about Windows tests @JerBub ;)

[janaobsteter]

Final decision: I'm moving the parallelisation from foreach package to future.apply. The code doesn't change much, from

x@colonies = foreach(colony = seq_len(nCol)) %dopar% {
      reQueen(
        x = x[[colony]],
        queen = queen[colony],
        simParamBee = simParamBee
      )
    }

to

  x@colonies = future.apply::future_lapply(X = seq_len(nCol),
                                           FUN = function(colony) {
    reQueen(
      x = x[[colony]],
      queen = queen[colony],
      simParamBee = simParamBee
    )
  })

Since lapply returns a list, there is not need to specify functions on how to combine the output.

The main reason for switching is that the default for future.apply is sequential mode - and for this one, the user doesn't have to specify anything. This also makes it safer to run on all OS, since it doesn't trigger forking, but also safer meaning that the user isn't opening any connections and leaving them open behind. If the user wants to change the behaviour, they just need to do

plan(multicore) - not recommended, this is forking
plan(multisession) - safe, PSOCK
plan(cluster)

For this, I'll first make a PR to merge this option (branch future) into the branch onePSOCK, and then @gregorgorjanc you can continue with reviewing the PR from onePSOCK to devel.

[gregorgorjanc]

Nice find @janaobsteter! I agree with the decision - it looks simpler and "safer" to use. I guess you can change future.apply::future_lapply to future_lapply.