C1-TensorFlow.qmd

---
output: html_document
editor_options:
  chunk_output_type: console
---

```{r}
#| echo: false
#| include: false
#| results: false
reticulate::use_condaenv("r-reticulate")
library(tensorflow)
tf
tf$abs(3.)
```

# Introduction to TensorFlow and Keras {#sec-tensorflowintro}

## Introduction to TensorFlow

One of the most commonly used frameworks for machine learning is **TensorFlow**. TensorFlow is an open source <a href="https://en.wikipedia.org/wiki/Linear_algebra" target="_blank" rel="noopener">linear algebra</a> library with focus on neural networks, published by Google in 2015. TensorFlow supports several interesting features, in particular automatic differentiation, several gradient optimizers and CPU and GPU parallelization.

These advantages are nicely explained in the following video:

```{r chunk_chapter3_43, eval=knitr::is_html_output(excludes = "epub"), results = 'asis', echo = F}
cat(
  '<iframe width="560" height="315"
  src="https://www.youtube.com/embed/MotG3XI2qSs"
  frameborder="0" allow="accelerometer; autoplay; encrypted-media;
  gyroscope; picture-in-picture" allowfullscreen>
  </iframe>'
)
```

To sum up the most important points of the video:

-   TensorFlow is a math library which is highly optimized for neural networks.
-   If a GPU is available, computations can be easily run on the GPU but even on a CPU TensorFlow is still very fast.
-   The "backend" (i.e. all the functions and all computations) are written in C++ and CUDA (CUDA is a programming language for NVIDIA GPUs).
-   The interface (the part of TensorFlow we use) is written in Python and is also available in R, which means, we can write the code in R/Python but it will be executed by the (compiled) C++ backend.

All operations in TensorFlow are written in C++ and are highly optimized. But don't worry, we don't have to use C++ to use TensorFlow because there are several bindings for other languages. TensorFlow officially supports a Python API, but meanwhile there are several community carried APIs for other languages:

-   R
-   Go
-   Rust
-   Swift
-   JavaScript

In this course we will use TensorFlow with the <a href="https://tensorflow.rstudio.com/" target="_blank" rel="noopener">https://tensorflow.rstudio.com/</a> binding, that was developed and published 2017 by the RStudio team. First, they developed an R package (reticulate) for calling Python in R. Actually, we are using the Python TensorFlow module in R (more about this later).

TensorFlow offers different levels of API. We could implement a neural network completely by ourselves or we could use Keras which is provided as a submodule by TensorFlow. Keras is a powerful module for building and training neural networks. It allows us building and training neural networks in a few lines of codes. Since the end of 2018, Keras and TensorFlow are completly interoperable, allowing us to utilize the best of both. In this course, we will show how we can use Keras for neural networks but also how we can use the TensorFlow's automatic differenation for using complex objective functions.

Useful links:

-   <a href="https://www.tensorflow.org/api_docs/python/tf" target="_blank" rel="noopener">TensorFlow documentation</a> (This is for the Python API, but just replace the "." with "\$".)
-   <a href="https://tensorflow.rstudio.com/" target="_blank" rel="noopener">Rstudio TensorFlow website</a>

### Data Containers

TensorFlow has two data containers (structures):

-   constant (tf\$constant): Creates a constant (immutable) value in the computation graph.
-   variable (tf\$Variable): Creates a mutable value in the computation graph (used as parameter/weight in models).

To get started with TensorFlow, we have to load the library and check if the installation worked.

```{r chunk_chapter3_44}
library(tensorflow)
library(keras)

# Don't worry about weird messages. TensorFlow supports additional optimizations.
exists("tf")

immutable = tf$constant(5.0)
mutable = tf$constant(5.0)
```

Don't worry about weird messages (they will only appear once at the start of the session).

### Basic Operations

We now can define the variables and do some math with them:

```{r chunk_chapter3_45}
a = tf$constant(5)
b = tf$constant(10)
print(a)
print(b)
c = tf$add(a, b)
print(c)
tf$print(c) # Prints to stderr. For stdout, use k_print_tensor(..., message).
k_print_tensor(c) # Comes out of Keras!
```

Normal R methods such as print() are provided by the R package "tensorflow".

The TensorFlow library (created by the RStudio team) built R methods for all common operations:

```{r chunk_chapter3_46}
`+.tensorflow.tensor` = function(a, b){ return(tf$add(a,b)) }
# Mind the backticks.
k_print_tensor(a+b)
```

