library("torch") # v0.13.0
library("torchvision") #v0.6.0
library("luz") #v0.4.0 makes it "easier" to use "torch"
library("magick") #2.8.5. 
## Linking to ImageMagick 6.9.11.60
## Enabled features: fontconfig, freetype, fftw, heic, lcms, pango, webp, x11
## Disabled features: cairo, ghostscript, raw, rsvg
## Using 8 threads
# Custom functions for viewing images and some of their details:
# view_random_image() and view_random_torch_sample()
source("utils.R")

NOTE: if your PC has CUDA compatible GPU, then most likely torch has been installed in a way

that it will run with GPU. Otherwise, it will run with CPU (slower). On Mac OS, if you later get a warning message about

MPS device, follow the instructions (involves setting a variable in .Renviron file)

1 Dataset

Here, we will use a dataset of 2000 image: arrows and vases (“pil” and “kar”)

# For running in Rstudio
view_random_image(data_dir="../torch_data/images/")
view_random_image(data_dir="../torch_data/images/", colour=TRUE)

# For knitting to Rmarkdown HTML document
# knitr::include_graphics(view_random_image(data_dir="../torch_data/images/", rmarkdown = TRUE))
# knitr::include_graphics(view_random_image(data_dir="../torch_data/images", rmarkdown = TRUE, colour=TRUE))

1.1 Metadata file

Columns used: “Gjenstand” and “Filnavn”

For coding simplicity in R, have ensured that each existing JPG image

is present in the metadata.csv file

metadata <- read.csv("../torch_data/metadata.csv")

print(head(metadata))
##     C.Nr Gjenstand Aksesjonsnr                    Filnavn Fotoark_Id
## 1 C36679       kar                         cf21076_06.jpg       3580
## 2 C27765       kar                         cf21046_03.jpg       3588
## 3 C21303       kar                         cf21706_01.jpg       3614
## 4 C21945       kar                         cf21417_11.jpg       3730
## 5 C16688       kar             cf03093_c15280_c16688a.jpg       3737
## 6 C24127       kar                         cf21417_17.jpg       3840
##       Motivbeskrivelse
## 1              Leirkar
## 2              Leirkar
## 3              Leirkar
## 4   Enkelt-gjenstander
## 5 To  leirkar med hank
## 6   Enkelt-gjenstander
# View(metadata)

Define torch::dataset with necessary preprocessed images and target labels (an R6 object)

create_musit_dataset <- torch::dataset(
  name = "musit_dataset",
  initialize = function(folder_path, metadata_path, colour=FALSE) {
    self$folder_path <- folder_path
    self$filenames <- list.files(folder_path, pattern = "\\.jpg$", full.names = TRUE)
    
    # Read the metadata CSV file
    self$metadata <- read.csv(metadata_path)
    
    # Debug: Check initial metadata and filenames
    cat("Metadata rows: ", nrow(self$metadata), "\n")
    cat("Filenames count: ", length(self$filenames), "\n")
    
    # Convert the string labels to factors, then to integers
    self$labels <- as.factor(self$metadata$Gjenstand)
    self$label_levels <- levels(self$labels)
    self$file_to_int_label <- setNames(as.integer(self$labels), self$metadata$Filnavn)
    
    # Different type of preprocessing if colour or not
    self$colour <- colour
    
  },
  
  .getitem = function(index) {
    img_path <- self$filenames[index]
    img <- magick::image_read(img_path)
    img <- torchvision::transform_resize(img, size = c(128, 128))
    
    if (!self$colour){
      channels = "gray"
      img <- magick::image_convert(img, colorspace = channels)
    } else {
      channels = "rgb"
    }
    
    # Convert to numeric array
    # NOTE: normally you would normalize (e.g. divide by 255), but
    # magick::image_data already does this
    img_array <- as.numeric(magick::image_data(img, channels = channels))
    
    # Reshape to (channels, height, width) from (width, channels, height)
    # as expected by torch.
    img_tensor <- torch_tensor(img_array)$permute(c(3, 1, 2))

    # Get the corresponding integer label for the image
    int_label <- self$file_to_int_label[basename(img_path)]
    
    # convert 1-2 to 0-1 (necessary for binary classification)
    y <- torch_tensor(as.integer(int_label) - 1)
    
    return(list(img = img_tensor, y = y))
  },
  
  .length = function() {
    length(self$filenames)
  },
  
  print_class_distribution = function() {
    print(table(self$labels))
  }
)

