mcmapply demonstration

demo.r

m <- matrix(1:12, nrow=3)

# Double every element
# ===================

# Method 1: nested for loop
for(i in seq_len(nrow(m))) {
  for(j in seq_len(ncol(m))) {
    m2[i, j] <- m[i, j] * 2
  }
}
m2


# Method 2: with mcmapply (same interface as mapply)
library(parallel)
m3 <- mcmapply(m, 
  FUN = function(x) {
    x * 2  
  }
)
# mcmapply returned a vector, so need to set the dims
dim(m3) <- dim(m)
m3


# Pairwise addition of columns, from two matrices
# ===============================================

# Two matrices
m <- matrix(1:12,  nrow=3)
p <- matrix(21:32, nrow=3)


# Preallocate a matrix to store results
res <- matrix(rep(0, ncol(m)^2), ncol = ncol(m))

# Method 1: nested for loop
for(i in seq_len(ncol(m))) {
  for(j in seq_len(i-1)) {
    res[i, j] <- sum(m[,i]) + sum(m[,j]) + sum(p[,i]) + sum(p[,j])
  }
}
res



# Method 2: with mcmapply (same interface as mapply)

# Build the i and j sequence vectors (there may well be a faster way,
# but this is relatively cheap)
# Preallocate vectors
n_pairs <- ncol(m) * (ncol(m)-1) / 2
i_vec <- rep(0, n_pairs)
j_vec <- rep(0, n_pairs)
idx <- 1
for(i in seq_len(ncol(m))) {
  for(j in seq_len(i-1)) {
    i_vec[idx] <- i
    j_vec[idx] <- j
    idx <- idx + 1
  }
}

library(parallel)
computed <- mcmapply(i_vec, j_vec, 
  FUN = function(i, j) {
    # This is where expensive operations should go
    sum(m[,i]) + sum(m[,j]) + sum(p[,i]) + sum(p[,j])
  }
)
# computed was a vector, but we need to put it in the correct-size matrix

# Preallocate a matrix to store results
res <- matrix(rep(0, ncol(m)^2), ncol = ncol(m))

# Copy the elements in
for(n in seq_along(i_vec)) {
  res[i_vec[n], j_vec[n]] <- computed[n]
}
res


# Method 3: mcmapply, with simplified logic

# Use all pairings of i and j
i_vec <- rep(seq_len(ncol(m)), times = ncol(m))
j_vec <- rep(seq_len(ncol(m)), each = ncol(m))
 
library(parallel)
computed <- mcmapply(i_vec, j_vec, 
  FUN = function(i, j) {
    # Just return 0 if not in the lower triangle
    if (i <= j) return(0)
    # This is where expensive operations should go
    sum(m[,i]) + sum(m[,j]) + sum(p[,i]) + sum(p[,j])
  }
)

# computed was a vector, but we need to put it in the correct-size matrix
res <- matrix(computed, ncol = ncol(m))