Their operators also automatically transform R numbers into constant tensors when attempting to add a tensor to an R number:

```{r chunk_chapter3_47}
d = c + 5  # 5 is automatically converted to a tensor.
print(d)
```

TensorFlow containers are objects, what means that they are not just simple variables of type numeric (class(5)), but they instead have so called methods. Methods are changing the state of a class (which for most of our purposes here is the values of the object). For instance, there is a method to transform the tensor object back to an R object:

```{r chunk_chapter3_48}
class(d)
class(d$numpy())
```

### Data Types

R uses dynamic typing, what means you can assign a number, character, function or whatever to a variable and the the type is automatically inferred. In other languages you have to state the type explicitly, e.g. in C:

```{C, eval=F}
int a = 5;
float a = 5.0;
char a = "a";
```

While TensorFlow tries to infer the type dynamically, you must often state it explicitly. Common important types:

-   float32 (floating point number with 32 bits, "single precision")
-   float64 (floating point number with 64 bits, "double precision")
-   int8 (integer with 8 bits)

The reason why TensorFlow is so explicit about types is that many GPUs (e.g. the NVIDIA GeForces) can handle only up to 32 bit numbers! (you do not need high precision in graphical modeling)

But let us see in practice what we have to do with these types and how to specifcy them:

```{r chunk_chapter3_49, eval=FALSE}
r_matrix = matrix(runif(10*10), 10, 10)
m = tf$constant(r_matrix, dtype = "float32") 
b = tf$constant(2.0, dtype = "float64")
c = m / b # Doesn't work! We try to divide float32/float64.
```

So what went wrong here? We tried to divide a float32 by a float64 number, but we can only divide numbers of the same type!

```{r chunk_chapter3_50, eval=TRUE}
r_matrix = matrix(runif(10*10), 10, 10)
m = tf$constant(r_matrix, dtype = "float64")
b = tf$constant(2.0, dtype = "float64")
c = m / b # Now it works.
```

We can also specify the type of the object by providing an object e.g. tf\$float64.

```{r chunk_chapter3_51}
r_matrix = matrix(runif(10*10), 10, 10)
m = tf$constant(r_matrix, dtype = tf$float64)
```

In TensorFlow, arguments often require exact/explicit data types: TensorFlow often expects integers as arguments. In R however an integer is normally saved as float. Thus, we have to use an "L" after an integer to tell the R interpreter that it should be treated as an integer:

```{r chunk_chapter3_52, eval=FALSE}
is.integer(5)
is.integer(5L)
matrix(t(r_matrix), 5, 20, byrow = TRUE)
tf$reshape(r_matrix, shape = c(5, 20))$numpy()
tf$reshape(r_matrix, shape = c(5L, 20L))$numpy()
```

Skipping the "L" is one of the most common errors when using R-TensorFlow!

### Exercises

::: {.callout-caution icon="false"}
#### Question: TensorFlow Operations

To run TensorFlow from R, note that you can access the different mathematical operations in TensorFlow via tf\$..., e.g. there is a tf\$math\$... for all common math operations or the tf\$linalg\$... for different linear algebra operations. Tip: type tf\$ and then hit the tab key to list all available options (sometimes you have to do this directly in the console).

An example: How to get the maximum value of a vector?

An example: How to get the maximum value of a vector?

```{r chunk_chapter3_task_5, eval=FALSE}
library(tensorflow)
library(keras)

x = 100:1
y = as.double(100:1)

max(x)  # R solution. Integer!
tf$math$reduce_max(x) # TensorFlow solution. Integer!

max(y)  # Float!
tf$math$reduce_max(y) # Float!
```

Rewrite the following expressions (a to g) in TensorFlow:

```{r chunk_chapter3_task_6, eval=TRUE}
x = 100:1
y = as.double(100:1)

# a)
min(x)

# b)
mean(x)

# c) Tip: Use Google!
which.max(x)

# d) 
which.min(x)

# e) Tip: Use Google! 
order(x)

# f) Tip: See tf$reshape.
m = matrix(y, 10, 10) # Mind: We use y here! (Float)
m_2 = abs(m %*% t(m))  # m %*% m is the normal matrix multiplication.
m_2_log = log(m_2)
print(m_2_log)

# g) Custom mean function i.e. rewrite the function using TensorFlow. 
mean_R = function(y){
  result = sum(y) / length(y)
  return(result)
}

mean_R(y) == mean(y)	# Test for equality.
```

`r hide("Click here to see the solution")`

