Debugging, Profiling and Optimization

.title[
# Debugging, Profiling and Optimization
]
.subtitle[
## RaukR 2022 • Advanced R for Bioinformatics
]
.author[
### <b>Ashfaq Ali (Slides Credit Marcin Kierczak)</b>
]
.institute[
### NBIS, SciLifeLab
]

---

---

## Contents

* [Overview](#overview)
* [Types of bugs](#types-of-bugs)
* [Handling Errors](#error-handling-1)
* [Debugging -- Errors and Warnings](#errors-and-warnings)
* [Debugging -- What are my Options?](#debugging-options)
* [Profiling](#profiling_proc_time)
* [Advanced profiling - **`Rprof`**](#Rprof)
* [Profiling **`profr`** package](#profr_package)
* [Profiling with **`profvis`** package](#profviz_package)
* [Optimizing your code](#optimizing_code1)
* [Copy-on-modify and memory allocation](#copy-on-modify)
* [Allocating memory](#memory-benchmark)
* [GPU](#gpu-1)
* [Parallelization using package **`parallel`**](#prallelization)

---
name: overview
## Topics of This Presentation
<br><br>
Code:
<br>

* **Debugging** -- my code does not run.

* **Profiling** -- now it does run but... out of memory!

* **Optimization** -- making things better.

---

### There are different types of bugs we can introduce:
* Syntax -- `prin(var1), mean(sum(seq((x + 2) * (y - 9 * b)))`

* Arithmetic -- `x/0` (not in R, though!) `Inf/Inf`

* Type -- `mean('a')`

* Logic -- everything works and produces seemingly valid output that is WRONG!

### To avoid bugs:
* Encapsulate your code in smaller units (functions), you can test.

* Use classes and type checking.

* Test at the boundaries, e.g. loops at min and max value.

* Feed your functions with test data that should result with a known output.

* Use *antibugging*: `stopifnot(y <= 75)`

---
name: arithmetic-bugs
## Arithmetic bugs

```r
(vec <- seq(0.1, 0.9, by=0.1))
vec == 0.7 
vec == 0.5
(0.5 + 0.1) - 0.6
(0.7 + 0.1) - 0.8 # round((0.7 + 0.1) , digits = 2) - 0.8
```

```
## [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
## [1] 0
## [1] -1.110223e-16
```
Beware of floating point arithmetics!
.small[

```r
head(unlist(.Machine))
head(unlist(.Platform))
```

```
##     double.eps double.neg.eps    double.xmin 
##   2.220446e-16   1.110223e-16  2.225074e-308 
##    double.xmax    double.base  double.digits 
##  1.797693e+308   2.000000e+00   5.300000e+01 
##                    OS.type                   file.sep 
##                     "unix"                        "/" 
##                 dynlib.ext                        GUI 
##                      ".so"                  "RStudio" 
##                     endian                    pkgType 
##                   "little" "mac.binary.big-sur-arm64"
```
]

---
name: error-handling-1

## Handling Errors

```r
input <- c(1, 10, -7, -2/5, 0, 'char', 100, pi)
for (val in input) {
  (paste0('Log of ', val, 'is ', log10(val)))
}
```

```
## Error in log10(val): non-numeric argument to mathematical function
```

One option is to use the `try` block:

```r
for (val in input) {
  val <- as.numeric(val)
  try(print(paste0('Log of ', val, ' is ', log10(val))))
}
```

```
## [1] "Log of 1 is 0"
## [1] "Log of 10 is 1"
## [1] "Log of -7 is NaN"
## [1] "Log of -0.4 is NaN"
## [1] "Log of 0 is -Inf"
## [1] "Log of NA is NA"
## [1] "Log of 100 is 2"
## [1] "Log of 3.14159265358979 is 0.497149872694133"
```

---
name: error-handling-2

## Handling Errors with `tryCatch`

```r
for (val in input) {
  val <- as.numeric(val)
  result <- tryCatch(log10(val), 
              warning = function(w) { print('Negative argument supplied. Negating.'); log10(-val) }, 
              error = function(e) { print('Not a number!'); NaN })
  print(paste0('Log of ', val, ' is ', result))
}
```

```
## [1] "Log of 1 is 0"
## [1] "Log of 10 is 1"
## [1] "Negative argument supplied. Negating."
## [1] "Log of -7 is 0.845098040014257"
## [1] "Negative argument supplied. Negating."
## [1] "Log of -0.4 is -0.397940008672038"
## [1] "Log of 0 is -Inf"
## [1] "Log of NA is NA"
## [1] "Log of 100 is 2"
## [1] "Log of 3.14159265358979 is 0.497149872694133"
```

---
name: errors-and-warnings

## Debugging -- Errors and Warnings
* An error in your code will result in a call to the `stop()` function that:

+ breaks the execution of the program (loop, if-statemetnt, etc.),
 
  + performs the action defined by the global parameter `error`.
  
* A warning just prints out the warning message (or reports it in another way).

* Global parameter `error` defines what R should do when an error occurs.

```r
options(error = )
```

* You can use `simpleError()` and `simpleWarning()` to generate errors and warnings in your code:

```r
f <- function(x) {
  if (x < 0) {
    x <- abs(x)
    w <- simpleWarning("Value less than 0. Taking abs(x)")
    w
  }
}
```

---
name: debugging-options

## Debugging -- What are my Options?

* Old-school debugging: a lot of `print` statements
  + print values of your variables at some checkpoints,
  + sometimes fine but often laborious,
  + need to remove/comment out manually after debugging.
  
* Dumping frames
  + on error, R state will be saved to a file,
  + file can be read into debugger,
  + values of all variables can be checked,
  + can debug on another machine, e.g. send dump to your colleague!

* Traceback
  + a list of the recent function calls with values of their params,

* Step-by-step debugging
  + execute code line by line within the debugger
  
---
name: debugging_dump_frames
## Option 1: Dumping Frames

```r
options(error = quote(dump.frames("testdump", TRUE)))

f <- function(x) {
    sin(x)
}
f('test')
```

```
## Error in sin(x): non-numeric argument to mathematical function
```

```r
options(error = NULL)
load("testdump.rda")
# debugger(testdump)
```

.smaller[<tt>Message:  Error in sin(x) : non-numeric argument to mathematical function <br>
Available environments had calls: <br>
1: f("test") <br>
 <br>
Enter an environment number, or 0 to exit   <br>
Selection: 1 <br>
Browsing in the environment with call: <br>
   f("test") <br>
Called from: debugger.look(ind) <br>
Browse[1]> x <br>
[1] "test" <br>
Browse[1]>  <br>
[1] "test" <br>
Browse[1]> 
</tt>
]
Last empty line brings you back to the environments menu.

---
name: debugging_traceback
## Option 2: Traceback

```r
f <- function(x) {
  log10(x)  
}

g <- function(x) {
  f(x)
}
g('test')
```

```
## Error in log10(x): non-numeric argument to mathematical function
```
<tt>
> traceback()<br>
2: f(x) at #2<br>
1: g("test")<br>
</tt>

`traceback()` shows what were the function calls and what parameters were passed to them when the error occured.

---
name: step_by_step-1
## Option 3: Debug step-by-step

Let us define a new function `h(x, y)`:

```r
h <- function(x, y) { 
  f(x) 
  f(y) 
  }
```
Now, we can use `debug()` to debug the function in a step-by-step manner:

```r
debug(h)
h('text', 7)
undebug(h)
```

---
name: step_by_step-2
## Option 3: Debug step-by-step cted.

`n` -- execute next line, `c` -- execute whole function, `q` -- quit debugger mode.

<tt>
> debug(h)<br>
> h('text', 7)<br>
debugging in: h("text", 7)<br>
debug at #1: {<br>
    f(x)<br>
    f(y)<br>
}<br>
</tt>
--
<tt>Browse[2]> x<br>
[1] "text"<br></tt>
--
<tt>Browse[2]> y<br>
[1] 7<br></tt>
--
<tt>Browse[2]> n<br>
debug at #2: f(x)<br></tt>
--
<tt>Browse[2]> x<br>
[1] "text"<br></tt>
--
<tt>Browse[2]> n<br>
Error in log10(x) : non-numeric argument to mathematical function<br>
</tt>

---

## Profiling -- `proc.time()`
Profiling is the process of identifying memory and time bottlenecks in your code.

```r
proc.time()
```

```
##     user   system  elapsed 
##  379.248   53.761 1980.224
```
* `user time` -- CPU time charged for the execution of user instructions of the calling process, 
* `system time` -- CPU time charged for execution by the system on behalf of the calling process,
* `elapsed time` -- total CPU time elapsed for the currently running R process.

```r
pt1 <- proc.time()
tmp <- runif(n =  10e5)
pt2 <- proc.time()
pt2 - pt1
```

```
##    user  system elapsed 
##   0.006   0.000   0.006
```
---
name: profiling_system_time

## Profiling -- `system.time()`

```r
system.time(runif(n = 10e6))
system.time(rnorm(n = 10e6))
```

```
##    user  system elapsed 
##   0.230   0.008   0.239 
##    user  system elapsed 
##   0.254   0.006   0.261
```

An alternative approach is to use `tic` and `toc` statements from the `tictoc` package.

```r
library(tictoc)
tic()
tmp1 <-runif(n = 10e6)
toc()
```

```
## 0.056 sec elapsed
```

---
name: profiling_system_time

## Profiling in Action

Let's see profiling in action! We will define four functions that fill a large vector in two different ways:
--

```r
fun_fill_loop1 <- function(n = 10e6, f) {
  result <- NULL
  for (i in 1:n) {
    result <- c(result, eval(call(f, 1)))
  }
  return(result)
}
```
--

```r
fun_fill_loop2 <- function(n = 10e6, f) {
  result <- vector(length = n)
  for (i in 1:n) {
    result[i] <- eval(call(f, 1))
  }
  return(result)
}
```
 
---
name: profiling_in_action

## Profiling in Action cted.

It is maybe better to use...
--
vectorization!
--

```r
fun_fill_vec1 <- function(n = 10e6, f) {
  result <- NULL
  result <- eval(call(f, n))
  return(result)
}
```
--

```r
fun_fill_vec2 <- function(n = 10e6, f) {
  result <- vector(length = n)
  result <- eval(call(f, n))
  return(result)
}
```

---
name: compare_loop_vec_sys_time

## Profiling our functions

```r
system.time(fun_fill_loop1(n = 10e4, "runif")) # Loop 1
system.time(fun_fill_loop2(n = 10e4, "runif")) # Loop 2
system.time(fun_fill_vec1(n = 10e4, "runif"))  # Vectorized 1
system.time(fun_fill_vec2(n = 10e4, "runif"))  # Vectorized 2
```

```
##    user  system elapsed 
##   7.750   3.240  11.067 
##    user  system elapsed 
##   0.168   0.015   0.184 
##    user  system elapsed 
##   0.001   0.001   0.000 
##    user  system elapsed 
##   0.001   0.000   0.001
```

The `system.time()` function is not the most accurate though. During the lab, we will experiment with package `microbenchmark`.

---
name: Rprof

## More advanced profiling
We can also do a bit more advanced profiling, including the memory profiling, using, e.g. `Rprof()` function.

And let us summarise:
.small[

```r
summary <- summaryRprof('profiler_test.out', memory='both')
datatable(summary$by.self, options=list(pageLength = 10, searching = F, info = F))
#knitr::kable(summary$by.self)
```

<div id="htmlwidget-117a767cea22b17234b4" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-117a767cea22b17234b4">{"x":{"filter":"none","vertical":false,"data":[["\"print.default\"","\"runif\"","\"eval\"","\"fun_fill_loop2\"","\"paste0\"","\"parent.frame\"","\"print\"","\"w$get_new\"","\"is.list\"","\"is.pairlist\""],[36.47,0.67,0.21,0.06,0.05,0.02,0.01,0.01,0.01,0.01],[97.2,1.79,0.56,0.16,0.13,0.05,0.03,0.03,0.03,0.03],[36.47,0.67,37.46,0.98,0.05,0.02,36.48,0.06,0.01,0.01],[97.2,1.79,99.84,2.61,0.13,0.05,97.23,0.16,0.03,0.03],[[36490.3],[1305.7],[38303.2],[1783.6],[8.6],[41.9],[36519.6],[17.1],[31.4],[25.6]]],"container":"<table class=\"display\">\n  <thead>\n    <tr>\n      <th> <\/th>\n      <th>self.time<\/th>\n      <th>self.pct<\/th>\n      <th>total.time<\/th>\n      <th>total.pct<\/th>\n      <th>mem.total<\/th>\n    <\/tr>\n  <\/thead>\n<\/table>","options":{"pageLength":10,"searching":false,"info":false,"columnDefs":[{"className":"dt-right","targets":[1,2,3,4,5]},{"orderable":false,"targets":0}],"order":[],"autoWidth":false,"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
]

---
name: profr_package

## Profiling -- `profr` package
There are also packages available that enable even more advanced profiling:

```r
library(profr)
Rprof("profiler_test2.out", interval = 0.01)
tmp <- table(sort(rnorm(1e5)))
Rprof(NULL)
profile_df <- parse_rprof('profiler_test2.out')
```
This returns a table that can be visualised:

<div id="htmlwidget-f5523aec530de00091dd" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-f5523aec530de00091dd">{"x":{"filter":"none","vertical":false,"data":[["1","2","3","4","5","6","7","8","9","10","11","12"],[1,2,2,2,3,3,3,3,4,4,4,4],[1,1,2,3,1,2,3,4,1,2,3,4],[1,1,1,1,1,1,1,1,1,1,1,1],["table","sort","factor","as.integer","rnorm","sort.default","unique","as.character","sort.int","unique.default","as.character","unique.default"],[0,0,0.04,0.6,0,0.02,0.04,0.4,0.02,0.04,0.06,0.38],[0.62,0.04,0.6,0.62,0.02,0.04,0.4,0.58,0.04,0.06,0.38,0.4],[1,1,1,1,1,1,1,1,1,1,1,1],[false,false,false,true,true,false,false,true,true,true,true,true],[0.62,0.04,0.56,0.02,0.02,0.02,0.36,0.18,0.02,0.02,0.32,0.02],["base","base","base","base","stats","base","base","base","base","base","base","base"]],"container":"<table class=\"display\">\n  <thead>\n    <tr>\n      <th> <\/th>\n      <th>level<\/th>\n      <th>g_id<\/th>\n      <th>t_id<\/th>\n      <th>f<\/th>\n      <th>start<\/th>\n      <th>end<\/th>\n      <th>n<\/th>\n      <th>leaf<\/th>\n      <th>time<\/th>\n      <th>source<\/th>\n    <\/tr>\n  <\/thead>\n<\/table>","options":{"pageLength":5,"searching":false,"info":false,"columnDefs":[{"className":"dt-right","targets":[1,2,3,5,6,7,9]},{"orderable":false,"targets":0}],"order":[],"autoWidth":false,"orderClasses":false,"lengthMenu":[5,10,25,50,100]}},"evals":[],"jsHooks":[]}</script>

---
name: profr_package_cted

## Profiling -- `profr` package cted.

We can also plot the results using -- `proftools` package-

```r
library("proftools") # depends on "graph" and "Rgraphviz" packages
profile_df2 <- readProfileData("profiler_test2.out")
plotProfileCallGraph(profile_df2, style = google.style, score = "total")
```

---
name: profviz_package

## Profiling with `profvis`

Yet another nice way to profile your code is by using Hadley Wickham's `profvis` package:

```r
library(profvis)
profvis({fun_fill_loop2(1e4, 'runif')
  fun_fill_vec2(1e4, 'runif')
  })
```

---
name: profvis_run

## Profiling with `profvis` cted.
<div id="htmlwidget-6739b97fb34a83c98047" style="width:100%;height:600px;" class="profvis html-widget"></div>
<script type="application/json" data-for="htmlwidget-6739b97fb34a83c98047">{"x":{"message":{"prof":{"time":[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,3,3,3,3,4,4,4,4],"depth":[18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,2,1,4,3,2,1,4,3,2,1],"label":["findCenvVar","getInlineInfo","tryInline","cmpCall","cmp","cmpSymbolAssign","h","tryInline","cmpCall","cmp","h","tryInline","cmpCall","cmp","genCode","cmpfun","compiler:::tryCmpfun","rmarkdown::render","is.pairlist","rmarkdown::render","eval","eval","fun_fill_loop2","rmarkdown::render","eval","eval","fun_fill_loop2","rmarkdown::render"],"filenum":[null,null,null,null,null,null,null,null,null,null,null,null,null,null,null,null,1,null,null,null,null,null,1,null,null,null,1,null],"linenum":[null,null,null,null,null,null,null,null,null,null,null,null,null,null,null,null,2,null,null,null,null,null,2,null,null,null,2,null],"memalloc":[89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,89.4757232666016,93.9495315551758,93.9495315551758,101.611190795898,101.611190795898,101.611190795898,101.611190795898,108.412048339844,108.412048339844,108.412048339844,108.412048339844],"meminc":[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4.47380828857422,0,7.66165924072266,0,0,0,6.80085754394531,0,0,0],"filename":[null,null,null,null,null,null,null,null,null,null,null,null,null,null,null,null,"<expr>",null,null,null,null,null,"<expr>",null,null,null,"<expr>",null]},"interval":10,"files":[{"filename":"<expr>","content":"library(profvis)\nprofvis({fun_fill_loop2(1e4, 'runif')\n  c()\n  })","normpath":"<expr>"}],"prof_output":"/Users/cob-aaf/Documents/GitHub/RaukR-2022/debugging/.Rproj.user/77B5F86E/profiles-cache/file16b2e46f28027.Rprof","highlight":{"output":["^output\\$"],"gc":["^<GC>$"],"stacktrace":["^\\.\\.stacktraceo(n|ff)\\.\\.$"]},"split":"h"}},"evals":[],"jsHooks":[]}</script>

---
name: optimizing_code1

## Optimizing your code

*We should forget about small efficiencies, say about 97% of the   time: premature optimization is the root of all evil. Yet we should not pass up our opportunities in that critical 3%. A good programmer will not be delulled into complacency by such reasoning, he will be wise to look carefully at the critical code; but only after that code has been identified.*
<br><br>
*-- Donald Knuth*

<div class="pull-left">
<img src="./assets/xkcd_automation.png" style="height:300px;">
<br>
.vsmall[source: http://www.xkcd/com/1319]
</div>

<div class="pull-right">
<img src="./assets/xkcd_is_it_worth_the_time_2x.png" style="height:300px;">
<br>
.vsmall[source: http://www.xkcd/com/1205]
</div>

---
name: optimization_types

## Ways to optimize the code

* write it in a more efficient way, e.g. use vectorization or `*apply` family instead of loops etc.,
* allocating memory to avoid copy-on-modify,
* use package `BLAS` for linear algebra,
* use `bigmemory` package,
* GPU computations,
* multicore support, e.g. `multicore`, `snow`
* use `data.table` or `tibble` instead of `data.frame`

---
name: copy-on-modify

## Copy-on-modify

```r
library(pryr)
order <- 1024
matrix_A <- matrix(rnorm(order^2), nrow = order)
matrix_B <- matrix_A
```
--
Check where the objects are in the memory:
--

```r
address(matrix_A)
address(matrix_B)
```

```
## [1] "0x2873f8000"
## [1] "0x2873f8000"
```
--
What happens if we modify a value in one of the matrices?
--

```r
matrix_B[1,1] <- 1
address(matrix_A)
address(matrix_B)
```

```
## [1] "0x2873f8000"
## [1] "0x287bfc000"
```

---
name: avoid-copy-on-modify-1

## Avoid copying by allocating memory

### No memory allocation

```r
f1 <- function(to = 3, silent=F) {
  tmp <- c()
  for (i in 1:to) {
    a1 <- address(tmp)
    tmp <- c(tmp, i)
    a2 <- address(tmp)
    if(!silent) { print(paste0(a1, " --> ", a2)) } 
  }
}
f1()
```

```
## [1] "0x11d01e8e0 --> 0x12849e648"
## [1] "0x12849e648 --> 0x12849e9c8"
## [1] "0x12849e9c8 --> 0x17253ee08"
```

---
name: avoid-copy-on-modify-1
## Avoid copying by allocating memory cted.

### With allocation

```r
f2 <- function(to = 3, silent = FALSE) {
  tmp <- vector(length = to, mode='numeric')
  for (i in 1:to) {
    a1 <- address(tmp)
    tmp[i] <- i
    a2 <- address(tmp)
    if(!silent) { print(paste0(a1, " --> ", a2)) }
  }
}
f2()
```

```
## [1] "0x173870898 --> 0x173870898"
## [1] "0x173870898 --> 0x173870898"
## [1] "0x173870898 --> 0x173870898"
```

---
name: memory-benchmark

## Allocating memory -- benchmark.

```r
library(microbenchmark)
benchmrk <- microbenchmark(f1(to = 1e3, silent = T), 
                           f2(to = 1e3, silent = T), 
                           times = 100L)
autoplot(benchmrk)
```

---
name: gpu-1
## GPU

```r
library(gpuR)
library(microbenchmark)
A = matrix(rnorm(1000^2), nrow=1000) # stored: RAM, computed: CPU
B = matrix(rnorm(1000^2), nrow=1000) 
gpuA = gpuMatrix(A, type = "float") # stored: RAM, computed: GPU
gpuB = gpuMatrix(B, type = "float")
vclA = vclMatrix(A, type = "float") # stored: GPU, computed: GPU
vclB = vclMatrix(B, type = "float")
bch <- microbenchmark(
  cpuC = A %*% B,
  gpuC = gpuA %*% gpuB,
  vclC = vclA %*% vclB, 
  times = 10L) 
```
.small[
More on [Charles Determan's Blog](https://www.r-bloggers.com/r-gpu-programming-for-all-with-gpur/).
]

---
name:gpu-2

## GPU cted.

```r
autoplot(bch)
```

---
name: prallelization

## Parallelization using package `parallel`
Easiest to paralllelize is `lapply`:

```r
result <- lapply(1:2, function(x) { c(x, x^2, x^3) } )
result
```

```
## [[1]]
## [1] 1 1 1
## 
## [[2]]
## [1] 2 4 8
```

```r
library(parallel)
num_cores <- detectCores() - 1
cl <- makeCluster(num_cores) # Init cluster
parLapply(cl, 1:2, function(x) { c(x, x^2, x^3)} )
stopCluster(cl)
```

```
## [[1]]
## [1] 1 1 1
## 
## [[2]]
## [1] 2 4 8
```

---
name: end-slide
class: end-slide, middle
count: false

# Thank you. Questions?

<p>R version 4.1.2 (2021-11-01)<br><p>Platform: aarch64-apple-darwin20 (64-bit)</p><p>OS: macOS Monterey 12.1</p><br>

Built on : <i class='fa fa-calendar' aria-hidden='true'></i> 09-Jun-2022 at <i class='fa fa-clock-o' aria-hidden='true'></i> 20:08:55

<b>2022</b> • [SciLifeLab](https://www.scilifelab.se/) • [NBIS](https://nbis.se/) • [RaukR](https://nbisweden.github.io/workshop-RaukR-2206/)