# For use in rest of script (several areas)
use_colour = TRUE


# Create instance of defined torch::dataset()
dataset <- create_musit_dataset(folder_path = "../torch_data/images/", 
                                metadata_path = "../torch_data/metadata.csv", colour=use_colour)
## Metadata rows:  2065 
## Filenames count:  2065

In this HTML document, the model uses colour



why we do not normalize in this specific context

view_random_image(data_dir="../torch_data/images/", plot_hist = TRUE)

view_random_image(data_dir="../torch_data/images/", colour=TRUE, plot_hist = TRUE)

viewing images directly from our torch defined dataset

dataset_grey <- create_musit_dataset("../torch_data/images/", "../torch_data/metadata.csv", colour=FALSE)
## Metadata rows:  2065 
## Filenames count:  2065
view_random_torch_sample(dataset=dataset_grey, label_levels=dataset_grey$label_levels)

dataset_colour <- create_musit_dataset("../torch_data/images/", "../torch_data/metadata.csv", colour=TRUE)
## Metadata rows:  2065 
## Filenames count:  2065
view_random_torch_sample(dataset=dataset_colour, label_levels=dataset_colour$label_levels,
                         colour=TRUE)

Total size

num_samples <- length(dataset)
print(num_samples)
## [1] 2065

Print class distribution

dataset$print_class_distribution()
## 
##  kar  pil 
##  775 1290

2 Splitting the dataset

train_ratio <- 0.7
valid_ratio <- 0.15
test_ratio <- 0.15

# Calculate number of samples per set
num_train <- floor(train_ratio * num_samples)
num_valid <- floor(valid_ratio * num_samples)
# Ensure the sum matches total samples
num_test <- num_samples - (num_train + num_valid)  

# Sanity check
cat(sprintf("Total: %d, Train: %d, Valid: %d, Test: %d\n", num_samples, num_train, num_valid, num_test))
## Total: 2065, Train: 1445, Valid: 309, Test: 311

Generate shuffled indices for the dataset

set.seed(42)  # For reproducibility, within same software+OS setup
indices <- sample(1:num_samples)

train_indices <- indices[1:num_train]
valid_indices <- indices[(num_train + 1):(num_train + num_valid)]
test_indices <- indices[(num_train + num_valid + 1):num_samples]

# Ensure no index is out of range
cat("Max index in train indices: ", max(train_indices), "\n")
## Max index in train indices:  2065
cat("Max index in valid indices: ", max(valid_indices), "\n")
## Max index in valid indices:  2062
cat("Max index in test indices: ", max(test_indices), "\n")
## Max index in test indices:  2061

2.1 Create Subsets

train_dataset <- dataset_subset(dataset, indices = train_indices)
valid_dataset <- dataset_subset(dataset, indices = valid_indices)
test_dataset <- dataset_subset(dataset, indices = test_indices)

2.2 Create DataLoaders

train_dl <- dataloader(train_dataset, batch_size = 32, shuffle = TRUE) # shuffled to reduce bias and variance during training
valid_dl <- dataloader(valid_dataset, batch_size = 32, shuffle = FALSE) #not shuffled to help consistency of evaluating performance
test_dl <- dataloader(test_dataset, batch_size = 32, shuffle = FALSE)
# Inspect number of batches
cat("Number of batches in training set: ", length(train_dl), "\n")
## Number of batches in training set:  46
cat("Number of batches in validation set: ", length(valid_dl), "\n")
## Number of batches in validation set:  10
cat("Number of batches in test set: ", length(test_dl), "\n")
## Number of batches in test set:  10

3 Neural network

3.1 Define model architecture

Note: modified model from MNIST



image/svg+xmlConvolution Pooling Convolution Fully Connected Pooling Convolution

Image taken from: https://commons.wikimedia.org/wiki/File:ConvolutionAndPooling.svg Andreas Maier, CC BY 4.0 https://creativecommons.org/licenses/by/4.0, via Wikimedia Commons