```{r chunk_chapter3_task_7, eval=TRUE, include=TRUE}
library(tensorflow)
library(keras)

x = 100:1
y = as.double(100:1)

# a)    min(x)
tf$math$reduce_min(x) # Integer!
tf$math$reduce_min(y) # Float!

# b)    mean(x)
# Check out the difference here:
mean(x)
mean(y)
tf$math$reduce_mean(x)  # Integer!
tf$math$reduce_mean(y)  # Float!

# c)    which.max(x)
tf$argmax(x)
tf$argmax(y)

# d)    which.min(x)
tf$argmin(x)

# e)    order(x)
tf$argsort(x)

# f)
# m = matrix(y, 10, 10)
# m_2 = abs(m %*% m)
# m_2_log = log(m_2)

# Mind: We use y here! TensorFlow just accepts floats in the following lines!
mTF = tf$reshape(y, list(10L, 10L))
m_2TF = tf$math$abs( tf$matmul(mTF, tf$transpose(mTF)) )
m_2_logTF = tf$math$log(m_2TF)
print(m_2_logTF)

# g)    # Custom mean function
mean_TF = function(y){
  result = tf$math$reduce_sum(y)
  return( result / length(y) )  # If y is an R object.
}
mean_TF(y) == mean(y)
```

`r unhide()`
:::

::: {.callout-caution icon="false"}
#### Question: Runtime

This exercise compares the speed of R to TensorFlow. The first exercise is to rewrite the following function in TensorFlow:

```{r chunk_chapter3_task_8, eval=TRUE}
do_something_R = function(x = matrix(0.0, 10L, 10L)){
  mean_per_row = apply(x, 1, mean)
  result = x - mean_per_row
  return(result)
}
```

Here, we provide a skeleton for a TensorFlow function:

```{r chunk_chapter3_task_9, eval=FALSE, purl=FALSE}
do_something_TF = function(x = matrix(0.0, 10L, 10L)){
   ...
}
```

We can compare the speed using the Microbenchmark package:

```{r chunk_chapter3_task_10, eval=FALSE, purl=FALSE}
test = matrix(0.0, 100L, 100L)
microbenchmark::microbenchmark(do_something_R(test), do_something_TF(test))
```

Try different matrix sizes for the test matrix and compare the speed.

Tip: Have a look at the the tf.reduce_mean documentation and the "axis" argument.

<br/>

Compare the following with different matrix sizes:

-   test = matrix(0.0, 1000L, 500L)
-   testTF = tf\$constant(test)

Also try the following:

```{r chunk_chapter3_task_11, eval=FALSE, purl=FALSE}
microbenchmark::microbenchmark(
   tf$matmul(testTF, tf$transpose(testTF)), # TensorFlow style.
   test %*% t(test)  # R style.
)
```

`r hide("Click here to see the solution")`

```{r chunk_chapter3_task_12, eval=TRUE, include=TRUE}
do_something_TF = function(x = matrix(0.0, 10L, 10L)){
  x = tf$constant(x)  # Remember, this is a local copy!
  mean_per_row = tf$reduce_mean(x, axis = 1L)
  result = tf$transpose(x) - mean_per_row 
  return(tf$transpose(result)) 
}
```

Short explanation: A 2 dimensional tensor in tensorflow can be seen as an array of arrays. E.g. a 20x30 (rows x columns) matrix is an array which consists of 20 arrays of length 30. If we write result = matrix - vector, tensorflow tries to subtract from every inner array of the matrix (so from every row) the vector. Therefore if the want to substract the first number of mean_per_row from the first column of the matrix, we have to transpose this matrix, do the subtraction and then of course transpose it back to the original dimensions. Alternativly one could transform the vector to an matrix of size 20x1: `result = x - tf$reshape(mean_per_row, shape = list(20L,1L))` and then omit the second transpose in the return statement. You could now also try to guess, which one of the two is faster and then check that!

```{r chunk_chapter3_task_13, eval=TRUE, include=TRUE}
test = matrix(0.0, 100L, 100L)
microbenchmark::microbenchmark(do_something_R(test), do_something_TF(test))

test = matrix(0.0, 1000L, 500L)
microbenchmark::microbenchmark(do_something_R(test), do_something_TF(test))
```

Why is R faster (the first time)?

-   

    a)  The R functions we used (apply, mean, "-") are also implemented in C.

-   

    b)  The problem is not large enough and TensorFlow has an overhead.

<br/>

