lecture5.Rmd

---
title: "Data science and analysis in Neuroscience"
author: "Kevin Allen"
date: "December 10, 2020"
output:
  ioslides_presentation: default
  beamer_presentation: default
  slidy_presentation: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(global.par = TRUE)
library(tidyverse)
library(knitr)
```

## A brief introduction to machine learning

1. Why should you care about machine learning?
2. Definition
4. Prediction versus inference
5. Supervised versus unsupervised
6. Regression versus classification
7. Instance-based versus model-based learning
8. Trainind and testing set
9. **Quizz!**
10. Example of a learning loop with a linear regression

## Objective

Our aim in the course is to understand what machine learning is and experiment with a few examples. 

## Why should you care about machine learning

It provides very useful tools to scientists.

* Track animal pose in a video, ([Deeplabcut](http://www.mousemotorlab.org/deeplabcut))
* Image segmentation of biological images, ([U-net](https://lmb.informatik.uni-freiburg.de/people/ronneber/u-net/))
* Clustering of spike waveforms ([Kilosort](https://github.com/MouseLand/Kilosort))
* 3D structure of proteins ([alphafold2](https://deepmind.com/blog/article/alphafold-a-solution-to-a-50-year-old-grand-challenge-in-biology))

Saves time, improves performance, and makes new experiments possible.

Very good libraries (packages) available.

## Definition of machine learning

Machine learning is the field of study that gives computers the ability to learn without being explicitly programmed.

-- Arthur Samuel, 1959

The computer learns from **input data** to achieve a specific objective.

Examples : A program learns to decide whether an email is spam or not based on training data. 

## Definition of machine learning

* Input: $X$ (can be single number or a matrix)
* Output: $Y$ (can be a single number or a matrix)
* Unknown function or model: $f()$
* Random error: $\epsilon$

<center>
$Y = f(X) + \epsilon$
</center>

<br>

Machine learning refers to a set of approaches for estimating the best parameters in $f$.

$f$ can be the equation of a line, a deep neural network, etc. It is your **model**.

## What is learning?

Learning can be defined as finding the best model parameters to solve a problem.

**Simple example**: Find the relationship between IQ and education with a linear regression model. Two parameters.

$$y = a*x + b $$ 

**Complex example**: Find a mouse in an image. Millions of parameters.

```{r, echo=FALSE,out.width = "350px"}
knitr::include_graphics("images/deep-neural-network.png")
```

## What is learning?

It is often possible to apply similar learning procedures for simple and complex models.

**Learning** or **training loop**

1. Start with random model parameters
2. Feed data with label to your model
3. Calculate the error of your model (loss).
4. Adjust the model parameters by a small amount to reduce the error.
5. Go back to 2.

## Prediction versus inference

Why do we want to estimate $f$?

<center>
$Y = f(X)$
</center>

### Prediction
* We focus on predicting $Y$.
* $f$ is treated as a black box (a useful black box)

### Inference
* **Understand** how $Y$ is affected as $X$ changes.
* Which predictors are associated with the response?
* Is the relation between $Y$ and each predictor adequately summarized using a linear equation?

## Supervised versus unsupervised

### Supervised
* The training set contains labeled data.
* For each observation of the predictors $X_{i}, i = 1,...,n$ there is a known response measurement $y_{i}$.
* Example: linear regression

### Unsupervised
* Uncovering hidden patterns from unlabeled data.
* For each observation $i = 1,...,n$, we observed a vector of measurements $X_{i}$, but no response $y_{i}$.
* Example: cluster analysis

## Regression versus classification

* If $Y$ is a continuous variable, then it is a regression task.
* If $Y$ is a categorical variable, then it is a classification task.

## Training and test sets

A **training set** is our observed data points that is used to estimate $f$. Our training set has $n$ observations.

A **test set** is used to test how accurate our model is. Not used for training!

Keeping data aside to test how well you model work is essential when using complex models. 

Complex models can learn to perform great on your training set but might generalize very poorly to new data. This is called **overfitting**.


## Time for a quizz!

[Link](https://docs.google.com/forms/d/e/1FAIpQLSfmL_igF1P0sZ_6aorGTE71pwNEa34oSWklG34y5vMPXvEYTQ/viewform?usp=sf_link)

or

https://tinyurl.com/y3jhxgr6

You have 5 minutes to complete the questions.

## How do computers learn?

Often using an iterative process (i.e. a loop).

1. Feed data to your model
2. Calculate the error (loss).
3. Adjust the model parameters by a small amount to minimize the loss.
4. Go back to 1.

In the next slides, we will implement this learning process.

## Our task: 

Write a learning loop to find the best parameters for a linear regression model.

We have bought a new thermometer for the lab but it does not show the units. 
We have an old thermometer that measure temperature in Celsius.
To understand the output of the new thermometer, we record pairs of readings from the old and new thermometers.

```{r data}
tc <- c(0.5, 14.0, 15.0, 28.0, 11.0, 8.0, 3.0, -4.0, 6.0, 13.0, 21.0)
tu <- c(35.7, 55.9, 58.2, 81.9, 56.3, 48.9, 33.9, 21.8, 48.4, 60.4, 68.4)
df <- data.frame(tc = tc,
           tu = tu)
