some functions for running BUGS in parallel

[Kery Marc]

Dear all,

 

With some help of Peter Solymos, Ken Kellner and Mathias Tobler I compiled the attached memo that shows you some functions to run  WinBUGS, OpenBUGS and JAGS in parallel. With jobs that take between hours and days, parallelization can make a real difference.

 

Best regards  -- Marc

_____________________________________________________________
 
Marc Kéry
Tel. ++41 41 462 97 93
marc...@vogelwarte.ch
www.vogelwarte.ch

Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
______________________________________________________________

Willkommen im neuen Besuchszentrum in Sempach! http://www.vogelwarte.ch/de/besuch/
Bienvenu au nouveau centre de visite à Sempach! http://www.vogelwarte.ch/fr/visite/

 

***  Hierarchical modeling books   ***
(4) Kéry & Royle (2017): AHM Volume 2,  Dynamic and advanced models, in prep.
(3) Kéry & Royle (2016): Applied hierarchical modeling in ecology, Academic Press, Volume 1, Prelude and Static Models, see www.mbr-pwrc.usgs.gov/pubanalysis/keryroylebook/
(2) Kéry & Schaub (2012): Bayesian Population Analysis using WinBUGS, Academic Press; see www.vogelwarte.ch/bpa
(1) Kéry (2010): Introduction to WinBUGS for Ecologists, Academic Press; see www.mbr-pwrc.usgs.gov/pubanalysis/kerybook 

***   Hierarchical modeling workshops ***
www.phidot.org/forum/viewforum.php?f=8

***   Hierarchical modeling Google Group mailing lists   ***
(1) unmarked: for questions specific to the R package unmarked
(2) hmecology: for material covered in Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2015); see groups.google.com/forum/?hl=en#!forum/hmecology

 

[Chris SUtherland]

Marc, Peter, Ken and Matthias,

 

Great resource!

Thanks!

--
*** Three hierarchical modeling email lists ***

(1) unmarked: for questions specific to the R package unmarked

(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2015)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.
To post to this group, send email to hmec...@googlegroups.com.
Visit this group at https://groups.google.com/group/hmecology.
To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/F801018864E940468323F405DF95D2DD17587999%40MailNT.vogelwarte.ch.
For more options, visit https://groups.google.com/d/optout.

[Guillaume Souchay]

Thanks guys for this !

I was aware of the parralel option in jagsUI but it's cool for BUGS model.

Thanks !!

Guillaume

[Simone Tenan]

Hi all,

I resume this thread to ask whether someone has experience of fitting a model to different datasets in parallel. I used the following example code with the mclapply() function of the parallel package to a list of datasets in a server with Ubuntu 20.04.1 LTS (GNU/Linux 5.4.0-64-generic x86_64).

Unfortunately, sometime some of the elements of the list produced by mclapply() ("outLs" in the code below) are empty.

E.g.

str(outLs[[11]])
 NULL

I wonder why, given that when I try fitting the model to the specific data set everything works.

Every alternative approach that you can suggest is welcome!

Thanks

Simone

# R packages
library(jagsUI)
library(parallel)

paral_fun <- function(x){
    #### Read in data #########################################
    CH <- as.matrix(read.table(file=x,sep=",",header=T))

    # Load models ##########################################
    source("modelsDD.R")

    # JAGS stuff ###########################################
    # Bundle data
    bugs.data <- list(...)

    # Initial values
    inits <- function(){list(...
                             )}

    # Parameters to be monitored
    parameters <- c("..."
                   )

    # MCMC settings
    n.adapt <-    50000
    n.burnin <-   100000
    n.iter <-     n.burnin+30000
    thin <-       3
    chains <-     3

    # run Model
    out <- jags(data = bugs.data,
                parameters.to.save = parameters,
                model = "model.txt",
                inits = inits,
                n.adapt = n.adapt,
                n.burnin = n.burnin,
                n.iter = n.iter,
                n.thin = thin,
                n.chains = chains,
                parallel=F)

    return(list(bugs.data=bugs.data,out=out))
    ###################################################
} # end of function

########################################################
# Execute the function across multiple nodes

outLs <- mclapply(FileList, paral_fun, mc.cores=30)

save(outLs, file="outLs.RData", ascii=TRUE)

[Ken Kellner]

Hi Simone,

My understanding is that mclapply() starts by assigning NULL to all elements of the output list. As it is running the job, if any of the parallel worker processes fail for some reason, nothing will be returned and the corresponding list element will remain NULL (apparently with no warning). Parallel workers could fail for a variety of reasons that depend on your exact hardware setup. Probably the most likely is that they run out of memory. The fact that you are seeing this happen randomly is good evidence to me that it is a memory issue. I'd consider trying fewer cores so that each core has more memory available.

Ken

[Simone Tenan]

Hi Ken,

Thanks for your reply. I monitored the amount of total memory used by all processes and it was far below the limit, but I don't know whether there is a memory limit each forked process can use irrespective of the total RAM available. I noticed that it works running a few processes (<10) in parallel and specifying a number of cores equal to the number of processes in mclapply(). In this case, however, I lose the advantage of running many processes in parallel. Then, in my first run I specified 30 cores for 74 processes; however, using 30 cores for 30 processes ends up with empty objects too. Do you know alternative approaches?  

