torch_Python_R_ResNet.qmd

---
title: "_ResNet_ : comparaison `{torch}` et `pytorch`"
lang: fr
execute:
  cache: true
output: 
  html_document:
   toc: true
   toc_float: true
---

## Intro

Le but est de comparer la syntaxe utilisée par les deux langages, seuls les exemples en `R` sont exécutés (difficultés de faire un quarto avec deux langages, voir [issue ici](https://github.com/quarto-dev/quarto-cli/discussions/4232) ou de faire tourner `pytorch` avec `{reticulate}`). L'exemple d'application est un _ResNet_.

Pour charger installer la bibliothèque, en python, plusieurs possibilités sont données sur le [site web officiel](https://pytorch.org/get-started/locally/). Vous pouvez utiliser `pip` ou `conda`, avec ou sans `cuda`. Pour `R`, cela a déjà été fait à la session 2022 de [finistR](https://stateofther.github.io/finistR2022/autodiff.html).

```{r, eval = FALSE}
install.packages(torch)
torch::install_torch() ## si vous avez un GPU compatible avec CUDA
## torch::install_torch(type = "cpu") ## sinon
```

Un bon ouvrage pour démarrer : [Deep Learning and Scientific Computing with `R` `torch`](https://skeydan.github.io/Deep-Learning-and-Scientific-Computing-with-R-torch/).

## Appel des fonctions

En `python`, un premier exemple serait déjà d'importer la bibliothèque et de créer un `tensor` :

```python
import torch
t1 = torch.tensor(1)
t1
```

Donnera en `R` :
```{r}
library(torch)
t1 <- torch_tensor(1)
t1
```

On comprend que la convention de nommage garde `torch_` en préfix de toutes les fonctions implémentées.

## Construction d'un réseau de neurone, exemple d'un ResNet

### Chargement des données

On va utiliser `torchvision` pour importer un jeu de données connu (des images annotées). On veut appliquer des transformations sur le jeu de test, mais hélas tout n'est pas encore implémenté. On commente les transformations pas encore implementées.

```python
transform = transforms.Compose([
    #transforms.Pad(4),
    #transforms.RandomHorizontalFlip(),
    #transforms.RandomCrop(32),
    transforms.ToTensor()])
```

```{r}
transform <- function(img) {
  # img <- torchvision::transform_pad(img, 4) # pas implémenté
  # img <- torchvision::transform_random_horizontal_flip(img) # pas implémenté
  # img <- torchvision::transform_random_crop(img, 32) # bug sur la taille des images
  img <- torchvision::transform_to_tensor(img)
  return(img)
}
```

Voici le code pour charger les données en `python` et `R` :
```python
num_samples = 1000

trainData = torchvision.datasets.CIFAR10(root="./data",
                                         train=True,
                                         transform=transform,
                                         download=True)

testData = torchvision.datasets.CIFAR10(root="./data",
                                        train=False,
                                        transform=transforms.ToTensor())

trainLoader = torch.utils.data.DataLoader(
    dataset=Subset(trainData, range(num_samples)),
    batch_size=256,
    shuffle=True)

testLoader = torch.utils.data.DataLoader(
    dataset=Subset(testData, range(num_samples)),
    batch_size=256,
    shuffle=False)
```

```{r}
num_samples = 1000

train_data <- torchvision::cifar10_dataset(
  root = "./data",
  train = TRUE,
  transform = transform,
  download = TRUE
)

test_data <- torchvision::cifar10_dataset(
  root = "./data",
  train = FALSE,
  transform = torchvision::transform_to_tensor
)

train_loader <- dataloader(
  dataset = dataset_subset(train_data, 1:num_samples),
  batch_size = 256,
  shuffle = TRUE
)

test_loader <- dataloader(
  dataset = dataset_subset(test_data, 1:num_samples),
  batch_size = 256,
  shuffle = FALSE
)
```

### Construction d'un block _residual_

Pour définir notre block _residual_, on crée une classe `torch.nn.Module` et on y définit deux méthodes : `__init__` et `forward` :
qui hérite de `torch.nn.Module`.

```python
def align(num_in, num_out, stride):
    if num_in != num_out or stride > 1:
        return nn.Sequential(
            nn.Conv2d(
                num_in, num_out, kernel_size=3, stride=stride, padding=1, bias=False
                ),
            nn.BatchNorm2d(num_out)
            )
    else:
        return lambda x: x

class ResBlock(nn.Module):
    def __init__(self, num_in, num_out, stride):
        super(ResBlock, self).__init__()
        self.align = align(num_in, num_out, stride)
        self.conv1 = nn.Conv2d(num_in, num_out, kernel_size=3,
                            stride=stride, padding=1, bias=False)
        self.bn1 = nn.BatchNorm2d(num_out)
        self.relu = nn.ReLU(inplace=True)
        self.conv2 = nn.Conv2d(num_out, num_out, kernel_size=3,
                            stride=1, padding=1, bias=False)
        self.bn2 = nn.BatchNorm2d(num_out)

    def forward(self, x):
        o = self.conv1(x)
        o = self.bn1(o)
        o = self.relu(o)
        o = self.conv2(o)
        o = self.bn2(o)
        o = o + self.align(x)
        o = self.relu(o)
        return o
```

De la même manière, on peut créer un objet `nn_module` en `R` et y spécifier `init` et `forward` :
```{r}
align <- function(num_in, num_out, stride) {
  if (num_in != num_out || stride > 1) {
    return(nn_sequential(
      nn_conv2d(
        num_in, num_out,
        kernel_size = 3, stride = stride, padding = 1,
        bias = FALSE
      ),
      nn_batch_norm2d(num_out)
    ))
  } else {
    return(function(x) x)
  }
}

res_block <- nn_module(
  initialize = function(num_in, num_out, stride) {
    self$align <- align(num_in, num_out, stride)
    self$conv1 <- nn_conv2d(num_in, num_out,
      kernel_size = 3,
      stride = stride, padding = 1, bias = FALSE
    )
    self$bn1 <- nn_batch_norm2d(num_out)
    self$relu <- nn_relu(inplace = TRUE)
    self$conv2 <- nn_conv2d(num_out, num_out,
      kernel_size = 3,
      stride = 1, padding = 1, bias = FALSE
    )
    self$bn2 <- nn_batch_norm2d(num_out)
  },
  forward = function(x) {
    o <- self$conv1(x)
    o <- self$bn1(o)
    o <- self$relu(o)
    o <- self$conv2(o)
    o <- self$bn2(o)
    o <- o + self$align(x)
    o <- self$relu(o)
    return(o)
  }
)
```

### Constructeur _ResNet_

Pour construire notre _ResNet_, on veut créer des block _residuals_ en chaîne. En `python`, on le fait de la manière suivante (toujours en utilisant `torch.nn.Module`) :

```python
def buildResBlocks(num_in, num_out, stride, num_blocks):
    blocks = [ResBlock(num_in, num_out, stride)]
    for _ in range(1, num_blocks):
        blocks.append(ResBlock(num_out, num_out, 1))
    return nn.Sequential(*blocks)

class ResNet(nn.Module):
    def __init__(self, num_classes):
        super(ResNet, self).__init__()
        self.blocks0 = nn.Sequential(
            nn.Conv2d(
                3, 16, kernel_size=3,
                stride=1, padding=1, bias=False
                ),
            nn.BatchNorm2d(16),
            nn.ReLU(inplace=True)
            )
        self.blocks1 = buildResBlocks(16, 16, 1, 2)
        self.blocks2 = buildResBlocks(16, 32, 2, 2)
        self.blocks3 = buildResBlocks(32, 64, 2, 2)
        self.avgpool = nn.AvgPool2d(8)
        self.fc = nn.Linear(64, num_classes)

    def forward(self, x):
        n = x.shape[0]
        o = self.blocks0(x)
        o = self.blocks1(o)
        o = self.blocks2(o)
        o = self.blocks3(o)
        o = self.avgpool(o)
        o = self.fc(o.reshape(n, -1))
        return o
```

En `R`, cela donne :

```{r}
build_res_blocks <- function(num_in, num_out, stride, num_blocks) {
  blocks <- list(res_block(num_in, num_out, stride))
  for (i in 2:num_blocks) {
    blocks[[i]] <- res_block(num_out, num_out, 1)
  }
  return(do.call(nn_sequential, blocks))
}

res_net <- nn_module(
  initialize = function(num_classes) {
    self$blocks0 <- nn_sequential(
      nn_conv2d(
        3,
        16,
        kernel_size = 3,
        stride = 1,
        padding = 1,
        bias = FALSE
      ),
      nn_batch_norm2d(16),
      nn_relu(inplace = TRUE)
    )
    self$blocks1 <- build_res_blocks(16, 16, 1, 2)
    self$blocks2 <- build_res_blocks(16, 32, 2, 2)
    self$blocks3 <- build_res_blocks(32, 64, 2, 2)
    self$avgpool <- nn_avg_pool2d(kernel_size = 8)
    self$fc <- nn_linear(64, num_classes)
  },
  forward = function(x) {
    n <- dim(x)[1]
    o <- self$blocks0(x)
    o <- self$blocks1(o)
    o <- self$blocks2(o)
    o <- self$blocks3(o)
    o <- self$avgpool(o)
    o <- torch_flatten(o, start_dim = 2)
    o <- self$fc(o)
    return(o)
  }
)
```

### Instanciation modèle et optimiseurs

Partie un peu plus rapide et simple, quasi identique dans les deux cas :

```python
device = "cpu"
model = ResNet(10).to(device)
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)
```

```{r}
device = "cpu"
model <- res_net$new(10)$to(device = device)
model
```

```{r}
optimizer <- optim_adam(model$parameters, lr = 0.001)
optimizer
```

### Apprentissage

```python
def train():
    for epoch in range(1, 11):
        for i, (x, y) in enumerate(trainLoader):
            (x, y) = x.to(device), y.to(device)
            o = model(x)
            loss = F.cross_entropy(o, y)

            optimizer.zero_grad()
            loss.backward()
            optimizer.step()
            if i % 100 == 0:
                 print("Epoch: {}\tLoss: {}".format(epoch, loss.item()))
```

```{r}
train <- function() {
  for (epoch in 1:10) {
    i <- 0
    coro::loop(for (batch in train_loader) {
      x <- batch[[1]]
      y <- batch[[2]]
      o <- model(x)
      loss <- nnf_cross_entropy(o, y)

      optimizer$zero_grad()
      loss$backward()
      optimizer$step()
      if (i %% 100 == 0) {
        cat(sprintf("Epoch: %d\tLoss: %.4f\n", epoch, loss$item()))
      }
      i <- i + 1
    })
  }
}
```

### Comparaison temps de calcul

```python
tic = time.time()
train()
print(time.time()-tic, "s")
```

```{r}
tic <- system.time(
  train()
)
tic
```

Sur mon ordinateur :
- `python` : 27.6s (elapsed)
- `R` : 33.1s

Des temps de calcul très comparables !

### Précision

Pour tester la précision :

```python
n, N = 0, 0
with torch.no_grad():
    for (x, y) in testLoader:
        (x, y) = x.to(device), y.to(device)
        o = model(x)
        _, ŷ = torch.max(o, 1)
        N += y.size(0)
        n += torch.sum(ŷ == y).item()
    print("Accuracy: {}".format(n/N))
```

```{r}
with_no_grad({
  n_tests_ok <- 0
  n_tests <- 0
  coro::loop(for (batch in test_loader) {
    x <- batch[[1]]
    y <- batch[[2]]
    o <- model(x)
    yest <- torch_max(o, dim = 2)[[2]]
    n_tests <- n_tests + y$shape
    n_tests_ok <- n_tests_ok + torch_sum(y == yest)$item()
  })
  cat("Accuracy", n_tests_ok / n_tests, "\n")
})
```

Les deux codes donnent des résultats semblables !