```

Write and run the code as we go!

## Always plot the data first

```{r plotfirst, fig.width=4, fig.height=3}
df %>% ggplot()+
  geom_point(mapping = aes(x = tu, y = tc))
```

When we use our new thermometer, we want to know what would the value be in Celsius.


## Note

There are easier ways with R to solve this problem (e.g. `lm`). 

Today we focus on how computers can estimate the best parameters in a model.

## Choosing a model to solve our problem

The scatter plot suggests that there is a linear relationship.

What transformation do we need to do to go from $t_u$ to $t_c$?

```{r model}
model <- function(tu, params){
  w = params[1]
  b = params[2]
  return(w * tu + b)
}
```

$w$ and $b$ are for weight and bias

They are the slope and intercept in this case.

Our task is to find the best value of $w$ and $b$.

## Using our model

Our model returns predictions.

```{r use_model}
params <- c(1.0,0.0) # set arbitrary parameter values
model(tu,params) # make predictions
tc # real values
```

## Measuring the error of our model

To improve our model, we need to estimate the error in the prediction of our model.

This is called a **loss function**.

We usually compare the prediction of the model to the known label.

## Loss function

For linear regressions, we usually use the sum of the squared difference between the predicted and observed values.

```{r lm,echo=FALSE,fig.width = 6, fig.height = 4.5}
dm<-df
# fit a linear model to the data
fit<-lm(tc~tu, data = dm)
# get predicted values and residuals
dm$predicted<-predict(fit)
dm$residuals<-residuals(fit)
# plot 
ggplot(data=dm, mapping = aes(x=tu,y=tc))+
  geom_segment(aes(xend = tu, yend = predicted),alpha = 0.2) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "red", formula = y~x) 
```

## Loss function

Let's write our loss function.

```{r loss}
loss_fn <- function(y,y_pred){
  squared_diffs = (y-y_pred)^2
  return(squared_diffs)
}
```

## Loss function

The loss increases exponentially as the difference between the predicted and observed value increases. This penalizes large prediction errors.

```{r plotloss, fig.width = 4, fig.height = 2.5}
data.frame(predicted = seq(-10,10,0.1),
          observed = 0) %>% 
  ggplot(mapping = aes(x =predicted-observed,
                       y = loss_fn(predicted,observed)))+
  geom_point()
```

## Loss function

We need to find the model parameters ($w$ and $b$) to minimize the loss function.

## Loss function

Small modification: return one value when we pass several pairs of number as input.

We will output the mean of the squared difference
```{r loss2}
loss_fn <- function(y,y_pred){
  squared_diffs = (y-y_pred)^2
  return(mean(squared_diffs))
}
loss_fn(seq(-6,6,0.1),0)
```


## Predictions and loss calculation

```{r model_loss}
tu # new thermometer values
params <- c(1.0,0.0) # arbitrary values
tp <- model(tu,params)
tp # model prediction
loss <- loss_fn(tp,tc)
print(paste("loss:",loss))
```
## How do we adjust $w$ and $b$ to reduce the loss?

We could test many random values for $w$ and $b$ and find the combination with the lowest loss.
This solution will not scale well to models with millions of parameters.

A more efficient approach is called **gradient descent**. 

* Calculate the rate of change of the loss with respect to each parameter.
* Modify each parameter in the direction of decreasing loss.

## Rate of change of the loss relative to $w$

```{r par,echo=F}
par(bty = 'n', mar=c(4,4,1,0.1), mgp=c(2,1,0)) 
```

Here is the loss for different values of $w$, while keeping $b$ at `r params[2]`.

```{r w_ratechange, echo = FALSE, fig.width = 3, fig.height = 3}
ws <- seq(-3,3,0.1)
get_loss <- function(w, model, params, tu, tc){
  params[1] <- w
  return(loss_fn(model(tu,params),tc))
}
all_losses<-unlist(lapply(ws,get_loss,model,params,tu,tc))
data.frame(w = ws,
           loss = all_losses) %>% 
  ggplot(mapping = aes(x = w, y = loss))+
  geom_point()
