Debugging, Profiling and Optimizing Code

Note

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 TA 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 to get the feel of how does it work.

1.3 Task: Floating-point Arithmetics.

Can you fix the code below so that it produces more reliable result?

Tip

Think in terms of system-specific representation \(\epsilon\).

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

Example Solution

# 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, e.g. proc.time or microbenchmark to see if there is any difference. Is the measurement reliable?

Example Solution

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 the normal 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.

Example Solution

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

Optional

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 the same spreadsheet (Debugging Lab sheet). How does this value compare to calling system.time within a loop in the previous chunk of code? Are the values similar?

Example Solution

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 a 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.

Note

We have noticed that for, e.g. M1 and M2 architectures on MacBooks it does not work well!

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)

Example Solution

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 a 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 the same spreadsheet (Debugging Lab sheet, Microbenchmark – precision)

Example Solution

precision <- microbenchmark::microtiming_precision()
mean(precision)
var(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?

Example Solution

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)

3 Optimize Your Code

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

3.1 Task: Optimize in different ways!

You are given a function that performs a moderately expensive calculation for one input value at a time:

simulate_value <- function(x, n = 5000) {
  set.seed(42)

  values <- rnorm(n, mean = x, sd = 1)
  result <- sum(sin(values)^2 + cos(values)^2)

  return(result)
}

Your goal is to make this calculation more convenient and faster when applied to many values.

Create a vector of inputs 1:100 run the function for every value of x using a for loop and measure how long does it take using microbenchmark (use at least 10 repeats).

Example Solution

x <- 1:1e2

microbenchmark::microbenchmark(
  loop = {
    results <- numeric(length(x))
    for (i in seq_along(x)) {
      results[i] <- simulate_value(x[i])
    }
  },
  times = 10
)

expr	time
loop	51326298
loop	51370285
loop	49558949
loop	50651082
loop	48490322
loop	48618514
loop	47712303
loop	49409405
loop	49061487
loop	52265964

Use base::Vectorize() to create a vectorized version of simulate_value() so that you can write: results <- simulate_vectorized(x).
Benchmark the vectorized version against the baseline.

Example Solution

simulate_vectorized <- Vectorize(simulate_value, "x")
microbenchmark::microbenchmark(
  loop = {
    results <- numeric(length(x))
    for (i in seq_along(x)) {
      results[i] <- simulate_value(x[i])
    }
  },

  vectorized = {
    results <- simulate_vectorized(x)
  },

  times = 10
)

expr	time
loop	44881753
vectorized	43585035
vectorized	47736130
vectorized	48114176
vectorized	51818774
vectorized	47590584
vectorized	44070344
loop	43895940
vectorized	42882271
loop	44815876
loop	43066957
vectorized	43204635
loop	52257234
vectorized	42584193
loop	45195239
loop	42020024
vectorized	40362841
loop	42336178
loop	43725688
loop	43239951

Try to distribute the calculations across multiple CPU cores using the future package and future.apply::future_lapply() function.

Example Solution

future::plan('multisession', workers=4)

microbenchmark::microbenchmark(
  loop = {
    results <- numeric(length(x))
    for (i in seq_along(x)) {
      results[i] <- simulate_value(x[i])
    }
  },

  vectorized = {
    results <- simulate_vectorized(x)
  },

  parallel = {
    results <- unlist(
      future.apply::future_lapply(x, simulate_value, future.seed = TRUE)
    )
  },

  times = 10
)

expr	time
vectorized	40930230
vectorized	40397016
parallel	185050043
loop	43506120
loop	45234701
vectorized	44393927
parallel	96136511
vectorized	41365596
parallel	89382612
vectorized	42226061
parallel	109936441
loop	41676407
loop	41163826
loop	42550986
loop	50108975
vectorized	42293196
loop	42260387
vectorized	40474093
parallel	83873815
loop	42970896
parallel	83593611
loop	41436546
vectorized	40791989
loop	42020536
vectorized	42333338
vectorized	46720738
parallel	119066254
parallel	116482238
parallel	104256127
parallel	97885083

Experiment with different size of x and different plans: ‘sequential’, ‘multicore’ and ‘multisession’.
Why do you thing vectorization had little effect here?
Is the paralellization always strictly better (faster)?

4 Session

Click here

sessionInfo()

R version 4.5.3 (2026-03-11)
Platform: x86_64-conda-linux-gnu
Running under: Ubuntu 26.04 LTS

Matrix products: default
BLAS/LAPACK: /home/roy/miniforge3/envs/r-4.5/lib/libopenblasp-r0.3.33.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Stockholm
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] digest_0.6.39        codetools_0.2-20     microbenchmark_1.5.0
 [4] fastmap_1.2.0        xfun_0.59            knitr_1.51          
 [7] parallel_4.5.3       htmltools_0.5.9      rmarkdown_2.31      
[10] cli_3.6.6            parallelly_1.48.0    future_1.75.0       
[13] compiler_4.5.3       globals_0.19.1       tools_4.5.3         
[16] future.apply_1.20.2  listenv_1.0.0        evaluate_1.0.5      
[19] yaml_2.3.12          otel_0.2.0           rlang_1.3.0         
[22] jsonlite_2.0.0       htmlwidgets_1.6.4