# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .r
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: R
#     language: R
#     name: ir
# ---

data <- read.csv("../data/epfl-bio-322a/training_data.csv")
test_set <- read.csv("../data/epfl-bio-322a/test_data.csv")

data <- data[, -3]
dim(data)
data <- na.omit(data)
dim(data)
data.response <- as.integer(data$VALENCE.PLEASANTNESS)
data.predictors <- data[, -2]
data.predictors$Intensity <- as.numeric(data.predictors$Intensity)

set.seed(127)
i <- sample(nrow(data), 3*nrow(data)/4)
datatrain <- as.matrix(data.predictors[i,])
datatest <- as.matrix(data.predictors[-i,])
ytrain <- data.response[i]
ytest <- data.response[-i]

dim(datatrain)
length(ytrain)

#data2 <- data[, 3:ncol(data)]
data2 <- data.predictors[, apply(data.predictors, 2, function(v) var(v)!=0)]

dim(data2)

library(caret)
#library(mlr)
#library(dplyr)
#library(tidyselect)
df2 = cor(data2)
hc = findCorrelation(df2, cutoff=0.99) 
hc = sort(hc)
names <- colnames(df2[, hc])
reduced_data <- data2[,!(colnames(data2) %in% names)]
#reduced_data <- select(datascale, -all_of(names)) #datascale[,-c(hc)]
dim(reduced_data)

head(test_set)

testset<-test_set[, 2:ncol(test_set)]
testset$Intensity<-as.numeric(testset$Intensity)
testset<-testset[,(colnames(testset) %in% names(reduced_data))]
head(testset)

set.seed(135)
i <- sample(nrow(reduced_data), 3*nrow(reduced_data)/4)
aftertrain <- reduced_data[i,]
aftertest <- reduced_data[-i,]
afterYtrain <- data.response[i]
afterYtest <- data.response[-i]

dim(aftertrain)
length(afterYtrain)

afterYtrain

trainCV<-as.data.frame(cbind(afterYtrain, aftertrain))

library(tree)
toytree <- tree(afterYtrain ~ .,trainCV, minsize = 1)
plot(toytree)
text(toytree, pretty=1, cex=0.5)

library(xgboost)

# Load the library 
library(randomForest) 

# Let’s create a baseline for comparison by using the recommend defaults for each parameter and mtry=floor(sqrt(ncol(x))) or mtry=7 and ntree=500.

# Create model with default paramters
control <- trainControl(method="repeatedcv", number=10, repeats=3)
seed <- 7
metric <- "RMSE"
set.seed(seed)
mtry <- ncol(aftertrain)/3
tunegrid <- expand.grid(.mtry=mtry)
rf_default <- train(afterYtrain~., data=trainCV, method="rf", metric=metric, tuneGrid=tunegrid, trControl=control)
print(rf_default)

# Random Search
seed<-7
metric<-"RMSE"
control <- trainControl(method="repeatedcv", number=10, repeats=3, search="random")
set.seed(seed)
mtry <- ncol(aftertrain)/3
rf_random <- train(afterYtrain~., data=trainCV, method="rf", metric=metric, tuneLength=15, trControl=control)
print(rf_random)
plot(rf_random)

# Algorithm Tune (tuneRF)
set.seed(seed)
bestmtry <- tuneRF(aftertrain, afterYtrain, stepFactor=1.5, improve=1e-5, ntree=500)
print(bestmtry)

# train model
control <- trainControl(method="repeatedcv", number=10, repeats=3)
tunegrid <- expand.grid(.mtry=c(1:15), .ntree=c(1000, 1500, 2000, 2500))
set.seed(seed)
custom <- train(Class~., data=dataset, method=customRF, metric=metric, tuneGrid=tunegrid, trControl=control)
summary(custom)
plot(custom)

tuneRF(reduced_data, data.response, ntreeTry=500, stepFactor=2, improve=0.05,trace=TRUE, plot=TRUE, doBest=TRUE)

