# You can write R code here and then click "Run" to run it on our platform

# The readr library is the best way to read and write CSV files in R
library(readr)
library(caret)
library(caretEnsemble)
library(Metrics)
library(caTools)
library(xgboost)
library(Matrix)
library(ROCR)
library(pROC)
library(gbm)

#Define Evaluation Metrics
print("Defining Evaluation Metrics")

#Rolling window: take window with definite size through the array
#param data: array-like
#param window_size: size
#return: the sequence of windows
#Example: data = array(1, 2, 3, 4, 5, 6), window_size = 4
#Then this function return array(array(1, 2, 3, 4), array(2, 3, 4, 5), array(3, 4, 5, 6))

rolling_window <- function(data,window_size){
  rw <- NULL
  for(i in 1:(length(data)-window_size+1)){
    rw[[i]] <- data[i:(i+window_size-1)]
  }
  rw
}

#Following function has been copied from R package "WeightedROC
# because it is not available on CRAN. and not available on Kaggle R Scripts
# for using on a local machine, users are recommended to install the relevant package 
# by running following command:
# devtools::install_github("tdhock/WeightedROC")
# Author: Toby Dylan Hocking
# library(WeightedROC)

WeightedROC <- function (guess, label, weight = rep(1, length(label))) 
{
  if (is.factor(label)) {
    label <- as.integer(label)
  }
  stopifnot(is.numeric(label))
  if (all(label %in% c(0, 1))) {
    label[label == 0] <- -1
  }
  if (all(label %in% c(1, 2))) {
    label[label == 1] <- -1
    label[label == 2] <- 1
  }
  stopifnot(label %in% c(-1, 1))
  stopifnot(is.numeric(guess))
  stopifnot(length(label) == length(guess))
  stopifnot(is.numeric(weight))
  stopifnot(length(label) == length(weight))
  #stopifnot(weight > 0)
  ord <- order(guess)
  y <- label[ord]
  w <- weight[ord]
  y.hat <- guess[ord]
  is.positive <- y == 1
  is.negative <- y == -1
  w.positive <- w.negative <- w
  w.positive[is.negative] <- 0
  w.negative[is.positive] <- 0
  cum.positive <- cumsum(w.positive)
  cum.negative <- cumsum(w.negative)
  is.end <- c(diff(y.hat) != 0, TRUE)
  n <- length(y)
  TPR <- c(1, 1 - cum.positive[is.end]/cum.positive[n])
  FPR <- c(1, 1 - cum.negative[is.end]/cum.negative[n])
  d <- data.frame(TPR, FPR)
  d
}


cvm <- function (subindices,total_events){
  #   Compute Cramer-von Mises metric.
  #   Compared two distributions, where first is subset of second one.
  #   Assuming that second is ordered by ascending
  #   :param subindices: indices of events which will be associated with the first distribution
  #   :param total_events: count of events in the second distribution
  #   :return: cvm metric
  target_distribution = (1:(1+total_events))/total_events
  subarray_distribution = cumsum((hist(x = subindices,breaks = seq(0.5,total_events+1.5,1),plot = FALSE))$counts)
  subarray_distribution = subarray_distribution/subarray_distribution[length(subarray_distribution)]
  mean((target_distribution-subarray_distribution)^2)
}


compute_cvm <- function(predictions,masses, n_neighbours = 200,step = 50){
  # Computing Cramer-von Mises (cvm) metric on background events: take average of cvms calculated for each mass bin.
  # In each mass bin global prediction's cdf is compared to prediction's cdf in mass bin.
  # :param predictions: array-like, predictions
  # :param masses: array-like, in case of Kaggle tau23mu this is reconstructed mass
  # :param n_neighbours: count of neighbours for event to define mass bin
  # :param step: step through sorted mass-array to define next center of bin
  # :return: average cvm value
  
  # First, reorder by masses
  predictions = predictions[order(masses)]
  # Second, replace probabilities with order of probability among other events
  predictions = order(order(predictions))
  # Now, each window forms a group, and we can compute contribution of each group to CvM
  cvms = c()
  rw = rolling_window(predictions,window_size = n_neighbours)
  for (window in rw[seq(1,length(rw),step)]){
    cvms = c(cvms,cvm(subindices = window,total_events = length(predictions)))
  }
  mean(cvms)
}

