Data - Performance Examples.Rmd

---
title: "Data - Performance Examples"
author: "Pablo Morala"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

####################################
# 1 - Load all needed packages and functions
####################################
library(nn2poly)
library(nn2poly.tools)
library(keras)
library(tensorflow)
library(cowplot)
library(ggplot2)
library(patchwork)

IEEE_width_inches <- 3.5
my_width <- IEEE_width_inches

my_seed <- 1
set.seed(my_seed)
tensorflow::tf$random$set_seed(my_seed) # Needed to have reproducible results with keras
```


In this Rmd file we are going to show how to use NN2Poly in 2 different situations with synthetically generated data and store the results to plot them in `Figures - Performance examples.Rmd`.


# Data generation

Data will be held the same in the following examples.  
```{r}
source("data_generation2.R")
source("plot_taylor_and_synaptic_potentials2.R")
```


Here we create the polynomial : $4x_1 - 3 x_2x_3$: 
```{r}
# Define the polynomial to be used in the desired form:
original_polynomial <- list()
original_polynomial$labels <- list(c(1), c(2,3))
original_polynomial$values <- c(4,-3)
```

Then we generate the data:

```{r}
# Define number of variables p and sample n
p <- 3
n_sample <- 500

data <- data_generation2(n_sample, p, original_polynomial, error_var = 0.1)

head(data)
```


```{r}

# Data scaling
maxs <- apply(data, 2, max) 
mins <- apply(data, 2, min) 
data <- as.data.frame(scale(data, center = mins + (maxs - mins) / 2, scale = (maxs - mins) / 2))

# Divide in train (0.75) and test (0.25)
index <- sample(1:nrow(data), round(0.75 * nrow(data)))
train <- data[index, ]
test <- data[-index, ]

train_x <- as.matrix(subset(train, select = -c(Y)))
train_y <- as.matrix(train$Y)

test_x <- as.matrix(subset(test, select = -c(Y)))
test_y <- as.matrix(test$Y)


# Scale the data in the desired interval and separate train and test
scale_method <- "-1,1"
data_scaled <- scale_data(data, scale_method)
aux <- divide_train_test(data_scaled, train_proportion = 0.75)
train <- aux$train
test <- aux$test

plot(data_scaled)
```



# Keras and nn2poly hyperparameters

These parameters will also be held constant over the networks

```{r}
# keras hyperparameters
my_loss <- "mse"
my_metrics <- "mse"
my_epochs <- 3000
my_batch <- 50
my_validation_split <- 0.2
my_verbose <- 0

# NN structure
L <- 3
h_neurons_at_each_layer <- 100
af <- "tanh"
my_optimizer <- optimizer_adam()
# h neurons with 1 at the end for regression
h_neurons_vector <- rep(h_neurons_at_each_layer,L)
h_neurons_vector <- c(h_neurons_vector,1)
# list with the af at each layer:
af_string_list <- vector(mode = "list", length = L+1)
for (i in 1:L){
  af_string_list[i] <- af
}
# Add linear at the end so we have regression setting.
af_string_list[L+1] <- "linear"

# nn2poly hyperparameter
forced_max_Q <- 3