```{r chunk_chapter3_task_14, eval=TRUE, include=TRUE}
test = matrix(0.0, 1000L, 500L)
testTF = tf$constant(test)

microbenchmark::microbenchmark(
  tf$matmul(testTF, tf$transpose(testTF)),  # TensorFlow style.
  test %*% t(test) # R style.
)
```

`r unhide()`
:::

::: {.callout-caution icon="false"}
#### Question: Linear Algebra

Google to find out how to write the following expressions in TensorFlow:

```{r chunk_chapter3_task_15, eval=TRUE}
A = matrix(c(1, 2, 0, 0, 2, 0, 2, 5, 3), 3, 3)

# i)
solve(A)  # Solve equation AX = B. If just A  is given, invert it.

# j)
diag(A) # Diagonal of A, if no matrix is given, construct diagonal matrix.

# k)
diag(diag(A)) # Diagonal matrix with entries diag(A).

# l)
eigen(A)

# m)
det(A)
```

`r hide("Click here to see the solution")`

```{r chunk_chapter3_task_16, eval=TRUE, include=TRUE}
library(tensorflow)
library(keras)

A = matrix(c(1., 2., 0., 0., 2., 0., 2., 5., 3.), 3, 3)
# Do not use the "L" form here!

# i)    solve(A)
tf$linalg$inv(A)

# j)    diag(A)
tf$linalg$diag_part(A)

# k)    diag(diag(A))
tf$linalg$diag(tf$linalg$diag_part(A))

# l)    eigen(A)
tf$linalg$eigh(A)

# m)    det(A)
tf$linalg$det(A)
```

`r unhide()`
:::

::: {.callout-caution icon="false"}
#### Question: Automatic differentation

TensorFlow supports automatic differentiation (analytical and not numerical!). Let's have a look at the function $f(x) = 5 x^2 + 3$ with derivative $f'(x) = 10x$. So for $f'(5)$ we will get $10$.

Let's do this in TensorFlow. Define the function:

```{r chunk_chapter3_task_17, eval=TRUE}
f = function(x){ return(5.0 * tf$square(x) + 3.0) }
```

We want to calculate the derivative for $x = 2.0$:

```{r chunk_chapter3_task_18, eval=TRUE}
x = tf$constant(2.0)
```

To do automatic differentiation, we have to forward $x$ through the function within the tf\$GradientTape() environment. We have also have to tell TensorFlow which value to "watch":

```{r chunk_chapter3_task_19, eval=TRUE}
with(tf$GradientTape() %as% tape,
  {
    tape$watch(x)
    y = f(x)
  }
)
```

To print the gradient:

```{r chunk_chapter3_task_20, eval=TRUE}
(tape$gradient(y, x))
```

We can also calculate the second order derivative $f''(x) = 10$:

```{r chunk_chapter3_task_21, eval=TRUE}
with(tf$GradientTape() %as% first,
  {
    first$watch(x)
    with(tf$GradientTape() %as% second,
      {
        second$watch(x)
        y = f(x)
        g = first$gradient(y, x)
      }
    )
  }
)

(second$gradient(g, x))
```

What is happening here? Think about and discuss it.

A more advanced example: *Linear regression*

In this case we first simulate data following $\boldsymbol{y} = \boldsymbol{X} \boldsymbol{w} + \boldsymbol{\epsilon}$ ($\boldsymbol{\epsilon}$ follows a normal distribution == error).

```{r chunk_chapter3_task_22, eval=TRUE}
set_random_seed(321L, disable_gpu = FALSE)	# Already sets R's random seed.

x = matrix(round(runif(500, -2, 2), 3), 100, 5)
w = round(rnorm(5, 2, 1), 3)
y = x %*% w + round(rnorm(100, 0, 0.25), 4)
```

In R we would do the following to fit a linear regression model:

```{r chunk_chapter3_task_23, eval=TRUE}
summary(lm(y~x))
```

Let's build our own model in TensorFlow. Here, we use now the variable data container type (remember they are mutable and we need this type for the weights ($\boldsymbol{w}$) of the regression model). We want our model to learn these weights.

The input (predictors, independent variables or features, $\boldsymbol{X}$) and the observed (response, $\boldsymbol{y}$) are constant and will not be learned/optimized.