roc_curve_splitted <- function(data_zero, data_one, sample_weights_zero, sample_weights_one){
  #   Compute roc curve
  #   :param data_zero: 0-labeled data
  #   :param data_one:  1-labeled data
  #   :param sample_weights_zero: weights for 0-labeled data
  #   :param sample_weights_one:  weights for 1-labeled data
  #   :return: roc curve
  labels = c(rep(0,length(data_zero)),rep(1,length(data_one)))
  weights = c(sample_weights_zero, sample_weights_one)
  data_all = c(data_zero, data_one)
  tpr.fpr = WeightedROC(guess = data_all,label = labels,weight = weights)
  tpr.fpr
}
compute_ks <- function(data_prediction, mc_prediction, weights_data, weights_mc){
  #   Compute Kolmogorov-Smirnov (ks) distance between real data predictions cdf and Monte Carlo one.
  #   :param data_prediction: array-like, real data predictions
  #   :param mc_prediction: array-like, Monte Carlo data predictions
  #   :param weights_data: array-like, real data weights
  #   :param weights_mc: array-like, Monte Carlo weights
  #   :return: ks value
  weights_data = weights_data/sum(weights_data)
  weights_mc = weights_mc/sum(weights_mc)
  #weights_data = weights_data - min(weights_data)+0.00001
  #weights_mc = weights_mc - min(weights_mc)+0.00001
  tpr.fpr = roc_curve_splitted(data_prediction,mc_prediction,weights_data,weights_mc)
  Dnm = max(abs(tpr.fpr[,1]-tpr.fpr[,2]))
  Dnm
}

roc_auc_truncated <- function(labels, predictions, tpr_thresholds=c(0.2, 0.4, 0.6, 0.8),
                              roc_weights=c(4, 3, 2, 1, 0)){
  #   labels = as.numeric(valResponse)
  #   predictions = predValxgbTree_caret[,2]
  #   tpr_thresholds=c(0.2, 0.4, 0.6, 0.8)
  #   roc_weights=c(4, 3, 2, 1, 0)
  #   Compute weighted area under ROC curve.
  #   :param labels: array-like, true labels
  #   :param predictions: array-like, predictions
  #   :param tpr_thresholds: array-like, true positive rate thresholds delimiting the ROC segments
  #   :param roc_weights: array-like, weights for true positive rate segments
  #   :return: weighted AUC
  tpr.fpr = WeightedROC(guess = predictions,label = labels)
  tpr.fpr = tpr.fpr[order(tpr.fpr$TPR,tpr.fpr$FPR),]
  area = 0
  tpr_thresholds = c(0,tpr_thresholds,1)
  for(index in (2:length(tpr_thresholds))){
    #     index = 6
    tpr.cut = sapply(tpr.fpr[,1],function(x){min(x,tpr_thresholds[index])})
    tpr.previous = sapply(tpr.fpr[,1],function(x){min(x,tpr_thresholds[index-1])})
    area = area+roc_weights[index-1]*(trapz(tpr.fpr$FPR,tpr.cut)-trapz(tpr.fpr$FPR,tpr.previous))
  }
  area = area/sum((tpr_thresholds[-1] - tpr_thresholds[-length(tpr_thresholds)])*roc_weights)
  area
}


# The competition datafiles are in the directory ../input
print("List the files we have available to work with:")
system("ls ../input")

# Read competition data files:
print("Read competition data files:")
train <- read_csv("../input/training.csv")
test <- read_csv("../input/test.csv")
check_correlation <- read_csv("../input/check_correlation.csv")
check_agreement <- read_csv("../input/check_agreement.csv")

#summary(check_agreement$weight)
#summary(check_agreement$weight-min(check_agreement$weight))
#summary(check_agreement$weight-min(check_agreement$weight+0.00001))

# Write to the log:
#cat(sprintf("Training set has %d rows and %d columns\n", nrow(train), ncol(train)))
#cat(sprintf("Test set has %d rows and %d columns\n", nrow(test), ncol(test)))
#Inspect data
#str(train)
#summary(train)
#str(read_csv("../input/sample_submission.csv"))
#table(train$signal)
train$signal = as.factor(train$signal)
levels(train$signal) <- make.names(levels(train$signal))

#give a column missing from right side input in camparison to left side
"%w/o%" <- function(x, y) x[!x %in% y] #--  x without y

missingcolnames = colnames(train) %w/o% colnames(test)


xx1 = colnames(train) %w/o% c("id","SPDhits","IPSig","IP","isolationb","isolationc",missingcolnames)
set.seed(244)
split = sample.split(train$signal, SplitRatio = 0.8)
predictor = (subset(train, split == TRUE))[,which(colnames(train) %in% xx1)]
response = (subset(train$signal, split == TRUE))
valPredictor = (subset(train, split == FALSE))[,which(colnames(train) %in% xx1)]
valResponse = (subset(train$signal, split == FALSE))
testpredictor = test[,which(colnames(test) %in% xx1)]

trainControlxgbTree <- trainControl(method = "boot",number = 4,repeats = 1,#summaryFunction = twoClassSummary,
                                    verboseIter = TRUE,
                                    preProcOptions = list(thresh = 0.999,ICAcomp = 111)
)
#summary(predictor)
#summary(response)

