Model for ranking data based on rater class membership \\( Z_i \\).
$$Z_i \sim \rr{Multinom}_N \left(\vec{p}_i\right)$$
$$\begin{align} & l(B \given \{\vec{Z}_i, C_{ijk}\} ) = \rr{const} + \\ & \sum_{i=1}^I \left( \sum_{j \in J_i} \sum_{k \in K_{ij}} \left[ \vec{x}_{jk}' B - \rr{log} \sum_{l = 1}^{L_j} e^{ \vec{x}_{jl}' B} \right] \right) \cdot \vec{Z}_i \end{align}$$
What could be parallelized?
Map-Reduce paradigm:
Communication speed becomes a big issue. Could be 15x or more slower than doing locally.
What parts of log-posterior fit this structure?
library(Rmpi)
master.initializeWorkers <- function (mpistt, lgr = NULL)
{
## fix the R exit function
master.fixLast()
## start MPI workers
mpi.spawn.Rslaves(nslaves = mpistt$nworkers)
## all set
invisible(0)
}
master.initializeWorkers()
In the homework example, let
$$\rr{log}\vec{a}_i = \sum_{j \in J_i} \sum_{k \in K_{ij}} \left[ \vec{x}_{jk}' B - \rr{log} \sum_{l = 1}^{L_j} e^{ \vec{x}_{jl}' B} \right]$$
With a prior on \\( Z_i \\):
$$Z_i \sim \rr{Multinom}_N \left(\vec{p}_i\right)$$
$$\begin{align} \rr{log}p(B, Z \given C_{ijk} ) = & \rr{const} + \\ & \sum_{i=1}^I \left( \rr{log} \vec{p}_i + \rr{log}\vec{a}_i\right) \cdot \vec{Z}_i \end{align}$$
Define \\( \theta = \{B, Z \} \\), and then, conveniently, \\( E(\theta) = - \rr{log} p(B, Z \given \rr{data}) \\), and
$$\Pi (X) = \prod_{i=1}^N {e^{- \beta_i E(\theta_i)} \over D(T_i)}$$
with \\( \beta_i = {1 \over T_i }\\) and \\( D(T_i) = \int e^{- \beta_i E(x_i)} dx_i\\).
\\( T_i \\)'s are temperatures for parallel chains.
$$\begin{align} &P_{\rr{replica}} (\theta_i \leftrightarrow \theta_{i+1}) = \\ & = \rr{min} \left( 1, e^{(\beta_i - \beta_{i+1})(E(\theta_i) - E(\theta_{i+1}))}\right)\end{align}$$
Note: replica exchange probability is symmetric.
## broadcast require/library command for all needed packages
mpi.bcast.cmd(cmd = library('MCMCpack'), rank = 0,
nonblock = FALSE)
mpi.bcast.cmd(cmd = library('rstream'), rank = 0,
nonblock = FALSE)
mpi.bcast.cmd(cmd = library('combinat'), rank = 0,
nonblock = FALSE)
## broadcast all needed functions (those whose names start
## with lpl, worker, or mpi)
funcNames <- ls(envir = .GlobalEnv)
lplfuns <- funcNames[grep("lpl|worker|mpi", funcNames)]
for (fun in lplfuns) {
mpi.bcast.NamedObj2slave(obj = eval(as.name(fun)), nam = fun)
}
We are literally populating the little R sessions running on workers.
master.setupRstreams <- function (mpistt, base, lgr = NULL)
{
## NOTE : this creates two global variables in each worker:
## stream -- a variable with all random number generators
## rns -- a list with all generated random number streams
streams <- sapply(1:mpistt$nworkers, function (ii)
{
master.makePackedRstreams(nworkers = mpistt$nworkers,
base = base,
seed = mpistt$seed,
lgr = lgr)
}, simplify = FALSE)
## send random streams to workers
for (ii in 1:mpistt$nworkers) {
mpi.send.Robj(obj = streams[[ii]], dest = ii, tag = 10)
}
## broadcast all workers to receive their streams
## broadcast all workers to generate random number streams
## and store them in the rns object (to be much used later)
mpi.bcast.cmd({
stream <- mpi.recv.Robj(source = mpi.any.source(),
tag = mpi.any.tag())
rns <- worker.generateRstreams(nn = stt$niter,
stream = stream, lgr = lgr)
})
invisible(0)
}
## rstream is faster than other libraries*
library(rstream)
dummy <- new("rstream.mrg32k3a", seed = 1:6)
rstream.sample(dummy, 1)
For details on the rstream library, see this paper.
taken from the reference paper.
Use message passing to communicate between workers.
Worker 1
for (el in elements) {
mpi.isend(el, type = 2, dest = 2, tag = 1)
}
mpi.isend(0, type = 2, dest = 2, tag = 2)
Worker 2
done <- FALSE
while (!done) {
x <- 0
mpi.irecv(x, type = 2, source = mpi.any.source(),
tag = mpi.any.tag())
mpi.wait()
messageinfo <- mpi.get.sourcetag()
workerid <- messageinfo[1]
tag <- messageinfo[2]
if (tag == 1) print(x) else done <- TRUE
}
## worker 1
mpi.send(1:10, type = 2, dest = 2, tag = 0)
## worker 2
x <- numeric(10)
mpi.irecv(x, type = 2, source = 1, tag = 0)
x
mpi.wait()
x
Can send and async receive vectors (unlike hw5 reference code).
mpi.printObj <- function (obj, nworkers, simplify = TRUE)
{
fun <- if (simplify) mpi.iparSapply else mpi.iparLapply
fun(1:nworkers, function (ee)
{
eval(substitute(obj))
})
}
Can use mpi.printObj to print objects from workers (super useful when debugging).
## Run the initial pre-loop commands of the algorithm - require, globals
master.initializeAlgo(mpistt = mpistt,
objs = list(datform = datform, datorig = datorig,
stt = stt),
lgr = lgr)
## Setup the random number streams on the workers -- important to call AFTER
## initializeAlgo because this uses a transferred worker function
master.setupRstreams(mpistt = mpistt, base = rsb, lgr = lgr)
## testing code
#dia <- mpi.printObj(obj = worker(), nworkers = mpistt$nworkers, simplify = FALSE)
#master.printWorkerStatuses(obj = dia, lgr = lgr)
## at this point the workers should be equipped with the needed functions and
## we should be good to go -- time to start the main loop
mpi.bcast.cmd(cmd = { wres <- worker() }, rank = 0, nonblock = FALSE)
## master function
res <- master(mpistt = mpistt, lgr = lgr)
## align the possibly switched temperatures, and for the clean object (draws
## at temperature 1)
resA <- master.alignTemperatures(res = res, lgr = lgr)
clean <- resA[[which(stt$temp == 1)]]
# close slaves and exit.
mpi.close.Rslaves()
master <- function (mpistt, lgr = NULL)
{
nworkers <- mpistt$nworkers
closedworkers <- 0
res <- list()
ii <- 0
while (closedworkers < nworkers && ii < mpistt$maxiter) {
ii <- ii + 1
# Receive a message from a worker
message <- mpi.recv.Robj(mpi.any.source(), mpi.any.tag())
messageinfo <- mpi.get.sourcetag()
workerid <- messageinfo[1]
tag <- messageinfo[2]
if (tag == 1) {
## The worker is giving an update on the progress
say(lgr, 1, message)
} else if (tag == 2) {
## The worker is done, and the message contains the resulting
## object
say(lgr, 1, "Worker ", workerid, " is done")
res[[workerid]] <- message
closedworkers <- closedworkers + 1
}
}
res
}
From worker.irecvObject:
## here define a global variable for receiving buffer
recvObj <<- 2012.2012
mpi.irecv(x = recvObj, type = 2, source = recSrc,
tag = mpi.any.tag())
Need to be constantly monitoring the recvObj buffer variable.
worker.irecvObject <- function (mpiReq, receiveFromSource, valueTag, lgr = NULL)
{
## receive object if there is no unprocessed object currently hanging
## around
if (mpiReq$allowReceiveObj) {
if (!mpiReq$objRecvInProgress) {
mpiReq$obj <- list()
}
## ensure to be receiving objects from a single source
if (!mpiReq$irecvActive) {
recSrc <- receiveFromSource
## here define a global variable for receiving buffer
recvObj <<- 2013.2013
mpi.irecv(x = recvObj, type = 2, source = recSrc,
tag = mpi.any.tag())
mpiReq$irecvActive <- TRUE
}
## check if the object has been received. If yes, populate obj with
## parsed values
if (mpiReq$irecvActive) {
testRecv <- mpi.testany(count = 1)
mpiReq$currentTag <- -1
if (testRecv$flag) {
mpiReq$objRecvInProgress <- TRUE
msginfo <- mpi.get.sourcetag()
src <- msginfo[1]
mpiReq$currentTag <- msginfo[2]
leninfo <- mpi.get.count(type = 2)
if (src != SELF) {
mpiReq$obj$source <- src
for (nam in names(valueTag)) {
if (valueTag[[nam]] == mpiReq$currentTag) {
temp <- recvObj + 0
mpiReq$obj[[nam]] <- worker.assign(temp)
break
}
}
}
mpiReq$irecvActive <- FALSE
}
}
}
mpiReq
}