```{r chunk_chapter3_task_24, eval=TRUE}
library(tensorflow)
library(keras)
set_random_seed(321L, disable_gpu = FALSE)	# Already sets R's random seed.

x = matrix(round(runif(500, -2, 2), 3), 100, 5)
w = round(rnorm(5, 2, 1), 3)
y = x %*% w + round(rnorm(100, 0, 0.25), 4)

# Weights we want to learn.
# We know the real weights but in reality we wouldn't know them.
# So use guessed ones.
wTF = tf$Variable(matrix(rnorm(5, 0, 0.01), 5, 1))

xTF = tf$constant(x)
yTF = tf$constant(y)

# We need an optimizer which updates the weights (wTF).
optimizer = tf$keras$optimizers$Adamax(learning_rate = 0.1)

for(i in 1:100){
  with(tf$GradientTape() %as% tape,
    {
      pred = tf$matmul(xTF, wTF)
      loss = tf$sqrt(tf$reduce_mean(tf$square(yTF - pred)))
    }
  )

  if(!i%%10){ k_print_tensor(loss, message = "Loss: ") }  # Every 10 times.
  grads = tape$gradient(loss, wTF)
  optimizer$apply_gradients(purrr::transpose(list(list(grads), list(wTF))))
}

k_print_tensor(wTF, message = "Resulting weights:\n")
cat("Original weights: ", w, "\n")
```

Discuss the code, go through the code line by line and try to understand it.

Additional exercise:

Play around with the simulation, increase/decrease the number of weights, add an intercept (you also need an additional variable in model).

`r hide("Click here to see the solution")`

```{r chunk_chapter3_task_25, eval=TRUE, include=TRUE, warning=FALSE}
library(tensorflow)
library(keras)
set_random_seed(321L, disable_gpu = FALSE)	# Already sets R's random seed.

numberOfWeights = 3
numberOfFeatures = 10000

x = matrix(round(runif(numberOfFeatures * numberOfWeights, -2, 2), 3),
           numberOfFeatures, numberOfWeights)
w = round(rnorm(numberOfWeights, 2, 1), 3)
intercept = round(rnorm(1, 3, .5), 3)
y = intercept + x %*% w + round(rnorm(numberOfFeatures, 0, 0.25), 4)

# Guessed weights and intercept.
wTF = tf$Variable(matrix(rnorm(numberOfWeights, 0, 0.01), numberOfWeights, 1))
interceptTF = tf$Variable(matrix(rnorm(1, 0, .5)), 1, 1) # Double, not float32.

xTF = tf$constant(x)
yTF = tf$constant(y)

optimizer = tf$keras$optimizers$Adamax(learning_rate = 0.05)

for(i in 1:100){
  with(tf$GradientTape(persistent = TRUE) %as% tape,
    {
      pred = tf$add(interceptTF, tf$matmul(xTF, wTF))
      loss = tf$sqrt(tf$reduce_mean(tf$square(yTF - pred)))
    }
  )

  if(!i%%10){ k_print_tensor(loss, message = "Loss: ") }  # Every 10 times.
  grads = tape$gradient(loss, list(wTF, interceptTF))
  optimizer$apply_gradients(purrr::transpose(list(grads, list(wTF, interceptTF))))
}

k_print_tensor(wTF, message = "Resulting weights:\n")
cat("Original weights: ", w, "\n")
k_print_tensor(interceptTF, message = "Resulting intercept:\n")
cat("Original intercept: ", intercept, "\n")
```

`r unhide()`
:::

## Introduction to PyTorch

PyTorch is another famous library for deep learning. Like TensorFlow, Torch itself is written in C++ with an API for Python. In 2020, the RStudio team released R-Torch, and while R-TensorFlow calls the Python API in the background, the R-Torch API is built directly on the C++ Torch library!

Useful links:

-   <a href="https://pytorch.org/docs/stable/index.html" target="_blank" rel="noopener">PyTorch documentation</a> (This is for the Python API, bust just replace the "." with "\$".)
-   <a href="https://torch.mlverse.org/" target="_blank" rel="noopener">R-Torch website</a>

To get started with Torch, we have to load the library and check if the installation worked.

```{r chunk_chapter3_53}
library(torch)
```

### Data Containers

Unlike TensorFlow, Torch doesn't have two data containers for mutable and immutable variables. All variables are initialized via the torch_tensor function:

```{r chunk_chapter3_54}
a = torch_tensor(1.)
```

To mark variables as mutable (and to track their operations for automatic differentiation) we have to set the argument 'requires_grad' to true in the torch_tensor function:

```{r chunk_chapter3_55}
mutable = torch_tensor(5, requires_grad = TRUE) # tf$Variable(...)
immutable = torch_tensor(5, requires_grad = FALSE) # tf$constant(...)
```

### Basic Operations

We now can define the variables and do some math with them:

```{r chunk_chapter3_56}
a = torch_tensor(5.)
b = torch_tensor(10.)
print(a)
print(b)
c = a$add(b)
print(c)
```

The R-Torch package provides all common R methods (an advantage over TensorFlow).

```{r chunk_chapter3_57}
a = torch_tensor(5.)
b = torch_tensor(10.)
print(a+b)
print(a/b)
print(a*b)
```

Their operators also automatically transform R numbers into tensors when attempting to add a tensor to a R number:

```{r chunk_chapter3_58}
d = a + 5  # 5 is automatically converted to a tensor.
print(d)
```

As for TensorFlow, we have to explicitly transform the tensors back to R:

```{r chunk_chapter3_59}
class(d)
class(as.numeric(d))
```

### Data Types

Similar to TensorFlow:

```{r chunk_chapter3_60, eval=FALSE}
r_matrix = matrix(runif(10*10), 10, 10)
m = torch_tensor(r_matrix, dtype = torch_float32()) 
b = torch_tensor(2.0, dtype = torch_float64())
c = m / b 
```

But here's a difference! With TensorFlow we would get an error, but with R-Torch, m is automatically casted to a double (float64). However, this is still bad practice!

During the course we will try to provide the corresponding PyTorch code snippets for all Keras/TensorFlow examples.

### Exercises

::: {.callout-caution icon="false"}
#### Question: Torch Operations

Rewrite the following expressions (a to g) in torch:

```{r chunk_chapter3_task_torch_6, eval=TRUE}
x = 100:1
y = as.double(100:1)

# a)
min(x)

# b)
mean(x)

# c) Tip: Use Google!
which.max(x)

# d) 
which.min(x)

# e) Tip: Use Google! 
order(x)

# f) Tip: See tf$reshape.
m = matrix(y, 10, 10) # Mind: We use y here! (Float)
m_2 = abs(m %*% t(m))  # m %*% m is the normal matrix multiplication.
m_2_log = log(m_2)
print(m_2_log)

# g) Custom mean function i.e. rewrite the function using TensorFlow. 
mean_R = function(y){
  result = sum(y) / length(y)
  return(result)
}

mean_R(y) == mean(y)	# Test for equality.
```

`r hide("Click here to see the solution")`

```{r chunk_chapter3_task_torch_7, eval=TRUE, include=TRUE}
library(torch)


x = 100:1
y = as.double(100:1)

# a)    min(x)
torch_min(x) # Integer!
torch_min(y) # Float!

# b)    mean(x)
# Check out the difference here:
mean(x)
mean(y)
torch_mean(torch_tensor(x, dtype = torch_float32()))  # Integer! Why?
torch_mean(y)  # Float!

# c)    which.max(x)
torch_argmax(x)
torch_argmax(y)

# d)    which.min(x)
torch_argmin(x)

# e)    order(x)
torch_argsort(x)

# f)
# m = matrix(y, 10, 10)
# m_2 = abs(m %*% m)
# m_2_log = log(m_2)

# Mind: We use y here! 
mTorch = torch_reshape(y, c(10, 10))
mTorch2 = torch_abs(torch_matmul(mTorch, torch_t(mTorch))) # hard to read!

# Better:
mTorch2 = mTorch$matmul( mTorch$t() )$abs()
mTorch2_log = mTorch$log()

print(mTorch2_log)

# g)    # Custom mean function
mean_Torch = function(y){
  result = torch_sum(y)
  return( result / length(y) )  # If y is an R object.
}
mean_Torch(y) == mean(y)
```

`r unhide()`
:::

::: {.callout-caution icon="false"}
#### Question: Runtime

1.  What is the meaning of "An effect is not significant"?
2.  Is an effect with three \*\*\* more significant / certain than an effect with one \*?

`r hide("Click here to see the solution")` This exercise compares the speed of R to torch The first exercise is to rewrite the following function in torch:

```{r chunk_chapter3_task_torch_8, eval=TRUE}
do_something_R = function(x = matrix(0.0, 10L, 10L)){
  mean_per_row = apply(x, 1, mean)
  result = x - mean_per_row
  return(result)
}
```

Here, we provide a skeleton for a TensorFlow function:

```{r chunk_chapter3_task_torch_9, eval=FALSE, purl=FALSE}
do_something_torch= function(x = matrix(0.0, 10L, 10L)){
   ...
}
```

We can compare the speed using the Microbenchmark package:

```{r chunk_chapter3_task_torch_10, eval=FALSE, purl=FALSE}
test = matrix(0.0, 100L, 100L)
microbenchmark::microbenchmark(do_something_R(test), do_something_torch(test))
```