```

We want adjust $w$ in steps so that we end up at the lowest point. 


## Derivative of the loss with respect to a parameter

The rate of change of the loss is the derivative of the loss with respect to a parameter.

The package called `torch` can calculate the derivative for you.

Torch is a c library that can calculate gradients. 

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

##  Torch

We need to create torch tensor (vector) with the argument `requires_grad = TRUE`.

```{r torch1}
params <- torch_tensor(c(1.0,0.0), requires_grad = TRUE)
params
class(params)
```
## Adjust the loos function for torch parameters

We need to slightly modify our loss function to work with torch_tensor
```{r adjustloss}
loss_fn <- function(y,y_pred){
  squared_diffs = (y-y_pred)^2
  return(squared_diffs$mean()) # mean is then with $ sign
}
```

## Calculate the gradients

```{r torch2}
params <- torch_tensor(c(1.0,0.0), requires_grad = TRUE)
loss <- loss_fn(model(tu,params),tc)
loss$backward() # calculate the gradients
params$grad
```

## Adjust our parameters

```{r adjust, results=F}
print(params) # print parameters

learning_rate <- 0.0001 # determine the size of the adjustment
 
with_no_grad({ # update the parameters
params$sub_(learning_rate * params$grad)
})

params$grad$zero_() # reset the gradients to 0
```


## So far we have...

* Input data with labels ($X$ and $Y$).
* A model with 2 parameters ($w$ and $b$).
* A loss function.
* torch to calculate the gradients of the parameters.

### We are ready to train our model!

## Training loop (minimal)


```{r trainingLoop}
trainingLoop <- function(n_epochs=10,learning_rate=0.0001
                         ,model,params,tu,tc){
  
  for (epoch in seq(1,n_epochs)){ # loop
    tp <- model(tu,params) # predict
    loss <- loss_fn(tp,tc) # loss
    loss$backward() # calculate the gradients
    with_no_grad({   # adjust parameters
      params$sub_(learning_rate * params$grad)
      })
    params$grad$zero_() # reset to 0.
    } # end of the loop
  
  return()
}
```

## Train the model (minimal)

```{r training}
params <- torch_tensor(c(1.0,0.0), requires_grad = TRUE)
trainingLoop(n_epochs = 10, learning_rate = 0.0001,
             model=model, params=params, tu = tu,tc=tc)
params
```

It is very hard to know what is going on!


## Training loop (with output information)

```{r trainingLoop1}
trainingLoop <- function(n_epochs = 10, learning_rate = 0.0001, 
                         model, params, tu, tc) {
  all_losses <- vector(length = n_epochs)
  for (epoch in seq(1, n_epochs)) {
    tp <- model(tu, params) # predict
    loss <- loss_fn(tp, tc) # loss
    loss$backward() # calculate the gradients
    with_no_grad({
      # adjust parameters
      params$sub_(learning_rate * params$grad)
    })
    if (epoch < 3 | epoch > n_epochs - 2)
      print_training_info(epoch, n_epochs, loss, params)
    params$grad$zero_() # reset to 0.
    all_losses[epoch] <- as.array(loss)
  }
  return(all_losses)
}
```


## Printing information in our loop

```{r defprint}
print_training_info <- function(epoch, n_epochs,loss, params){
  loss_r <- round(as_array(loss),3)
  params_r <- round(as_array(params),3)
  grads_r <- round(as_array(params$grad),3)
  
  print(paste("Epoch:", epoch, "of", n_epochs , ", Loss:", loss_r))
  print(paste("Gradients:", grads_r[1], grads_r[2], 
              " Parameters:", params_r[1],params_r[2]))
}
```

## Train the model

```{r training1}
params <- torch_tensor(c(1.0,0.0), requires_grad = TRUE)
losses <- trainingLoop(n_epochs = 10, learning_rate = 0.0001,
                       model=model, params=params, tu = tu,tc=tc)
