```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(reticulate)
use_virtualenv("~/.virtualenvs/r-reticulate")
#py_install("keras")
if (!require('tidyverse')) install.packages("tidyverse")            # data wrangling
if (!require('text2vec')) install.packages("text2vec")              # word embeddings
if (!require('umap')) install.packages("umap")                      # visualization
if (!require('caret')) install.packages("caret")                    # for evaluation metrics (machine learning)
if (!require('keras3')) install.packages("keras3")                  # for machine learning
if (!require('quanteda')) install.packages("quanteda")              # functions for text processing
if (!require("smotefamily")) install.packages("smotefamily")        # to handle imbalanced data -> generates synthetic data
#if (!require("torch")) install.packages("torch")
#install_torch("lantern")
```


# Word embeddings and party manifestos
## 1. Tidying text
This is an essential task for word embeddings as well since reducing the number of tokens increase computational efficiency and gives us more meaningful relationships between words. 

```{r}
df <- readRDS("./data_prep.RDS")
df <- subset(df, issue!="NA")
df$id <- 1:nrow(df)

# pre-process text: remove punctuation, numbering, superfluous whitespace, linebreaks 
df$text_prep <- tolower(df$sentence_context) %>% 
  paste0(" ", .) %>%
  str_replace_all(., "\\d+", " ") %>% 
  str_replace_all(., "[[:punct:]]", " ") %>% 
  str_replace_all(., "[[:symbol:]]", " ") %>%  
  str_replace_all(., "\n", " ") 

# tokenize corpus into unigrams
toks <- quanteda::tokens(df$text_prep) %>% tokens_remove(quanteda::stopwords(language="en"))
toks[[1]][1:5]

# further pruning of vocabulary -> only terms that appear at least 5 times 
# first dfm transformation to extract those features before removing others from tokenized object
features <- dfm(toks, verbose = FALSE) %>% dfm_trim(min_termfreq = 5) %>% featnames()

toks_sub <- tokens_select(toks, features, padding=T)


# add meta information to dfm 
docvars(toks, "text_id") <- df$id
docvars(toks, "country") <- df$country
docvars(toks, "issue") <- df$issue
docvars(toks, "parfam") <- df$parfam_n
docvars(toks, "year") <- df$year

# for embedding regression, let's add a variable that distinguishes between pre- and post 2015
df$prepost <- ifelse(df$year<2015, "pre", "post")
table(df$prepost)

docvars(toks, "prepost") <- df$prepost
```

# Supervised classification with embeddings

We are moving closer to a more accurate representation of text. In principle, we've got more information, so this should also help us in the predicting task we have started on Friday. Does it? 

Don't be afraid if you cannot follow every step here. This is meant to be a showcase session. At home, go through it line by line. Most of it will be easy to understand. 

The first part is just to ensure that we use exactly the same sample as yesterday. For your own application, this is not really important.
```{r}
embeddings <- readRDS("./embeddings_mat.RDS")

# Classification of so many categories is difficult, let's start simply. We want to distinguish statements on welfare from statements on the national way of life (immigration, patriotism, etc.)
df_sub <- subset(df, issue %in% c("504", "505", "601", "602"))

# for comparison, let's take the exact same ids as on Friday
df_sub2 <- readRDS("pred_df.RDS")  

df_sub <- subset(df_sub, id %in% df_sub2$id)
```


Now this is new: every input vector needs to be of the same size. There are two options to proceed: either we take the longest sequence as our default length or we take a length informed by summary measures (like the mean or median) as a benchmark. We decide for the latter since taking the longest sequence makes computation much slower. But be aware: this is a feature we could/should tweak if we want to improve performance.
```{r}
# count tokens
df_sub$count <- str_count(df_sub$text_prep, "\\w+")
summary(df_sub$count)

# on the basis of summary statistics, we need to decide how to prune sequences to n number of vectors -> here, we use the mean length of sequences but this is arbitrary and a potential hyperparameter
max_len <- round(median(df_sub$count,na.rm=T))
```


Now, we can create a new vocabulary (because we have a subset of our data frame by now) and prune our sequences to *max_len*.
```{r}
tokens <- word_tokenizer(df_sub$text_prep)
it <- itoken(tokens)
vocab <- create_vocabulary(it, stopwords = c(stopwords(language="en"),"also", "s", "t", "d"))
vocab <- subset(vocab, term_count>5 & term_count<20000)

# sequencing
vectorize_layer <- layer_text_vectorization(
  max_tokens = nrow(vocab),          
  output_sequence_length = max_len    
)
vectorize_layer %>% adapt(df_sub$text_prep)
features <- vectorize_layer(df_sub$text_prep)
features <- as.array(features)
head(features)
```