# q taylor with 1 at the end for regression
q_taylor_at_each_layer <- 8 # We can set a high value as we are going
# to limit it with a forced max Q
q_taylor_vector <- rep(q_taylor_at_each_layer,L)
q_taylor_vector <- c(q_taylor_vector,1)
```



# Function to generate the desired graphs:

```{r}
perform_example_from_train_test<- function(train,
                                           test,
                                           af_string_list,
                                           h_neurons_vector,
                                           q_taylor_vector,
                                           forced_max_Q,
                                           my_max_norm,
                                           my_optimizer,
                                           my_loss,
                                           my_metrics,
                                           my_epochs,
                                           my_batch,
                                           my_validation_split,
                                           my_verbose,
                                           all_partitions) {

  # Obtain parameters:
  p <- ncol(train) - 1
  
  # Divide again in x and y
  train_x <- as.matrix(subset(train, select = -c(p+1)))
  train_y <- as.matrix(subset(train, select = c(p+1)))
  
  test_x <- as.matrix(subset(test, select = -c(p+1)))
  test_y <- as.matrix(subset(test, select = c(p+1)))
  
  # Build the nn
  nn <- build_keras_model(
    p,
    af_string_list,
    h_neurons_vector,
    my_max_norm
  )
  
  # Compile the model
  compile(nn,
          loss = my_loss,
          optimizer = my_optimizer,
          metrics = my_metrics
  )
  
  # Fit the model
  history <- fit(nn,
                 train_x,
                 train_y,
                 verbose = my_verbose,
                 epochs = my_epochs,
                 validation_split = my_validation_split,
                 batch_size = my_batch
  )
  plot_history <- plot(history)

  # Obtain the predicted values with the NN to compare them
  prediction_NN <- predict(nn, test_x)
  
  # plot that comparison between NN and original
  plot_NN_performance <- plot_NN_PR_comparison(unname(test_y), prediction_NN)
  
  plot_NN_performance <- plot_NN_performance +
    ggplot2::labs(x = "Original Y") +
    theme_half_open()
  
  # MSE between NN and Poly
  n_test <- length(test_y)
  MSE_NN_vs_original<- sum((prediction_NN - test_y)^2) / n_test
  
  # Extract the weights:
  keras_weights <- keras::get_weights(nn)
  
  n <- length(keras_weights)
  
  if(my_max_norm[[1]]=="no_constraint"){
    
    n2 <- n/2
    
    nn_weights <- vector(mode = "list", length = n2)
    for (i in 1:n2){
      nn_weights[[i]] <- rbind(keras_weights[[2*i]], keras_weights[[2*i-1]])
    }
    
  } else {
    nn_weights <- keras_weights[1:(n - 2)]
    
    nn_weights[[n - 1]] <- rbind(keras_weights[[n]], keras_weights[[n - 1]])
  }
  
  # use the  nn2poly algorithm
  all_layers_coeffs <- nn2poly::nn2poly_algorithm(
    weights_list = nn_weights,
    af_string_list = af_string_list,
    q_taylor_vector = q_taylor_vector,
    forced_max_Q = forced_max_Q
  )
  
  final_poly <- all_layers_coeffs[[length(all_layers_coeffs)]]
  
  
  prediction_poly <- as.vector(eval_poly(x = test_x, poly = final_poly))
    
  # MSE between NN and Poly
  MSE_NN_vs_poly <- sum((prediction_NN - prediction_poly)^2) / n_test
  
  # Plot the PR vs NN predictions
  plot_PR_vs_NN <- plot_NN_PR_comparison(prediction_poly, prediction_NN) +
    theme_half_open()
  
  # Plot the taylor expansion at each layer:
  plot_taylor <- plot_taylor_and_synpatic_potentials2(
    data = train,
    weights_list = nn_weights,
    af_string_list = af_string_list,
    q_taylor_vector = q_taylor_vector,
    forced_max_Q = forced_max_Q,
    my_max_norm
  )
  
  output <- vector(mode = "list", length = 0)
  output$MSE_NN_vs_poly <- MSE_NN_vs_poly
  output$MSE_NN_vs_original <- MSE_NN_vs_original
  output$plot_NN_performance <- plot_NN_performance
  output$plot_PR_vs_NN <- plot_PR_vs_NN
  output$plot_taylor <- plot_taylor
  output$plot_history <- plot_history
  output$final_poly <- final_poly
  output$train_x <- train_x
  
  
  return(output)
}

```




# Example 1


Tanh AF, l1-norm constraint, 3 hidden layers, 50 neurons per layer.

```{r, warning=FALSE}

my_max_norm <- list("l1_norm", 1)

example1 <- perform_example_from_train_test(
          train = train,
          test = test,
          af_string_list = af_string_list,
          h_neurons_vector = h_neurons_vector,
          q_taylor_vector = q_taylor_vector,
          forced_max_Q = forced_max_Q,
          my_max_norm = my_max_norm,
          my_optimizer = my_optimizer,
          my_loss = my_loss,
          my_metrics = my_metrics,
          my_epochs = my_epochs,
          my_batch = my_batch,
          my_validation_split = my_validation_split,
          my_verbose
        )
```


```{r}
example1
```



```{r}
# Remove history  plot as it is not used in the Figures
example1$plot_history <- NULL
saveRDS(example1,"temporal/nn_performance_example_l1")
```



# Example 2


Same as before, no constraints

```{r, warning=FALSE}

my_max_norm <- list("no_constraint", 1)


example2 <- perform_example_from_train_test(
          train = train,
          test = test,
          af_string_list = af_string_list,
          h_neurons_vector = h_neurons_vector,
          q_taylor_vector = q_taylor_vector,
          forced_max_Q = forced_max_Q,
          my_max_norm = my_max_norm,
          my_optimizer = my_optimizer,
          my_loss = my_loss,
          my_metrics = my_metrics,
          my_epochs = my_epochs,
          my_batch = my_batch,
          my_validation_split = my_validation_split,
          my_verbose
        )
```


```{r}
example2
```


```{r}
# Remove history  plot as it is not used in the Figures
example2$plot_history <- NULL
saveRDS(example2,"temporal/nn_performance_example_nc")
```



# Repeating example 1, 10 diferent times to have average coefficients:


```{r}
my_max_norm <- list("l1_norm", 1)
my_epochs <- 3000

n_sims <- 10
simulations <- vector(mode="list", length = n_sims)
simulations_plot <- vector(mode="list", length = n_sims)

for (sims in 1:n_sims){
  print(sims)
  result <- perform_example_from_train_test(
          train = train,
          test = test,
          af_string_list = af_string_list,
          h_neurons_vector = h_neurons_vector,
          q_taylor_vector = q_taylor_vector,
          forced_max_Q = forced_max_Q,
          my_max_norm = my_max_norm,
          my_optimizer = my_optimizer,
          my_loss = my_loss,
          my_metrics = my_metrics,
          my_epochs = my_epochs,
          my_batch = my_batch,
          my_validation_split = my_validation_split,
          my_verbose
        )
  
  simulations[[sims]] <- result$final_poly
  simulations_plot[[sims]] <- result$plot_PR_vs_NN
}



```



```{r}
# Remove history  plot as it is not used in the Figures
simulations
saveRDS(simulations,"temporal/nn_performance_example_10_repetitions")
```