losses
```
## Normalizing our input data

The gradients had a very different magnitude.

When using gradient descent, it is usually a good idea to normalize our input data.

```{r norm}
summary(tu)
tun <- tu * 0.1
summary(tun)
```
## Training with normalize data

```{r trainnorm,cache=F}
params <- torch_tensor(c(1.0,0.0), requires_grad = TRUE)
losses <- trainingLoop(n_epochs = 100, learning_rate = 0.01,
                       model=model, params=params, tu = tun,tc=tc)
```

## Plot the losses

We are still improving after 100 epochs.

```{r plotlosses2,fig.height = 3.5, fig.width = 4}
plot(losses,type="l",xlab="epoch",ylab="loss")
```

## Training for 3000 epochs

```{r trainnorm1, cache=F}
params <- torch_tensor(c(1.0,0.0), requires_grad = TRUE)
losses <- trainingLoop(n_epochs = 3000, learning_rate = 0.01,model=model, params=params, tu = tun,tc=tc)
```
Expected parameters: `w = 5.5556` and `b = -17.7778`. Not bad!

## Plot the losses

```{r plotlosses}
plot(losses,type="l",xlab="epoch",ylab="loss")
```

## Plot the output of the trained model

```{r plotmodel, fig.height = 2.5, fig.width=3}
df2 <- df
df2$tp <- as_array(model(tun,params)) 
ggplot(data=df2)+
  geom_point(mapping = aes(x = tu, y = tc))+
  geom_point(mapping = aes(x = tu, y = tp),color = "red")+
  xlab("Fahrenheit")+
  ylab("Celcius")

```


## Plot the output of the trained model

```{r plotmodel1, fig.height = 2.5, fig.width=3}
df2 <- df
df2$tp <- as.array(model(tun,params)) 
ggplot(data=df2)+
  geom_point(mapping = aes(x = tu, y = tc))+
  geom_point(mapping = aes(x = tu, y = tp),color = "red")+
  geom_smooth(mapping = aes(x = tu, y = tc), method = "lm", formula =y~x) +
  xlab("Fahrenheit")+
  ylab("Celcius")

```

## In summary

* We had a dataset and asked how to transform $X$ in order to obtain $Y$. 
* We chose a linear regression model with 2 parameters.
* We defined a loss function.
* We use "torch" to calculate the gradients.
* We trained our model by minimizing the loss function.

Very similar procedures are used to train more complex models (deep neural network).

## Linear regression

* In real life, use the function `lm()` to find the regression line (best fit). 


## Compare our results to lm()

```{r best_fit4}
myFittedModel<-lm(formula = tc~tun, data=df)
myFittedModel
```
Use `summary(myFittedModel)` to know if the model is significant. 


## Best libraries for machine learning

This is where python has the upper hand on R.

* Scikit-Learn (python, no deep neural networks)
* PyTorch (python, deep neural network)
* TensorFlow (python, deep neural network)

But R is catching up.

* caret package
* reticulate package
* torch package
* [datacamp](https://www.datacamp.com/tracks/machine-learning-scientist-with-r)

## Online courses

There are several excellent machine learning courses online (datacamp)

Many of them use python as language (especially for deep neural network).


## Good books

* [Deep Learning with PyTorch](https://pytorch.org/assets/deep-learning/Deep-Learning-with-PyTorch.pdf)

* [Hands-On Machine Learning with Scikit-Learn and TensorFlow](https://www.amazon.de/Hands-Machine-Learning-Scikit-Learn-TensorFlow/dp/1491962291)

* [Deep Learning for Coders with fastai and PyTorch](https://www.fast.ai/)

* [An Introduction to Statistical Learning: With Applications in R](https://www.academia.edu/36691506/An_Introduction_to_Statistical_Learning_Springer_Texts_in_Statistics_An_Introduction_to_Statistical_Learning) 

## For next week

* Read a book chapter: [The Machine Learning Landscape](https://www.oreilly.com/library/view/hands-on-machine-learning/9781491962282/ch01.html) (Hands on machine learning with sklearn and tensorflow, Chapter 1)

* Read a Nature Neuroscinece paper: [Deeplabcut](http://orga.cvss.cc/wp-content/uploads/2019/05/Mathis-etal-2018-NatureNeuroscience.pdf) Mathis et al., 2018.

* Have a look at the [DeepLabCut repository](https://github.com/AlexEMG/DeepLabCut) (https://github.com/AlexEMG/DeepLabCut)