In machine learning, we want to train a model based on annotated data and evaluate its performance on test data. We always need to split our data into two subsets. 

**Note**: The best way is actually to split into three subsets: train, test and validation (left out) because train and test data can be mixed in the repeated training sequence.

```{r}
# create training data (let's use same IDs as on Friday)
index_train <- readRDS("./train_ids.RDS") %>% as.numeric
index_test <- readRDS("./test_ids.RDS")  %>% as.numeric

# extract features and labels (topics)
row_names_train <- as.numeric(row.names(df_sub)[with(df_sub, id %in% index_train)])
row_names_test <- as.numeric(row.names(df_sub)[with(df_sub, id %in% index_test)])


x_train <- features[row_names_train,]
x_test <- features[row_names_test,]

# define the labels of interest. Here, distinguish statements on welfare state from migration
df_sub$labels <- ifelse(df_sub$issue %in% c("504", "505"), 1, 0)
y_train <- subset(df_sub, id %in% index_train) 
y_train <- unlist(y_train[,"labels"])
y_test <- subset(df_sub, id %in% index_test) 
y_test <- unlist(y_test[,"labels"])
```


*Tuning your code*
In case you want to balance the data, there are different approaches. You could undersample majority classes, oversample minority classes or use SMOTE which creates synthetic minority classes.
Here, you'll see how SMOTE would work
```{r}
data_train <- data.frame(as.array(x_train), y_train)
smote_data <- SLS(data_train, target = data_train$y_train, K = 5, C = 5)
table(smote_data$data$y_train) # worked!

x_train_balanced <- as.matrix(smote_data$data[, 1:max_len])
y_train_balanced <- smote_data$data$y_train
```


Let's save the data for later usage (the transformer model).
```{r}
train <- data.frame(text_prepared = df_sub[df_sub$id %in% index_train,7], labels = y_train)
train$labels_text <- ifelse(train$labels==0, "migration", "welfare")
write.csv(train, "./data_transformers/training.csv")
test <- data.frame(text_prepared = df_sub[df_sub$id %in% index_test,7], labels = y_test)
test$labels_text <- ifelse(test$labels==0, "migration", "welfare")
write.csv(test, "./data_transformers/test.csv")
```


Now, we are at the point to define our model, to compile it, and to run it! 

A simple model (only using information from our embeddings)
```{r}
model <- keras_model_sequential() %>%
  layer_embedding(
    input_dim = dim(embeddings)[1],  # Assuming the number of words in your vocabulary matches the rows of t_mat
    input_length = max_len,
    output_dim = dim(embeddings)[2],  # Assuming the dimensionality of your embeddings
    weights = list(embeddings),
    trainable = F
  ) %>%
  layer_flatten() %>%
  layer_dense(units = 1, activation = "sigmoid")  

summary(model)

# Step 4: Compile your model
model %>% compile(
  loss = "binary_crossentropy",
  optimizer = optimizer_adam(learning_rate = 0.001),
  metrics = c('accuracy', metric_precision(), metric_recall())
)


stop_if_no_improvement <- callback_early_stopping(
  monitor = "val_loss",
  patience = 5,
  restore_best_weights = TRUE
)



# Step 5: Train your model
history <- model %>% fit(
  x = x_train,
  y = y_train,
  epochs = 20,
  batch_size = 128,
  validation_data = list(x_test, y_test),
  callbacks = list(stop_if_no_improvement)
)

print(paste0("Our best epoch is epoch ", (best_epoch <- which(history$metrics$val_loss == min(history$metrics$val_loss)))))
print(paste0("The smallest loss in our validation set is ", (loss <- history$metrics$val_loss[best_epoch] %>% round(3))))
print(paste0("The accuracy in our validation set is ", acc <- history$metrics$val_acc[best_epoch] %>% round(3)))
```
Okay, we see that there is quite a large error. The loss in our validation set is 0.23. Even the accuracy is not that great. 9% of our cases are wrongly predicted. That's more than it was on Friday with Naive Bayes!

If we look at the confusion matrix: 
```{r}
# predict
pred_results <- as.data.frame(predict(model, x_test))

# classify cases with probability > 0.5 as 1 
pred_results$pred <- ifelse(pred_results$V1>=0.5, 1, 0)

# combine predictions with true annotated data
pred_results$true <- y_test

confusionMatrix(as.factor(pred_results$pred), as.factor(pred_results$true))
```
That's not a good performance. Accuracy is ok but most sentences on migration were misclassified (, thus the low sensitivity).