Thanks!!

Simone

[Peter Solymos]

Simone,

Is there a particular reason why you are using mclapply()? Forking is known to have mysterious problems that are hard to troubleshoot. See e.g. the Warning on the ?mclapply home page:

"It is strongly discouraged to use these functions in GUI or embedded environments, because it leads to several processes sharing the same GUI which will likely cause chaos (and possibly crashes). Child processes should never use on-screen graphics devices."

You can use socket clusters and parLapply instead. Your code also does not seem to set random numbers for the parallel processes, which can be potentially dangerous (i.e. identical samples). Neither mclapply nor parLapply takes care of the parallel numbers for JAGS, and JAGS's initialization is independent of how seed is set in the R process. It is a pet peeve of mine, I know, and others don't seem to care that much, but a great deal of thought was put into this by Martyn Plummer, and it is good practice to follow defaults put in place by a software's creator.

Some refactoring was done in rjags to expose the initialization part, which is used in dclone::parallel.inits that is used in a couple of downstream functions for convenience. So it is safe to do something like this:

library(dclone)
ncl <- 3
cl <- makePSOCKcluster(ncl)
clusterEvalQ(cl, library(dclone))
parLoadModule(cl, "glm")
pm <- jags.parfit(cl, data, "beta", model, n.chains=ncl)
stopCluster(cl)

Alternatively, if you want to monitor and possibly sample further:

cl <- makePSOCKcluster(ncl)
parJagsModel(cl, name="res", file=model, data=data,
    n.chains = ncl, n.adapt=1000)
parUpdate(cl, "res", n.iter=1000) # burnin
m <- parCodaSamples(cl, "res", jpara, n.iter=2000)
# you can keep sampling here until you shut down the cluster
stopCluster(cl)

Cheers,

Peter

--
*** Three hierarchical modeling email lists ***

(1) unmarked: for questions specific to the R package unmarked

(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture

(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2020)

---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/359f7b39-9a01-4466-ac59-4c955f411c39n%40googlegroups.com.

[Simone Tenan]

Hi all,

Peter Solymos kindly suggested a way to tackle the problem I had of independent parallel jobs (thanks Peter!). I paste here below his reply.

Regards,

Simone

----

If the parallel jobs are independent you should be able to specify the data as a list, and then passing the list elements as data to the parLapply function, RNG should be fine because you can initiate the chains from the same call. If the data sets are really different in sample size the DAG graph size will be different and computing time will vary so you can use parLapplyLB which uses load balancing and fills in idle nodes.

library(dclone)

## your models
mod_fun <- function() {
    for (i in 1:N) {
        Y[i] ~ dnorm(mu[i], tau)
        mu[i] <- alpha + beta * (x[i] - x.bar)
    }
    x.bar <- mean(x[])
    alpha ~ dnorm(0.0, 1.0E-4)
    beta ~ dnorm(0.0, 1.0E-4)
    sigma <- 1.0/sqrt(tau)
    tau ~ dgamma(1.0E-3, 1.0E-3)
}

## data generation
set.seed(1234)
N <- 100
alpha <- 1
beta <- -1
sigma <- 0.5
x <- runif(N)
linpred <- crossprod(t(model.matrix(~x)), c(alpha, beta))

Y1 <- rnorm(N, mean = linpred, sd = sigma)
Y2 <- rnorm(N*2, mean = -linpred, sd = 2*sigma)
## list of data for the model
Data <- list(
  data1=list(N = N, Y = Y1, x = x),
  data2=list(N = N*2, Y = Y2, x = c(x, x))
)
## what to monitor
par <- c("alpha", "beta", "sigma")

test <- jags.fit(Data[[1]], par, mod_fun)
test <- jags.fit(Data[[2]], par, mod_fun)

## set up cluster
cl <- makeCluster(length(Data))

clusterEvalQ(cl, library(dclone))
clusterExport(cl, c("par", "mod_fun"))

Mods <- parLapply(cl, Data, function(d) {
  jags.fit(d, par, mod_fun, n.chains = 3, n.update=100, n.iter=100)
})

stopCluster(cl)

## you can also use a function if you'd like to save the output
## instead of passing it back to the main process:
## dname is the data set name from the list
## this is safer because you won't loose results that are already done

par_fun <- function(dname) {
  out <- try({
    m <- jags.fit(Data[[dname]], par, mod_fun,
                  n.chains = 3, n.update=100, n.iter=100)
    save(m, file=paste0(dname, ".RData"))
  })
  !inherits(out, "try-error")
}


cl <- makeCluster(length(Data))
clusterEvalQ(cl, library(dclone))
clusterExport(cl, c("par", "mod_fun", "Data"))
## this is where files will be saved
clusterEvalQ(cl, getwd())

OK <- parLapply(cl, names(Data), par_fun)

stopCluster(cl)