Calculating the connectance of a matrix iteratively

Question

I have a performance issue in R. I have a function that iterates over a dataframe for different levels of "site" and "method". The function samples 1:1000 interactions (rows), converts these to matrices and calculates a value of connectance for each matrix. This is repeated n times. The function runs exactly as I want it, returning a dataframe with connectance values for 1 to 1000 interactions when repeated a small number of times. The problem is when I increase the number of repetitions (say to 100) the function runs progressively slower.

df <- read.table(text = "bird_sp plant_sp value site method
 1 species_a plant_a 1 a m
 2 species_a plant_a 1 a m
 3 species_b plant_b 1 a m
 4 species_b plant_b 1 a m
 5 species_c plant_c 1 a m
 6 species_a plant_a 1 b m
 7 species_a plant_a 1 b m
 8 species_b plant_b 1 b m
 9 species_b plant_b 1 b m
 10 species_c plant_c 1 b m
 11 species_a plant_a 1 a f
 12 species_a plant_a 1 a f
 13 species_b plant_b 1 a f
 14 species_b plant_b 1 a f
 15 species_c plant_c 1 a f
 16 species_a plant_a 1 b f
 17 species_a plant_a 1 b f
 18 species_b plant_b 1 b f
 19 species_b plant_b 1 b f
 20 species_c plant_c 1 b f", header = TRUE)
xDegrees <-function(df, size, numRep){
 #Loading required library
 require(bipartite)
 #Creating vector of unique combinations
 df <- within(df, {SiteMethod <- paste(site, method, sep = ":")})
 #Creating empty dataframe
 connectMatrix <- as.data.frame(matrix(rep(0,4), ncol = 4))
 colnames(connectMatrix) <- c("Site","Method","Size","connectance")
 #Beginning of matrix
 k <- 1
 #Beginning subsetting loop
 for(i in 1:length(unique(df$SiteMethod))){
 #subsetting dataset
 dfSub <- subset(df, SiteMethod == unique(df$SiteMethod)[i])
 #Storing site for matrix
 site <- as.character(dfSub[1,]$site)
 #Storing method for matrix
 method <- as.character(dfSub[1,]$method)
 for(l in 1:numRep){
 #Beginning calculation loop
 for(j in 1:length(size)){
 #show progress
 print(paste("S:M", i,j , "completed", sep = " "))
 #The size being calculated
 subSize <- size[j]
 #generate random samples and convert to matrices
 rows <- sample(1:nrow(dfSub), subSize, replace=T)
 intlist <- dfSub[rows,]
 mat <- with(intlist, tapply(value, list(plant_sp, bird_sp), sum))
 mat[is.na(mat)] <- 0
 #network level function to calculate connectance
 con <- networklevel(mat, index = c("connectance"))
 #Stitch matrix together
 connectMatrix[k,] <- c(site, method, subSize, con)
 #Update row
 k <- k + 1
 }
 }
 }
#Return complete matrix
return(connectMatrix)
}
#run the function. 1:1000 interactions, 100 reps
stuff <- xDegrees(df, size = 1:1000, numRep = 100)

Any ideas on how to speed this up?

Answer 1

Consider by to slice your dataframes by the site and method factors and then pass dataframe subsets into an adjusted xDegrees function. Other than the group slice, the other major change is combining for loops with rep and passing result into one sapply call. Below of course is not tested since you do not provide a reproducible example.

# Loading required library
require(bipartite)
xDegrees <-function(dfSub, size, numRep){
 # Retrieving current site
 s <- as.character(dfSub$site[[1]])
 # Retrieving current method
 m <- as.character(dfSub$method[[1]])
 # Build large vector of all random runs
 sample_iter <- rep(seq_along(size), numRep) 
 # Generate random samples and convert to matrices
 con_vec <- sapply(sample_iter, function(i) {
 rows <- sample(1:nrow(dfSub), i, replace=T))
 intlist <- dfSub[rows,]
 mat <- with(intlist, tapply(value, list(plant_sp, bird_sp), sum))
 mat[is.na(mat)] <- 0
 # network level function to calculate connectance
 con <- networklevel(mat, index = c("connectance"))
 })
 # size is a multiple of con_vec (whose length = size * numRep)
 return(data.frame(method = m, site = s, subsize = size, con = con_vec))
}
df_List <- by(df, df[,c("site", "method")], FUN=function(d)
 xDegrees(d, size = 1:1000, numRep = 100))
finald_df <- do.call(rbind, df_List)

Answer 2

Without installing bipartite and getting into the exact details of what the actual calculations are looking to do, there are a few things that jump out at me that might help you.

Growing the Object (connectMatrix)

It looks like you're in Circle 2 - Growing Objects of the R Inferno.
Specifically, the line where you've got:

#Stitch matrix together 
connectMatrix[k,] <- c(site, method, subSize, con)

In that scenario, every time you add another row to the connectMatrix, you're expanding and fragmenting the memory you're using, which as you've observed starts out okay but then starts to really drag. Here's another resource on that same issue, and a much less intricate Stack Overflow question.

nested for loops

You've got three layers of nested for loops

for(i in 1:length(unique(df$SiteMethod))){
 #subsetting dataset
 ...
 ...
 ...
 for(l in 1:numRep){
 #Beginning calculation loop
 for(j in 1:length(size)){
 ...
 ...
 }
 }

The outside for loop jumps out at me as something the plyr package might be good for -- essentially it looks like you're splitting the main dataframe into subsets, then doing things to it (the meat of your second and third for loops). Basically, plyr was designed for splitting-applying-and combining, and there's a really great paper on how to implement it.

You may be able to do two (or maybe even three) layers of splitting-applying-combining.

If the data lend themselves to grouping operations, those tend to help, too, but it doesn't seem like it based on the comments in the code.

Parallel processing

If, after you've optimized for not growing the object and chunked the code for splitting/combining (or another method of getting to the same outcome), you're still not at the level of performance you're looking for, you can journey into parallel application.

Basically, parallel application works well when you're using apply or plyr approaches and you can convert over to mcmapply or one of the other parallel approaches.

For example, we use mcmapply in production at work to take a 30 minute process and trim it down to about 6 minutes. Parallel processing tends to be worth it (in my personal experience) only once you've done everything else you can to get to the best speed you can without parallelizing - it tends to be more work to get it right, and, as a bonus, once you've done the other things it tends to make it easier to convert the process to run on multiple cores.

Finally, Efficient R is a good general resource on optimizing R processes.

• C8H10N4O2: upvoted for good links