net <- nn_module(
  # Desired name
  "CNN-classification",
  
  # In the "initialize" method, you set up the architecture of the neural network
  # (i.e. define layers and parameters of the model)
  
  initialize = function(colour = FALSE, dropout_rate1 = 0.25, dropout_rate2 = 0.5) {

    num_channels <- ifelse(!colour, 1, 3)

    # 1st Convolutional Layer:
    self$conv1 <- nn_conv2d(num_channels, 32, 3, 1)  

    # 2nd Convolutional layer:
    self$conv2 <- nn_conv2d(32, 64, 3, 1) 

    # Randomly dropping features to prevent overfitting
    self$dropout1 <- nn_dropout2d(dropout_rate1)
    self$dropout2 <- nn_dropout(dropout_rate2)

    
    # Fully connected Layers:
    # (see max pooling / flatteing in forward pass)
    self$fc1 <- nn_linear(64 * 62 * 62, 128)
    
    
    # Output layer
    self$fc2 <- nn_linear(128, 1)  
  },
  
  # In the "forward" method, essentially, you define how the data is process through
  # the model defined in the "initialize" method.
  forward = function(x) {
    x %>%                                  
      self$conv1() %>%                     
      nnf_relu() %>%                       # Apply ReLU activation (replaces negative values with 0)
      self$conv2() %>%                     
      nnf_relu() %>%                       
      nnf_max_pool2d(2) %>%                
      self$dropout1() %>%                  
      torch_flatten(start_dim = 2) %>%     
      self$fc1() %>%                       
      nnf_relu() %>%                       # Apply ReLU activation
      self$dropout2() %>%                  # Apply regular dropout (after flattening)
      self$fc2()                           
  }
)

3.2 Optional: using torch to help identify output size of each layer

1st conv layer

num_channels <- ifelse(!use_colour, 1, 3)
x <- torch_randn(c(1, num_channels, 128, 128))
conv1 <- nn_conv2d(num_channels, 32, 3, 1)  
x1 <- conv1(x)
x1$size()
## [1]   1  32 126 126

Here, we can see that the number of output channels is 32

ReLU doesn’t change dimensions, but let’s apply it to be thorough

x1 <- nnf_relu(x1)

Second conv layer

conv2 <- nn_conv2d(x1$size()[2], 64, 3, 1)
x2 <- conv2(x1)
print("After conv2:")
## [1] "After conv2:"
print(x2$size())  # Should be [1, 64, 124, 124]
## [1]   1  64 124 124

ReLU again (size does not change)

x2 <- nnf_relu(x2)

Max pooling

x3 <- nnf_max_pool2d(x2, 2)
print("After max_pool2d:")
## [1] "After max_pool2d:"
print(x3$size())  # Should be [1, 64, 62, 62]
## [1]  1 64 62 62

After flattening

x4 <- torch_flatten(x3, start_dim = 2)
print("After flattening:")
## [1] "After flattening:"
print(x4$size())  # Should be [1, 64 * 62 * 62]
## [1]      1 246016
print(64 * 64* 62)
## [1] 253952

First fully connected layer

fc1 <- nn_linear(64 * 62 * 62, 128)
x5 <- fc1(x4)
print("After fc1:")
## [1] "After fc1:"
print(x5$size())  # Should be [1, 128]
## [1]   1 128

Final fully connected layer

fc2 <- nn_linear(128, 1)
x6 <- fc2(x5)
print("After fc2:")
## [1] "After fc2:"
print(x6$size())  # Should be [1, 1]
## [1] 1 1

3.3 Model fitting

num_training_rounds = 11
exists_model <- file.exists("musit-binary-colour-epochs-10.pt")
if(!exists_model){
  fitted <- net %>%
  luz::setup(
    
    # Choose loss function appropriate to the task (e.g. binary classification).
    loss = nnf_binary_cross_entropy_with_logits,
    # If multiclass, loss = nnf_cross_entropy,
    
    # How to update weights of network during training
    optimizer = optim_adam,
    metrics = list(
      # if multiclass, luz_metric_accuracy()
      luz_metric_binary_accuracy_with_logits() 
    )
  ) %>%
  # Set your `net` params in set_hparams
  set_hparams(colour = use_colour) %>% 
  # set_opt_hparams(lr = 0.0001) %>%
  
  # Specify additional training parameters
  fit(train_dl,
      epochs = num_training_rounds, 
      # Also common to evaluate on the test data during model fitting
      valid_data = valid_dl,
      callbacks = list(
        luz_callback_early_stopping(patience = 5)
      )
  )
}