The pre-trained embeddings treat each word separately and disregard structure sentence and inter-word relationships which might help to classify our data correctly. 

Therefore, we usually add hidden layers which can do that. Let's tune our deep learning model with these layers.
1. Dense layers provide a weight to each output.
2. Convolutional layers build n-grams of words which might actually help us to retrieve more contextual meaining. 
3. Our pre-trained embeddings can actually be fine-tuned by setting "trainable" to TRUE -> this makes them more domain specific. 
```{r}
# define model
model <- keras_model_sequential() %>%
  layer_embedding(
    input_dim = dim(embeddings)[1],  # Assuming the number of words in your vocabulary matches the rows of t_mat
    input_length = max_len,
    output_dim = dim(embeddings)[2],  # Assuming the dimensionality of your embeddings
    weights = list(embeddings),
    trainable = T
  ) %>%
  layer_conv_1d(filters = 128, kernel_size = 5, activation = 'relu') %>%
  layer_global_max_pooling_1d() %>%
  layer_dense(units = 128, activation = "relu") %>% # for dense layers between input and output, we often rely on "rectified linear units" that are non-linear and do not lead to vanishing gradient (weights become too small)
  layer_dropout(rate = 0.4) %>%  # dropout to reduce overfitting (nodes are randomly held out)
  layer_dense(units = 1, activation = "sigmoid")  

summary(model)

# Step 4: Compile your model
model %>% compile(
  loss = "binary_crossentropy",
  optimizer = optimizer_adam(learning_rate = 0.001),
  metrics = c('accuracy', metric_precision(), metric_recall())
)


stop_if_no_improvement <- callback_early_stopping(
  monitor = "val_loss",
  patience = 5,
  restore_best_weights = TRUE
)



# Step 5: Train your model
history <- model %>% fit(
  x = x_train,
  y = y_train,
  epochs = 10,
  batch_size = 128,
  validation_data = list(x_test, y_test),
  callbacks = list(stop_if_no_improvement)
)

print(paste0("Our best epoch is epoch ", (best_epoch <- which(history$metrics$val_loss == min(history$metrics$val_loss)))))
print(paste0("The smallest loss in our validation set is ", (loss <- history$metrics$val_loss[best_epoch] %>% round(3))))
print(paste0("The accuracy in our validation set is ", acc <- history$metrics$val_acc[best_epoch] %>% round(3)))
#saveRDS(model, "./model_e7b16d300.RDS")
```
Let's predict our test data with the more sophisticated deep learning model.
```{r}
# now predict probability for every speech
pred_results <- as.data.frame(predict(model, x_test))
```

...and show evaluation metrics.
```{r}
# transform probability into label: <0.5 -> 0, >0.5=1
pred_results$pred <- ifelse(pred_results$V1>=0.5, 1, 0)

# add true category, text and party
pred_results$true <- y_test
pred_results$issue_cmp <- df_sub$issue[index_test]
pred_results$sentence <- df_sub$sentence[index_test]
pred_results$party <- df_sub$partyname[index_test]

(confusion_matrix <- table(pred_results$pred, pred_results$true))

confusionMatrix(as.factor(pred_results$pred), as.factor(pred_results$true))
```
Now, the performance is way better than it was when using only quanteda's built-in functions. 

*Try it out*
Try building your own neural network. Can you beat the predictive quality?
```{r}

```




```{r}
# define model
model <- keras_model_sequential() %>%
  layer_embedding(
    input_dim = dim(embeddings)[1],  # Assuming the number of words in your vocabulary matches the rows of t_mat
    input_length = max_len,
    output_dim = dim(embeddings)[2],  # Assuming the dimensionality of your embeddings
    weights = list(embeddings),
    trainable = T
  ) %>%
  layer_conv_1d(filters = 64, kernel_size = 3, activation = 'relu') %>%
  layer_max_pooling_1d(pool_size = 2) %>%
  layer_dropout(rate = 0.2) %>%  
  layer_conv_1d(filters = 128, kernel_size = 5, activation = 'relu') %>%
  layer_max_pooling_1d(pool_size = 2) %>%
  layer_dropout(rate = 0.2) %>%  
  layer_conv_1d(filters = 256, kernel_size = 5, activation = 'relu') %>%
  layer_global_max_pooling_1d() %>%
  layer_dropout(rate = 0.2) %>%  
  layer_dense(units = 128, activation = "relu") %>% 
  layer_dropout(rate = 0.2) %>%  
  layer_dense(units = 1, activation = "sigmoid")  

summary(model)

# Step 4: Compile your model
model %>% compile(
  loss = "binary_crossentropy",
  optimizer = optimizer_adam(learning_rate = 0.001),
  metrics = c('accuracy', metric_precision(), metric_recall())
)


stop_if_no_improvement <- callback_early_stopping(
  monitor = "val_loss",
  patience = 5,
  restore_best_weights = TRUE
)



# Step 5: Train your model
history <- model %>% fit(
  x = x_train,
  y = y_train,
  epochs = 10,
  batch_size = 128,
  validation_data = list(x_test, y_test),
  callbacks = list(stop_if_no_improvement)
)

print(paste0("Our best epoch is epoch ", (best_epoch <- which(history$metrics$val_loss == min(history$metrics$val_loss)))))
print(paste0("The smallest loss in our validation set is ", (loss <- history$metrics$val_loss[best_epoch] %>% round(3))))
print(paste0("The accuracy in our validation set is ", acc <- history$metrics$val_acc[best_epoch] %>% round(3)))
#saveRDS(model, "./model_e7b16d300.RDS")
```