# Create random forest for regression 
pleasantness.rf <- randomForest(data.response ~ ., data = reduced_data, mtry = 1218,ntree=1000,importance = TRUE, na.action = na.omit) 

# Print regression model 
print(pleasantness.rf) 

# +
# Output to be present as PNG file  
png(file = "randomForestRegression.png") 
  
# Plot the error vs the number of trees graph 
plot(pleasantness.rf) 
  
# Saving the file 
dev.off() 
# -

which.min(pleasantness.rf$mse)

# We can see that we need at least 927 trees to have the lowest error. We will now rerun the model and add an argument called “ntrees” to indicating the number of trees we want to generate.

# +
#set.seed(123)
#pleasantness.rf.883<-randomForest(afterYtrain~.,data = aftertrain,ntree=883)
#pleasantness.rf.883
# -

varImpPlot(pleasantness.rf)

importance(pleasantness.rf)

pred.rf <- data.frame(predict(pleasantness.rf, testset, predict.all = T))
sapply((pred.rf - data.response)^2, mean, na.rm = TRUE) #mean((pred.rf - afterYtest)^2, na.rm=TRUE)

pred.rf

write.csv(data.frame(Id=1:length(pred.rf$aggregate), VALENCE.PLEASANTNESS=pred.rf$aggregate), "predictions_random_forest.csv", row.names=F)

# +
set.seed(123)

ctrl_parameters <- trainControl(method = "cv", 10, savePredictions=TRUE)
parametersGrid <-  expand.grid(max_depth=seq(1,10,2),eta = seq(0.0001,0.01,0.002),
                            nrounds=seq(9000,10000,50),colsample_bytree=c(0.5,0.7),gamma=0, min_child_weight=2, subsample=0.5)
parametersGrid
# -

modelxgboost <- train(afterYtrain~., 
                  data = trainCV,
                  method = "xgbTree",
                  trControl = ctrl_parameters,
                  tuneGrid=parametersGrid, metric="RMSE")
modelxgboost

boostAvant <- xgboost(as.matrix(aftertrain), label = afterYtrain,
                 eta = 0.0001,
                 max.depth = 2,
                 nrounds = 500,
                 objective = "reg:squarederror")

responses_predicted <- predict(boostAvant, as.matrix(aftertest))
responses_predicted

xgb.plot.importance(xgb.importance(model = boostAvant))

write.csv(data.frame(Id = 1:length(responses_predicted), VALENCE.PLEASANTNESS = responses_predicted), 'boostedBIS.csv', row.names = F)

# +
# play around with settings of xgboost - eXtreme Gradient Boosting (Tree) library
# https://github.com/tqchen/xgboost/wiki/Parameters
# max.depth - maximum depth of the tree
# nrounds - the max number of iterations
# -

library(Metrics)
smallestError <- 100
i<-0
for (depth in seq(1,10,1)) {
        for (rounds in seq(1,20,1)) {
            for(eta in seq(0.0001,0.01,0.01)) {
                print(i)
                i<-i+1
                # train
                bst <- xgboost(data = as.matrix(aftertrain),
                               label = afterYtrain, neta=eta,
                               max.depth=depth, nround=rounds,
                               objective = "reg:squarederror", verbose=0)
                gc()
               
                # predict
                predictions <- predict(bst, as.matrix(aftertest), outputmargin=TRUE)
                aftertest<-unlist(aftertest, recursive = TRUE, use.names = TRUE)
                err <-  rmse(as.numeric(aftertest), as.numeric(predictions))
               
                if (err < smallestError) {
                        print(i)
                        smallestError = err
                        print(paste(depth,rounds,err))
                }
            }
        }
} 

# +
cv <- 30
cvDivider <- floor(nrow(aftertrain) / (cv+1))