NOTE: no built in functions to plot training curves (i.e. loss)

Can extract fitted$metrics$<valid or train>[[<epoch number>]]$loss

Save model 

type <- ifelse(use_colour, "colour", "greyscale")
# luz_save(fitted, paste0("musit-binary-", type, "-epochs-", num_training_rounds, ".pt"))
fitted <- luz_load(paste0("musit-binary-", type, "-epochs-10.pt"))

if (exists_model)
  print("The number of epochs in the loaded model is 10")
## [1] "The number of epochs in the loaded model is 10"

Note: if you use the shared, predefined model, there is a chance that some of the test set (and validation set)

will overlap with the training set. This is not ideal as it will inflate the model evaluation results,

but is OK for the purposes of this workshop (you can, on your own time, run a model with 10 or more epochs,

save it, and use that.

4 Model evaluation

4.1 Predict categories in test set

Note: the predictions are with respect to the “second” level. In our case, pil (“arrow”).

Can verify the “levels” with our custom function dataset$print_class_distribution()

preds <- predict(fitted, test_dl)

Convert logits to probabilities using the sigmoid function

probs <- torch_sigmoid(preds)
# probs <- preds %>% nnf_softmax(dim = -1) # or `dim = 2` # Multiclass

convert to traditional array/matrix

preds <- probs$to(device = "cpu")$squeeze() %>% as.numeric()
# preds <- probs %>% torch_argmax(dim = -1) %>% as.numeric() # multiclass

Reinitialize test “dataloader” to be the whole dataset

test_dl <- dataloader(test_dataset, batch_size = length(test_indices))

targets <- (test_dl %>% dataloader_make_iter() %>% dataloader_next())$y %>% as.numeric()

# Alternatively, for very large test set, can loop over original test_dl,
# get the targets, and perform predictions for each batch, before collecting
# final result

Apply desired threshold

threshold = 0.5 # can adjust as needed (will affect model performance metrics)
predicted_class = ifelse(preds >= threshold, 1, 0)
# predicted_class=max.col(as.matrix(probs)) # multiclass

df_res <- data.frame(predicted_prob = preds, predicted_class = predicted_class, targets = targets)

Additional, for later plotting and tabulation

df_res$correct <- ifelse(df_res$predicted_class == df_res$targets, TRUE, FALSE)
df_res$Filnavn <- dataset$metadata$Filnavn[test_indices]

# Get original class labels
original_labels <- dataset$label_levels

Show random result

NOTE: here, need to specify colour parameter used in creating torch::dataset() (create_musit_dataset)

view_random_torch_sample(dataset=test_dataset, colour=use_colour, is_test_set=TRUE, 
                         label_levels = original_labels,
                         df_res = df_res)

4.2 Model performance

4.2.1 Confusion matrix

first parameter of table becomes the ROWS

i.e. columns are the predictions, rows are the observations

conf_mat <- table(df_res$targets, df_res$predicted_class)

Use these to rename the confusion matrix

colnames(conf_mat) <- paste0(original_labels)
rownames(conf_mat) <- paste0(original_labels)

conf_mat
##      
##       kar pil
##   kar  92  19
##   pil   9 191

Number of wrong results :

sum(conf_mat) - sum(diag(conf_mat))
## [1] 28
# or 
# nrow(df_res) - sum(df_res$correct)

4.2.2 Accuracy

sum(diag(conf_mat))/sum(conf_mat)
## [1] 0.9099678

4.2.3 Precision

diag(conf_mat)/colSums(conf_mat)
##       kar       pil 
## 0.9108911 0.9095238

4.2.4 Recall

diag(conf_mat)/rowSums(conf_mat)
##       kar       pil 
## 0.8288288 0.9550000

Show random incorrect result

# NOTE: here, need to specify colour parameter used in creating torch::dataset() (create_musit_dataset())
view_random_torch_sample(dataset=test_dataset, colour=use_colour, is_test_set=TRUE, 
                         label_levels = original_labels,
                         df_res = df_res, show_wrong_result = TRUE)

5 Evaluation form

Or click on the link:

https://nettskjema.no/a/464099