# Multi-class prediction 
```{r}
index_train <- readRDS("./train_ids.RDS") %>% as.numeric
index_test <- readRDS("./test_ids.RDS")  %>% as.numeric

# extract features and labels (topics)
row_names_train <- as.numeric(row.names(df_sub)[with(df_sub, id %in% index_train)])
row_names_test <- as.numeric(row.names(df_sub)[with(df_sub, id %in% index_test)])

x_train <- features[[row_names_train-1]]
x_test <- features[[row_names_test-1]]

# define the labels of interest. Here, distinguish statements on welfare state from migration
df_sub$labels[df_sub$issue=="504"] <- 0
df_sub$labels[df_sub$issue=="505"] <- 1
df_sub$labels[df_sub$issue=="601"] <- 2
df_sub$labels[df_sub$issue=="602"] <- 3

y_train <- subset(df_sub, id %in% index_train) 
y_train <- unlist(y_train[,"labels"]) 
y_test <- subset(df_sub, id %in% index_test) 
y_test <- unlist(y_test[,"labels"])

y_train <- to_categorical(y_train, num_classes = 4)
y_test <- to_categorical(y_test)

# define model
model2 <- keras_model_sequential() %>%
  layer_embedding(
    input_dim = dim(embeddings)[1],  # Assuming the number of words in your vocabulary matches the rows of t_mat
    input_length = max_len,
    output_dim = dim(embeddings)[2],  # Assuming the dimensionality of your embeddings
    weights = list(embeddings),
    trainable = T
  ) %>%
  layer_conv_1d(filters = 128, kernel_size = 5, activation = 'relu') %>%
  layer_global_max_pooling_1d() %>%
  layer_dense(units = 128, activation = "relu") %>% # for dense layers between input and output, we often rely on "rectified linear units" that are non-linear and do not lead to vanishing gradient (weights become too small)
  layer_dropout(rate = 0.4) %>%  # dropout to reduce overfitting (nodes are randomly held out)
  layer_dense(units = 4, activation = "softmax")  # this is the decisive parameter for multiclass prediction; often softmax is used as activation, the number of units must equal the number of categories

summary(model2)

# Step 4: Compile your model
model2 %>% compile(
  loss = "categorical_crossentropy", # for multiclass, we need to set a different loss function
  optimizer = optimizer_adam(learning_rate=0.0001),
  metrics = c("accuracy")
)


stop_if_no_improvement <- callback_early_stopping(
  monitor = "val_loss",
  patience = 5,
  restore_best_weights = TRUE
)


# Step 5: Train your model
history <- model2 %>% fit(
  x = x_train,
  y = y_train,
  epochs = 15,
  batch_size = 128,
  validation_data = list(x_test, y_test),
  callbacks = list(stop_if_no_improvement)
)

print(paste0("Our best epoch is epoch ", (best_epoch <- which(history$metrics$val_loss == min(history$metrics$val_loss)))))
print(paste0("The smallest loss in our validation set is ", (loss <- history$metrics$val_loss[best_epoch] %>% round(3))))
print(paste0("The accuracy in our validation set is ", acc <- history$metrics$val_acc[best_epoch] %>% round(3)))
```

```{r}
# now predict probability for every speech
pred_results2 <- as.data.frame(predict(model2, x_test))
colnames(pred_results2) <- c("504", "505", "601", "602")
pred_results2$pred <- colnames(pred_results2)[apply(pred_results2,1,which.max)]

# add true category, text and party
pred_results2$true <- df_sub$issue[df_sub$id %in% index_test]

confusionMatrix(as.factor(pred_results2$pred), as.factor(pred_results2$true))
```