# JDS1123 Code
###############################################################################
# Using stressor the first time we recommend referencing the stressor.pdf
#  (vignette) on how to install the correct version of python. Once python is
#  installed and added to the PATH variable. You can install stressor with:
#
#  install.packages("stressor")
#
###############################################################################
library(stressor)
library(tidyverse)

# Figure 1 Created with Tikzpictures.

# Figure 2 Also created with Tikzpictures.


# Figure 3 Clustering example with cross validation
set.seed(43421)
vec1 <- rnorm(100)
vec2 <- rnorm(100)
mat <- matrix(vec1, vec2, nrow = 100, ncol = 2)

clus <- stressor::cv_cluster(mat, 5, 5)
mat2 <- cbind(mat, clus)
matdf <- as.data.frame(mat2)
matdf$clus <- factor(matdf$clus)

# Visualization code needed for SCV figure. (Figure 3)
ggplot(matdf, aes(x = V1, y=V2, col=clus, pch = clus)) +
  geom_point(size = 3, alpha = .6) +
  scale_color_brewer(palette = "Dark2") +
  scale_x_continuous(breaks = seq(-4, 3), limits = c(-4, 3)) +
  scale_y_continuous(breaks = c(-2, -1, 0, 1, 2, 3), limits = c(-2, 3)) +
  labs(col = "Groups", pch = "Groups")


# Table 3 code Crop Suitability example
#  NOTE: If this is your first time using the stressor package. We strongly recommend looking at
#   the stressor.pdf which is the vignette for the stressor package and contains crucial steps to
#   be able to run the stressor package, especially the create_virtualenv() function.
apmc <- readRDS("apmc.RDS")
stressor::create_virtualenv()
ml_ap <- mlm_classification(AP ~ ., train_data = apmc, seed = 43421)
mlm_apmc_cv <- cv(ml_ap, data = apmc, n_folds = 10)
mlm_apmc_scv <- cv(ml_ap, data = apmc, n_folds = 10, k_mult = 5)

# The above code can be skipped with this load statement
# load("apmc_table.Rdata")
ml_cv_accuracy <- score(apmc$AP, ml_cv, metrics = "Accuracy")
ml_scv_accuracy <- score(apmc$AP, ml_scv, metrics = "Accuracy")

# Table 4 code Corn Yield example
corn_yield <- readRDS("corn_yield.RDS")
stressor::create_virtualenv()
mlm_yield <- mlm_regressor(YIELD ~ ., corn_yield, seed = 43421)
mlm_yield_cv <- cv(mlm_yield, corn_yield, n_folds = 10)
mlm_yield_scv <- cv(mlm_yield, corn_yield, n_folds = 10, k_mult = 5)
mlm_yield_scv_latlon <- cv(mlm_yield, corn_yield, n_folds = 10, k_mult = 5,
                       grouping_formula = ~ lat + lon)

# The above code can be skipped with this load statement
# load("yield.Rdata")

yield_cv_rmse <- score(corn_yield$YIELD, mlm_yield_cv, metrics = "RMSE")
yield_scv_rmse <- score(corn_yield$YIELD, mlm_yield_scv, metrics = "RMSE")
yield_scv_latlon_rmse <- score(corn_yield$YIELD, mlm_yield_scv_latlon,
                               metrics = "RMSE")

# Figure 4 code Demonstrating the creation of the data and visualization.

# Replication function
validation_rep <- function(formula, t_data, n_folds = 10, k_mult = NULL,
                           grouping_formula = NULL, times = 10) {
  rmse_list <- vector(mode = "list", length = times)
  data <- model.frame(formula = formula, data = t_data)
  vv <- attr(terms(formula(data)), which = "variables")
  rr <- as.character(vv[[2]]) # The response variable name
  for (i in  1:times) {
    print(i)
    mlm_object <- mlm_regressor(formula, data)
    if (!is.null(grouping_formula)) { # statement needed for SCV_latlon
      data <- t_data
    }
    cv_method <- cv(object = mlm_object, data = data, n_folds = n_folds,
                    k_mult = k_mult, grouping_formula = grouping_formula)
    rmse_list[[i]] <- score(data[, rr], cv_method, metrics = "RMSE")
  }
  rmse_list
}

# This next chunk of code can take up to a day and a half to complete. Included
#  in the data folder is data_joint.RDS.
set.seed(43421)
cv_20 <- validation_rep(formula, corn_yield_train, times = 20)
scv_20 <- validation_rep(formula, corn_yield_train, n_folds = 10, k_mult = 5,
                         times = 20)
scv_latlon_20 <- validation_rep(formula, corn_yield_train, n_folds = 10,
                                k_mult = 5, grouping_formula = ~ lat + lon,
                                times = 20)

data_cv <- unlist(cv_20)
models <- unlist(lapply(cv_20, row.names))

data_scv <- unlist(scv_20)
models2 <- unlist(lapply(scv_20, row.names))

data_scv_latlon <- unlist(scv_latlon_20)
models3 <- unlist(lapply(cv_20, row.names))

rmse <- c(data_cv, data_scv, data_scv_latlon)
models <- c(models, models2, models3)
Method <- c(rep(c("CV", "SCV", "SCV_latlon"), each = 360))

data_joint <- data.frame(models = models, rmse = rmse, Method = Method)
data_joint$models <- as.factor(data_joint$models)
data_joint$Method <- as.factor(data_joint$Method)

# data_joint <- readRDS("data_joint.RDS") # Uncomment if you want to load.

cblindcol <- RColorBrewer::brewer.pal(3, "Dark2")
data_joint$models <- with(data_joint, reorder(models, rmse, median))

# Visualization Code for Figure 4 Side by Side box plot comparison.
ggplot(data_joint, aes(x = models, y = rmse, fill = Method)) +
  geom_boxplot(coef = 9999) +
  scale_y_continuous(name = "RMSE", breaks = seq(15, 70, 5),
                     limits = c(15, 70)) +
  scale_fill_brewer(palette = "Dark2") +
  theme(axis.text.x = element_text(angle = -60, vjust = .6, hjust = 0),
        axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14),
        legend.position = c(0, 1),
        legend.justification = c(0, 1),
        legend.background = element_blank()) +
  xlab("Models")
