The objective of this lab is to improve your coding skills by focusing on code debugging, benchmarking and optimization. Below, you will find a number of tasks connected to the topics covered in the Debugging, profiling and optimization lecture. Some tasks extend lectures content and require you to find some more information online. Please, note that while we are providing example solutions to many tasks, these are only examples. If you solve a task in a different way it does not matter your solution is wrong. In fact, it may be better than our solution. If in doubt, ask PI for help. We are here for you!

1 Debugging

1.1 Task: Code Correctness

Which of the following chunks of code are correct and which contain errors? Identify these errors.

1.1.1 Chunk 1

input <- sample(1:1000, size = 1000, replace = T)
currmin <- NULL
for (i in input) {
  if (input > currmin) {
    currmin <- input
    print(paste0("The new minimum is: ", currmin))
  }
}

1.1.2 Chunk 2

input <- sample(1:1000, size = 1000, replac = T)
currmin <- NULL
for (i in input) {
  if (input < currmin) {
    currmin <- input
    print(paste0("The new minimum is: ", currmin))
  }
}

1.1.3 Chunk 3

for (cnt in 1:100) {
  if (cnt > 12) {
    print("12+")
  } else {
    print("Not 12+")
  }
}

1.1.4 Chunk 4

result <- logical(10)
input <- sample(1:10, size = 10, replace = T)
for (i in 0:length(input)) {
  if (input[i] >= 5) {
    result[i] <- TRUE
  }
}

1.2 Task: Debugger.

Play with debugger as described in lecture slides.

1.3 Task: Floating-point Arithmetics.

Can you fix the code below so that it produces more reliable result?
Hint: think in terms of system-specific representation \(\epsilon\).

Put the value of your double \(\epsilon\) into this spreadsheet (Best Coding Practises Lab sheet).

vec <- seq(0.1, 0.9, by=0.1)
vec == 0.7

# One way is to use epsilon
# Check machine's floating point representation
vec <- seq(0.1, 0.9, by=0.1)

# Make a custom function that uses machines' epsilon for comparing
# values
is_equal <- function(x, y) {
  isEqual <- F
  if (abs(x - y) < unlist(.Machine)['double.eps']) {
    isEqual <- T
  }
  isEqual
}

# Some tests
0.7 == 0.6 + 0.1
is_equal(0.7, 0.6 + 0.1)
0.7 == 0.8 - 0.1
is_equal(0.7, 0.8 - 0.1)

# Now you can use the is_equal to fix the code!

2 Profiling

2.1 Task: Filling A Large Matrix.

Create a 10 000 x 10 000 matrix and fill it with random numbers (from 1 to 42), first row by row and later column by column. Use proc.time to see if there is any difference. Is the measurement reliable? Record the values you got in this spreadsheet:

N <- 10e3 * 10e3

# By row
t1 <- proc.time()
M <- matrix(sample(1:42, size = N, replace = T), nrow = sqrt(N), byrow = T)
t2 <- proc.time()
(t2 - t1)

# By column
t1 <- proc.time()
M <- matrix(sample(1:42, size = N, replace = T), nrow = sqrt(N), byrow = F)
t2 <- proc.time()
(t2 - t1)

2.2 Task: Timing Reliability.

In the lecture slides, you have seen how to time sampling from Gaussian distribution:

system.time(rnorm(n = 10e6))

Is such single measurement reliable? Run the code 100 times, plot and record the mean and the variance of the elapsed time. Put these values (elapsed.time mean and variance) into this spreadsheet (Best Coding Practises Lab sheet).

timing <- double(100)
for (i in 1:100) {
  st <- system.time(rnorm(n = 10e6))
  timing[i] <- st[3]
}
boxplot(timing) 
mean(timing)
var(timing)

An alternative approach or, more exactly, an alternative notation that achieves the same as the previous chunk of code but in a more compact way makes use of the replicate, a wrapper function around sapply that simplifies repeated evaluation of expressions. The drawback is you do not get the vector of the actual timing values but the results of calling system.time are already averaged for you. Try to read about the replicate and use it to re-write the code above. Put the elapsed.time into this spreadsheet (Best Coding Practises Lab sheet). How does this value compare to calling system.time within a loop in the previous chunk of code? Are the values similar?

st2 <- system.time(replicate(n = 100, rnorm(n = 10e6)))