print(" Building xgbTree caret model")
#gtrain = xgb.DMatrix(data = sparse.model.matrix(~., data = predictor),label = response)
#gval = xgb.DMatrix(data = sparse.model.matrix(~., data = valPredictor),label = valResponse)

xgbTree_caret <- train(x = predictor,
                       y = response,
                       method = "xgbTree",
                       tuneGrid = expand.grid(nrounds = 10*(1:15),
                                              eta = c(0.276),
                                              max_depth = c(3)),                            
                       trControl = trainControlxgbTree,
                       #watchlist = list(val = xgval,train = xgtrain),
                       #eval_metric = "error",
                       min_child_weight = 10,
                       #early.stop.round = 50,
                       #printEveryN = 20,
                       subsample = 0.7,
                       #verbose = 1,
                       colsample_bytree =0.8,
                       base_score = 0.5,
                       nthread =8#,
                       #preProcess = c("center","scale")
                       #                        ,tuneLength = 100
)

gbm_formula=as.formula(paste0("signal ~ ", paste(xx1,collapse = " + ")))
gb=gbm(gbm_formula,         # formula
       data=train,                   # dataset
       
       distribution="gaussian",     # see the help for other choices
       n.trees=1000,                # number of trees
       shrinkage=0.05,              # shrinkage or learning rate,
       # 0.001 to 0.1 usually work
       interaction.depth=3,         # 1: additive model, 2: two-way interactions, etc.
       bag.fraction = 0.5,          # subsampling fraction, 0.5 is probably best
       train.fraction = 0.5,        # fraction of data for training,
       # first train.fraction*N used for training
       n.minobsinnode = 10,         # minimum total weight needed in each node
       cv.folds = 3,                # do 3-fold cross-validation
       keep.data=TRUE,              # keep a copy of the dataset with the object
       verbose=FALSE)

print("GB and XGB both trained")

plot(xgbTree_caret)
print("plot complete")

predTrainxgbTree_caret = predict.train(xgbTree_caret,newdata = predictor, type = "prob")
print("Confusion Matrix of prediction on training set")
table(response, predTrainxgbTree_caret[,2]>0.5)

auc_train = as.numeric(performance(prediction(predTrainxgbTree_caret[,2],response),"auc")@y.values)
print("Plot ROC for training set")

plot.roc(response, predTrainxgbTree_caret[,2],print.thres = seq(00.05,0.95, by = 0.05))
print("AUC for prediction on training set")
auc_train# # 

print("Weighted AUC for prediction on training set")
roc_auc_truncated(labels = as.numeric(response),predictions = predTrainxgbTree_caret[,2])

predValxgbTree_caret = predict(xgbTree_caret,newdata = valPredictor, type = "prob")
print("Confusion Matrix of prediction on validation set")
table(valResponse, predValxgbTree_caret[,2]>0.5)
auc_val_xgb = as.numeric(performance(prediction(predValxgbTree_caret[,2],valResponse),"auc")@y.values)
print("Plot ROC for training set")
plot.roc(valResponse, predValxgbTree_caret[,2],print.thres = seq(00.05,0.95, by = 0.05))
print("AUC for prediction on validation set")
auc_val_xgb# # 
print("Weighted AUC for prediction on validation set")
roc_auc_truncated(labels = as.numeric(valResponse),predictions = predValxgbTree_caret[,2])

print("Check Agreement Test")
agreement_probs = predict(xgbTree_caret,newdata = check_agreement,type = "prob")
ks = compute_ks(agreement_probs[which(check_agreement$signal == 0),2],
                agreement_probs[which(check_agreement$signal == 1),2],
                check_agreement$weight[which(check_agreement$signal==0)],
                check_agreement$weight[which(check_agreement$signal==1)])
print(paste0("KS METRIC ",ks," ks < 0.09 = ",ks < 0.09))

print("Check CvM Test")
correlation_probs =predict(xgbTree_caret,check_correlation,type = "prob")
cvm = compute_cvm(predictions = correlation_probs[,2],masses = check_correlation$mass)
print (paste0("CvM metric ",cvm, " cvm < 0.002 = ",cvm < 0.002))

PredTest_xgbTree_caret = predict(xgbTree_caret,newdata = testpredictor, type = "prob")
PredTest_xgbTree_caret = data.frame(Id = test$id, prediction = PredTest_xgbTree_caret[,2])
write.csv(PredTest_xgbTree_caret, "Submission_xgbTree_caret.csv", row.names= FALSE)
system("gzip --force Submission_xgbTree_caret.csv")

print("Modelling Complete")

# Generate output files with write_csv(), plot() or ggplot()
# Any files you write to the current directory get shown as outputs