Try different matrix sizes for the test matrix and compare the speed.

Tip: Have a look at the the torch_mean documentation and the "dim" argument.

<br/>

Compare the following with different matrix sizes:

-   test = matrix(0.0, 1000L, 500L)
-   testTorch = torch_tensor(test)

Also try the following:

```{r chunk_chapter3_task_torch_11, eval=FALSE, purl=FALSE}
microbenchmark::microbenchmark(
   torch_matmul(testTorch, testTorch$t()), # Torch style.
   test %*% t(test)  # R style.
)
```

`r hide("Click here to see the solution")`

```{r chunk_chapter3_task_torch_12, eval=TRUE, include=TRUE}
do_something_torch = function(x = matrix(0.0, 10L, 10L)){
  x = torch_tensor(x)  # Remember, this is a local copy!
  mean_per_row = torch_mean(x, dim = 1)
  result = x - mean_per_row
  return(result)
}
```

```{r chunk_chapter3_task_torch_13, eval=TRUE, include=TRUE}
test = matrix(0.0, 100L, 100L)
microbenchmark::microbenchmark(do_something_R(test), do_something_torch(test))

test = matrix(0.0, 1000L, 500L)
microbenchmark::microbenchmark(do_something_R(test), do_something_torch(test))
```

Why is R faster (the first time)?

-   

    a)  The R functions we used (apply, mean, "-") are also implemented in C.

-   

    b)  The problem is not large enough and torch has an overhead.

<br/>

```{r chunk_chapter3_task_torch_14, eval=TRUE, include=TRUE}
test = matrix(0.0, 1000L, 500L)
testTorch = torch_tensor(test)

microbenchmark::microbenchmark(
   torch_matmul(testTorch, testTorch$t()), # Torch style.
   test %*% t(test)  # R style.
)
```

`r unhide()`
:::

::: {.callout-caution icon="false"}
#### Question: Linear Algebra

Google to find out how to write the following tasks in torch:

```{r chunk_chapter3_task_torch_15, eval=TRUE}
A = matrix(c(1, 2, 0, 0, 2, 0, 2, 5, 3), 3, 3)

# i)
solve(A)  # Solve equation AX = B. If just A  is given, invert it.

# j)
diag(A) # Diagonal of A, if no matrix is given, construct diagonal matrix.

# k)
diag(diag(A)) # Diagonal matrix with entries diag(A).

# l)
eigen(A)

# m)
det(A)
```

`r hide("Click here to see the solution")`

```{r chunk_chapter3_task_torch_16, eval=TRUE, include=TRUE}
library(torch)

A = matrix(c(1., 2., 0., 0., 2., 0., 2., 5., 3.), 3, 3)
# Do not use the "L" form here!

# i)    solve(A)
linalg_inv(A)

# j)    diag(A)
torch_diag(A)

# k)    diag(diag(A))
torch_diag(A)$diag()

# l)    eigen(A)
linalg_eigh(A)

# m)    det(A)
linalg_det(A)
```

`r unhide()`
:::

::: {.callout-caution icon="false"}
#### Question: Automatic differentation

Torch supports automatic differentiation (analytical and not numerical!). Let's have a look at the function $f(x) = 5 x^2 + 3$ with derivative $f'(x) = 10x$. So for $f'(5)$ we will get $10$.

Let's do this in torch Define the function:

```{r chunk_chapter3_task_torch_17, eval=TRUE}
f = function(x){ return(5.0 * torch_pow(x, 2.) + 3.0) }
```

We want to calculate the derivative for $x = 2.0$:

```{r chunk_chapter3_task_torch_18, eval=TRUE}
x = torch_tensor(2.0, requires_grad = TRUE)
```

To do automatic differentiation, we have to forward $x$ through the function and call the \$backward() method of the result:

```{r chunk_chapter3_task_torch_19, eval=TRUE}
y = f(x)
y$backward(retain_graph=TRUE )

```

To print the gradient:

```{r chunk_chapter3_task_torch_20, eval=TRUE}
x$grad
```

We can also calculate the second order derivative $f''(x) = 10$:

```{r chunk_chapter3_task_torch_21, eval=TRUE}
x = torch_tensor(2.0, requires_grad = TRUE)
y = f(x)
grad = torch::autograd_grad(y, x, retain_graph = TRUE, create_graph = TRUE)[[1]] # first
(torch::autograd_grad(grad, x, retain_graph = TRUE, create_graph = TRUE)[[1]]) # second

```