2.3 Task: Microbenchmarking.

While system.time might be sufficient most of the time, there is also an excellent package microbenchmark that enables more accurate time profiling, aiming at microsecond resolution that most of modern operating systems offer. Most of the benchmarking the microbenchmark does is implemented in low-overhead C functions and also the package makes sure to: * estimate granularity and resolution of timing for your particular OS, * warm up your processor before measuring, i.e. wake the processor up from any idle state or likewise.

Begin by installing the microbenchmark package.

2.3.1 Checking System Time.

Check the current value of the platform’s timer.

microbenchmark::get_nanotime()

Modify the code below so that it uses the current value of platform’s timer:

timing <- double(100)
for (i in 1:100) {
  st <- system.time(rnorm(n = 10e6))
  timing[i] <- st[3]
}
boxplot(timing)

Put the mean and the variance into this spreadsheet (Best Coding Practises Lab sheet, Microbenchmark – loop)

library(microbenchmark)
timing <- double(100)
for (i in 1:100) {
  nanotime_start <- get_nanotime()
  rnorm(n = 10e6)
  nanotime_stop <- get_nanotime()
  timing[i] <- nanotime_stop - nanotime_start
}
mean(timing)
var(timing)
boxplot(timing)

2.3.2 Microtiming Precision.

There is an experimental function in the microbenchmark package that helps the package estimate granularity and resolution of your particular timing subsystem. According to the documentation, the function measures the overhead of timing a C function call rounds times and returns all non-zero timings observed.
Run the microtiming_precision function and put the mean and the variance of the resulting vector into this spreadsheet (Best Coding Practises Lab sheet, Microbenchmark – precision)

precision <- microbenchmark::microtiming_precision()
mean(precision)
var(precision)

Run the function one time without assigning its value to a variable and consult the documentation. Compare the output of running the function without assigning the value to a variable, the values stored in the variable by the function upon assignment and the value specified in the documentation.

# In version 1.4-4 of the package, all three ways give different results!
microbenchmark::microtiming_precision()
precision <- microbenchmark::microtiming_precision()
?microbenchmark::microtiming_precision

2.3.3 The Microbenchmark Way.

Finally, let’s benchmark our rnorm example using microbenchmark:

microbenchmark the rnorm(n = 10e6) expression,
plot the results using both ggplot2 and a boxplot (read the microbenchmark package documentation),
look at the summary of the benchmark,
how long does it take to dispatch a simple function that does nothing compared to evaluating a constant and adding two integers?

require(microbenchmark)
require(ggplot2)
mb <- microbenchmark(rnorm(n = 10e6))
autoplot(mb)
boxplot(mb)
summary(mb)
f <- function() {}
mb2 <- microbenchmark(f(), pi, 2+2)
summary(mb2)
autoplot(mb2)

2.4 Task: More Advanced Profiling.

Now, we will use a even more sophisticated approach to profiling.

2.4.1 The `Rprof` way.

Write three functions that fill by row a \(N \times N\) matrix \(M\) with randomly generated numbers from a vector given as argument bag, allow for passing random seed value as function argument with the default value of 42. After filling the matrix with values, add to each and every element of \(M\) the number of column the element is in and return such matrix from the function. Functions should:
fill_alloc) – use memory allocation prior to loop in which the matrix is being filled and allocate memory using init value passed as argument and by default set to NULL,
fill_noalloc – not use memory allocation prior to the loop,
fill_noloop should not the loop for filling the matrix in.

** NOTE! ** do not perform addition of column number in the same

Following this and using rnorm(1000, mean = 0, sd = 1):

use Rprof to profile the functinos using the same seed and N=100,
use Rprof to check whether there is a difference between initializing the matrix using NULL and 0 in fill_alloc,
what happens if \(N = 10\) compared to \(N = 20\) to \(N = 100\)?

fill_noloop <- function(N, bag, seed = 42) {
  set.seed(seed)
  values <- sample(bag, size = N^2, replace = T)
  M <- matrix(data = values, nrow = N, byrow = T)
  for (col_num in 1:N) {
    M[, col_num] <- M[, col_num] + col_num
  }
  return(M)
}