smallestError <- 100
for (depth in seq(1,10,1)) { 
        for (rounds in seq(1,20,1)) {
            for(eta in seq(0.0001,0.01,0.01)) {
                totalError <- c()
                indexCount <- 1
                for (cv in seq(1:cv)) {
                        # assign chunk to data test
                        dataTestIndex <- c((cv * cvDivider):(cv * cvDivider + cvDivider))
                        dataTest <- aftertrain[dataTestIndex,]
                        # everything else to train
                        dataTrain <- aftertrain[-dataTestIndex,]
                       
                        bst <- xgboost(data = as.matrix(aftertrain),
                                       label = aftertrain,
                                       max.depth=depth, nround=rounds,
                                       objective = "reg:squarederror")
                        predictions <- predict(bst, as.matrix(aftertest), outputmargin=TRUE)
                       
                        err <- rmse(as.numeric(aftertest), as.numeric(predictions))
                        totalError <- c(totalError, err)
                }
                if (mean(totalError) < smallestError) {
                        smallestError = mean(totalError)
                        print(paste(depth,rounds,smallestError))
                } 
            }
        }
} 

# +
# Test both models out on full data set
# assign everything else to test
bst <- xgboost(data = as.matrix(reduced_data),
               label = reduced_data,
               max.depth=4, nround=19, objective = "reg:squarederror", verbose=0)
pred <- predict(bst, as.matrix(testset), outputmargin=TRUE)
rmse(as.numeric(testset), as.numeric(pred))

bst <- xgboost(data = as.matrix(reduced_data),
               label = reduced_data,
               max.depth=3, nround=20, objective = "reg:squarederror", verbose=0)
pred <- predict(bst, as.matrix(testset), outputmargin=TRUE)
rmse(as.numeric(testset), as.numeric(pred))
# -

boostAvant <- xgboost(as.matrix(datatrain), label = ytrain,
                 eta = 0.001,
                 max.depth = 5,
                 nrounds = 10000,
                 objective = "reg:squarederror")



trainCV<-cbind(afterYtrain, aftertrain)

# Helper packages
library(dplyr)
library(tidymodels)
library(vip)         # for variable importance plots

# +
tune_spec <- 
  decision_tree(
    cost_complexity = tune(),
    tree_depth = tune()
  ) %>% 
  set_engine("rpart") %>% 
  set_mode("regression")

tune_spec
# -

tree_grid <- grid_regular(cost_complexity(),
                          tree_depth(),
                          levels = 5)
#The function grid_regular() is from the dials package. It chooses sensible values to try for each hyperparameter; here, we asked for 5 of each. Since we have two to tune, grid_regular() returns 5 
#×5 = 25 different possible tuning combinations to try in a tidy tibble format.

tree_grid

tree_grid %>% 
  count(tree_depth)

set.seed(234)
cell_folds <- vfold_cv(trainCV)

# +
set.seed(345)

tree_wf <- workflow() %>%
  add_model(tune_spec) %>%
  add_formula(afterYtrain ~ .)

tree_res <- 
  tree_wf %>% 
  tune_grid(
    resamples = cell_folds,
    grid = tree_grid
    )

tree_res
# -

tree_res %>% 
  collect_metrics()

tree_res %>%
  collect_metrics() %>%
  mutate(tree_depth = factor(tree_depth)) %>%
  ggplot(aes(cost_complexity, mean, color = tree_depth)) +
  geom_line(size = 1.5, alpha = 0.6) +
  geom_point(size = 2) +
  facet_wrap(~ .metric, scales = "free", nrow = 2) +
  scale_x_log10(labels = scales::label_number()) +
  scale_color_viridis_d(option = "plasma", begin = .9, end = 0)

tree_res %>%
  show_best("rmse")

# +
best_tree <- tree_res %>%
  select_best("rmse")

best_tree

# +
final_wf <- 
  tree_wf %>% 
  finalize_workflow(best_tree)

final_wf

StarcoderdataJupyterScriptsDedupFiltered