What is happening here? Think about and discuss it.

A more advanced example: *Linear regression*

In this case we first simulate data following $\boldsymbol{y} = \boldsymbol{X} \boldsymbol{w} + \boldsymbol{\epsilon}$ ($\boldsymbol{\epsilon}$ follows a normal distribution == error).

```{r chunk_chapter3_task_torch_22, eval=TRUE}
set_random_seed(321L, disable_gpu = FALSE)	# Already sets R's random seed.

x = matrix(round(runif(500, -2, 2), 3), 100, 5)
w = round(rnorm(5, 2, 1), 3)
y = x %*% w + round(rnorm(100, 0, 0.25), 4)
```

In R we would do the following to fit a linear regression model:

```{r chunk_chapter3_task_torch_23, eval=TRUE}
summary(lm(y~x))
```

Let's build our own model in TensorFlow. Here, we use now the variable data container type (remember they are mutable and we need this type for the weights ($\boldsymbol{w}$) of the regression model). We want our model to learn these weights.

The input (predictors, independent variables or features, $\boldsymbol{X}$) and the observed (response, $\boldsymbol{y}$) are constant and will not be learned/optimized.

```{r chunk_chapter3_task_torch_24, eval=TRUE}
library(torch)
torch::torch_manual_seed(42L)

x = matrix(round(runif(500, -2, 2), 3), 100, 5)
w = round(rnorm(5, 2, 1), 3)
y = x %*% w + round(rnorm(100, 0, 0.25), 4)

# Weights we want to learn.
# We know the real weights but in reality we wouldn't know them.
# So use guessed ones.
wTorch = torch_tensor(matrix(rnorm(5, 0, 0.01), 5, 1), requires_grad = TRUE)

xTorch = torch_tensor(x)
yTorch = torch_tensor(y)

# We need an optimizer which updates the weights (wTF).
optimizer = optim_adam(params = list(wTorch), lr = 0.1)

for(i in 1:100){
  pred = xTorch$matmul(wTorch)
  loss = (yTorch - pred)$pow(2.0)$mean()$sqrt()

  if(!i%%10){ print(paste0("Loss: ", as.numeric(loss)))}  # Every 10 times.
  loss$backward()
  optimizer$step() # do optimization step
  optimizer$zero_grad() # reset gradients
}
cat("Inferred weights: ", round(as.numeric(wTorch), 3), "\n")
cat("Original weights: ", w, "\n")
```

Discuss the code, go through the code line by line and try to understand it.

Additional exercise:

Play around with the simulation, increase/decrease the number of weights, add an intercept (you also need an additional variable in model).

`r hide("Click here to see the solution")`

```{r chunk_chapter3_task_torch_25, eval=TRUE, include=TRUE}
library(torch)
torch::torch_manual_seed(42L)

numberOfWeights = 3
numberOfFeatures = 10000

x = matrix(round(runif(numberOfFeatures * numberOfWeights, -2, 2), 3),
           numberOfFeatures, numberOfWeights)
w = round(rnorm(numberOfWeights, 2, 1), 3)
intercept = round(rnorm(1, 3, .5), 3)
y = intercept + x %*% w + round(rnorm(numberOfFeatures, 0, 0.25), 4)

# Guessed weights and intercept.
wTorch = torch_tensor(matrix(rnorm(numberOfWeights, 0, 0.01), numberOfWeights, 1), requires_grad = TRUE)
interceptTorch = torch_tensor(matrix(rnorm(1, 0, .5), 1, 1), requires_grad = TRUE) # Double, not float32.

xTorch = torch_tensor(x)
yTorch = torch_tensor(y)

# We need an optimizer which updates the weights (wTF).
optimizer = optim_adam(params = list(wTorch, interceptTorch), lr = 0.1)

for(i in 1:100){
  pred = xTorch$matmul(wTorch)$add(interceptTorch)
  loss = (yTorch - pred)$pow(2.0)$mean()$sqrt()

  if(!i%%10){ print(paste0("Loss: ", as.numeric(loss)))}  # Every 10 times.
  loss$backward()
  optimizer$step() # do optimization step
  optimizer$zero_grad() # reset gradients
}
cat("Inferred weights: ", round(as.numeric(wTorch), 3), "\n")
cat("Original weights: ", w, "\n")

cat("Inferred intercept: ", round(as.numeric(interceptTorch), 3), "\n")
cat("Original intercept: ", intercept, "\n")
```

`r unhide()`
:::