fill_noalloc <- function(N, bag, seed = 42) {
  set.seed(seed)
  values <- sample(bag, size = N^2, replace = T)
  M <- NULL
  cnt = 1
  for (row in 1:N) {
    row_tmp <- c()
    for (col in 1:N) {
      row_tmp <- c(row_tmp, values[cnt])
      cnt <- cnt + 1
    }
    M <- rbind(M, row_tmp)
  }
  for (col_num in 1:N) {
    M[, col_num] <- M[, col_num] + col_num
  }
  return(M)
}

fill_alloc <- function(N, bag, seed = 42, init = NA) {
  set.seed(seed)
  values <- sample(bag, size = N^2, replace = T)
  M <- matrix(rep(init, times=N^2), nrow = N, byrow = T)
  cnt = 1
  for (row in 1:N) {
    for (col in 1:N) {
      M[row, col] <- values[cnt]
      cnt <- cnt + 1
    }
  }
  for (col_num in 1:N) {
    M[, col_num] <- M[, col_num] + col_num
  }
  return(M)
}

summary <- summaryRprof('profiler_test_fillers.out', memory='both')
summary$by.self

# answers to the remaining questions are not given

2.5 Optimization

Have a look at our answers.

How can you optimize the fill_alloc even further (call the optimized version fill_alloc_opt)?

fill_alloc_opt <- function(N, bag, seed = 42, init = NA) {
  set.seed(seed)
  values <- sample(bag, size = N^2, replace = T)
  M <- matrix(rep(init, times=N^2), nrow = N, byrow = T)
  cnt = 1
  for (row in 1:N) {
    for (col in 1:N) {
      M[row, col] <- values[cnt] + col
      cnt <- cnt + 1
    }
  }
  return(M)
}

Optimize the fill_noloop to fill_noloops that does not use any loops at all.

fill_noloops <- function(N, bag, seed = 42) {
  values <- sample(bag, size = N^2, replace = T)
  inc <- rep(x = 1:N, times = N)
  M <- matrix(data = values + inc, nrow = N, byrow = T)
  return(M)
}

2.6 Using the `profr` package.

Install and load the profr package.
Use profr to profile fill_noloop, fill_noloops and fill_alloc_opt.

Please use the following fucntion, if you were not able to run profr::ggplot.profr.

ggplot.profr <- function(data, ..., minlabel = 0.1, angle=0) {
  if (!requireNamespace("ggplot2", quietly = TRUE))
    stop("Please install ggplot2 to use this plotting method")
  data$range <- diff(range(data$time))

  # quiet R CMD check note
  start <- NULL
  end <- NULL
  time <- NULL

  ggplot2::ggplot(as.data.frame(data)) +
    ggplot2::geom_rect(
      ggplot2::aes(xmin = start, xmax = end, ymin = level - 0.5, ymax = level + 0.5),
      fill = "grey95", colour = "black", size = 0.5) +
    ggplot2::geom_text(
      ggplot2::aes(start + range / 60, level, label = f),
      data = subset(data, time > max(time) * minlabel),
      size = 4, angle = angle, hjust = 0) +
    ggplot2::scale_y_continuous("time") +
    ggplot2::scale_x_continuous("level")
}

library(profr)
Rprof("profr_noloop.out", interval = 0.01)
fill_noloop(1000, rnorm(1000), seed = 42)
Rprof(NULL)
profile_noloop_df <- parse_rprof('profr_noloop.out')

Rprof("profr_noloops.out", interval = 0.01)
fill_noloops(100, rnorm(1000), seed = 42)
Rprof(NULL)
profile_noloops_df <- parse_rprof('profr_noloops.out')

Rprof("profr_alloc_opt.out", interval = 0.01)
fill_alloc_opt(10, rnorm(1000), seed = 42)
Rprof(NULL)
profile_alloc_opt_df <- parse_rprof('profr_alloc_opt.out')

profr::ggplot.profr(profile_noloop_df)
profr::ggplot.profr(profile_noloops_df)
profr::ggplot.profr(profile_alloc_opt_df)

2.7 Using the `profvis` package.

Install and load the profvis package.
Use profvis to profile fill_noloop, and fill_alloc functions.

3 Optimize Your Code

In this section, we will deal with some selected ways to optimize your code.

3.1 Task: Fix and Optimize This!

Given is a function:

optimize_me <- function(N = 1000, values = c(1:1e4)) {
  N = 10; values = c(1:1e4)
  dat1 <- matrix(size = N^2)
  for (i in 1:N) {
    for (j in 1:N) {
      dat1[i, j] <- sample(values, 1)
    }
  }
  dat0 <- dat1
  dat1[lower.tri(dat1)] <- t(dat1)[lower.tri(dat1)]
  
  dat2 <- NULL 
  for (i in 1:N) {
    i_tmp <- c()
    for (j in 1:N) {
      i_tmp <- c(i_tmp, sample(values, 1))
    }
    dat2 <- rbind(dat2, i_tmp)
  }
  dat2[lower.tri(dat2)] <- t(dat2)[lower.tri(dat2)]
 
  M <- dat2
  for (i in 1:N) {
    for (j in 1:N) {
      M[i, j] <- dat1[i, j] * dat2[i, j]
    }
  }
  for (i in 1:N) {
    for (j in 1:N) {
      M[i, j] <- M[i, j] + values[3]
    }
  }
  N <- M %*% dat0
  result <- apply(N, 2, mean)
  return(result)
}

What does it do, step-by-step?
Profile it.
Is dat1 <- matrix(size = N^2) better than dat1 <- matrix(NA, nrow=N, ncol=N)?
Can you optimize something using BLAS?
Can you optimize by using apply somewhere?
Can you optimize apply further?
What else can you optimize. Do it. Report speed gain and memory gain compared to the original version in this spreadsheet (Best Coding Practises Lab sheet, Optimization gains).

4 Session info

## R version 4.1.2 (2021-11-01)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Monterey 12.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] bsplus_0.1.3      fontawesome_0.2.2 captioner_2.2.3   bookdown_0.26    
## [5] knitr_1.38       
## 
## loaded via a namespace (and not attached):
##  [1] lubridate_1.8.0 digest_0.6.29   R6_2.5.1        jsonlite_1.8.0 
##  [5] magrittr_2.0.3  evaluate_0.15   stringi_1.7.6   rlang_1.0.2    
##  [9] cli_3.2.0       rstudioapi_0.13 jquerylib_0.1.4 bslib_0.3.1    
## [13] generics_0.1.2  rmarkdown_2.14  tools_4.1.2     stringr_1.4.0  
## [17] xfun_0.30       yaml_2.3.5      fastmap_1.1.0   compiler_4.1.2 
## [21] htmltools_0.5.2 sass_0.4.1

Built on: 09-Jun-2022 at 20:14:21.

2022 • SciLifeLab • NBIS • RaukR

Debugging, profiling and optimizing the code

RaukR 2022 • Advanced R for Bioinformatics

Marcin Kierczak

1 Debugging

1.1 Task: Code Correctness

1.1.1 Chunk 1

1.1.2 Chunk 2

1.1.3 Chunk 3

1.1.4 Chunk 4

1.2 Task: Debugger.

1.3 Task: Floating-point Arithmetics.

2 Profiling

2.1 Task: Filling A Large Matrix.

2.2 Task: Timing Reliability.

2.3 Task: Microbenchmarking.

2.3.1 Checking System Time.

2.3.2 Microtiming Precision.

2.3.3 The Microbenchmark Way.

2.4 Task: More Advanced Profiling.

2.4.1 The `Rprof` way.

2.5 Optimization

2.6 Using the `profr` package.

2.7 Using the `profvis` package.

3 Optimize Your Code

3.1 Task: Fix and Optimize This!

4 Session info

Debugging, profiling and optimizing the code

RaukR 2022 • Advanced R for Bioinformatics

Marcin Kierczak

1 Debugging

1.1 Task: Code Correctness

1.1.1 Chunk 1

1.1.2 Chunk 2

1.1.3 Chunk 3

1.1.4 Chunk 4

1.2 Task: Debugger.

1.3 Task: Floating-point Arithmetics.

2 Profiling

2.1 Task: Filling A Large Matrix.

2.2 Task: Timing Reliability.

2.3 Task: Microbenchmarking.

2.3.1 Checking System Time.

2.3.2 Microtiming Precision.

2.3.3 The Microbenchmark Way.

2.4 Task: More Advanced Profiling.

2.4.1 The Rprof way.

2.5 Optimization

2.6 Using the profr package.

2.7 Using the profvis package.

3 Optimize Your Code

3.1 Task: Fix and Optimize This!

4 Session info

2.4.1 The `Rprof` way.

2.6 Using the `profr` package.

2.7 Using the `profvis` package.