# 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)

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

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

check_agreement <- read_csv(file = '../input/check_agreement.csv')
check_correlation <- read_csv(file = '../input/check_correlation.csv')


# 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)))

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

# Matthew Bascom
# http://www.kaggle.com/Bascom

# Tau Lepton Flavour Violation Final Script.R

# Competition: "Flavours of Physics: Finding Tau -> 3mu"
# http://www.kaggle.com/c/flavours-of-physics

# Sunday, October 4, 2015

##############
# Objectives # 
##############

# The objective of this competition was to predict a rare decay 
# phenomenon observed at CERN's Large Hadron Collider in 
# Geneva, Switzerland. 

# The Kolmogorov-Smirnov Test and Cramer-von
# Mises Test were used to verify model design requirements. The 
# classifier model must illustrate non-correlation to tau-particle 
# mass. This is verified using the CvM-test. 
# The model must be indifferent to real events and events 
# generated in CERN's Monte Carlo simulations of proton collisions. 
# This is verified using the KS-test.

# Final evaluation is performed with the weighted area under the 
# Receiver Operating Characteristic (ROC) curve. The weighted ROC
# metric is used to evaluate the classifier model's performance, 
# incorporating the predictor's True-Positive and False-Positive Rates.
# The weighted AUC assigns different weights to different to different
# bins of signal efficiency.

###########################
# Summary and Conclusions #
###########################

# For this competition I utilized a logit generalized linear model 
# for generating binomial event probabilities. My choice of model was 
# influenced in part by minimizing the AIC criteria. You may explore
# the features of the GLM by loading the 'LHCb.GLM.10.12.RData' file.

# I used 27 of the 47 available variables, and several variable 
# interaction terms. Though using min_ANNmuon was explicitly forbidden 
# for fitting the model, I filtered the data to account for this background 
# before fitting the model, so that min_ANNmuon > .4, as is suggested in [1].

# Future analysis may include exploring different
# interaction combinations and, perhaps, fewer terms.

# The Kolmogorov-Smirnov Test Statistic ensures that my model is 
# not sensitive to differences between simulation and real event data. 

# The null hypothesis of the KS-test is that the model is sensitive 
# to features not perfected modeled in simulation data, such as 
# system data that discerns muons from other particles.
# If KS < .09, the null hypothesis may be rejected. 
# My model produces a KS value of 0.0252. 
# The null is rejected.

# The Cramer-von Mises test statistic checks to ensure that the 
# classifier model is uncorrelated to tau-particle mass. As stated in 
# Blake, et. al. (2015), "such correlations can cause an artificial 
# signal-like mass peak or lead to incorrect background estimations."

# The null hypothesis is that the classifier model correlates to tau mass. 
# If CvM < .002, the null hypothesis may be rejected. 
# My model produces a CvM value of 0.000951. 
# The null is rejected.

# Model AUC and Weighted AUC on Training Data:
# [1] 0.9067686	# AUC
# [1] 0.9699447	# Weighted AUC

#############
# Procedure #
#############

# Load libraries.

library(stats)
library(ROCR)
library(caTools)					
#source('WeightedROC II.R')
						# devtools::install_github("tdhock/WeightedROC")
						# Author: Toby Dylan Hocking
						# library(WeightedROC)

LHCv.GLM.10.12 <- glm(formula = signal ~ p0_pt + p0_p + p0_eta + p1_pt + p1_p +
    p1_eta + p2_pt + p2_p + p2_eta + FlightDistance + IPSig +
    VertexChi2 + pt + DOCAone + DOCAtwo + IP_p1p2 + DOCAthree +
    IP_p0p2 + isolationa + isolationb + isolationc + isolationd +
    isolatione + isolationf + iso + FlightDistanceError + dira +
    p0_pt:p0_p + p0_pt:p0_eta + p0_p:p0_eta + p1_pt:p1_p + p1_p:p1_eta +
    p2_pt:p2_p + p2_p:p2_eta + FlightDistance:IPSig + FlightDistance:VertexChi2 +
    IPSig:VertexChi2 + FlightDistance:pt + IPSig:pt + VertexChi2:pt +
    DOCAtwo:IP_p1p2 + DOCAthree:IP_p0p2 + isolationa:isolationb +
    isolationa:isolationc + isolationb:isolationc + isolationa:isolationd +
    isolationb:isolationd + isolationc:isolationd + isolationa:isolatione +
    isolationb:isolatione + isolationc:isolatione + isolationd:isolatione +
    isolationa:isolationf + isolationb:isolationf + isolationc:isolationf +
    isolationd:isolationf + isolatione:isolationf + isolationa:iso +
    isolationb:iso + isolationc:iso + isolationd:iso + isolatione:iso +
    isolationf:iso + FlightDistanceError:dira + p0_pt:p0_p:p0_eta +
    FlightDistance:IPSig:pt + isolationa:isolationb:isolationc +
    isolationa:isolationb:isolationd + isolationa:isolationc:isolationd +
    isolationb:isolationc:isolationd + isolationa:isolationb:isolatione +
    isolationa:isolationc:isolatione + isolationb:isolationc:isolatione +
    isolationa:isolationd:isolatione + isolationb:isolationd:isolatione +
    isolationc:isolationd:isolatione + isolationa:isolationb:isolationf +
    isolationa:isolationc:isolationf + isolationb:isolationc:isolationf +
    isolationa:isolationd:isolationf + isolationb:isolationd:isolationf +
    isolationc:isolationd:isolationf + isolationa:isolatione:isolationf +
    isolationb:isolatione:isolationf + isolationc:isolatione:isolationf +
    isolationd:isolatione:isolationf + isolationa:isolationb:iso +
    isolationa:isolationc:iso + isolationb:isolationc:iso + isolationa:isolationd:iso +
    isolationb:isolationd:iso + isolationc:isolationd:iso + isolationa:isolatione:iso +
    isolationb:isolatione:iso + isolationc:isolatione:iso + isolationd:isolatione:iso +
    isolationa:isolationf:iso + isolationb:isolationf:iso + isolationc:isolationf:iso +
    isolationd:isolationf:iso + isolatione:isolationf:iso + isolationa:isolationb:isolationc:isolationd +
    isolationa:isolationb:isolationd:isolatione + isolationb:isolationc:isolationd:isolatione +
    isolationa:isolationb:isolationc:isolationf + isolationa:isolationb:isolationd:isolationf +
    isolationa:isolationc:isolationd:isolationf + isolationb:isolationc:isolationd:isolationf +
    isolationa:isolationb:isolatione:isolationf + isolationa:isolationc:isolatione:isolationf +
    isolationa:isolationd:isolatione:isolationf + isolationb:isolationd:isolatione:isolationf +
    isolationc:isolationd:isolatione:isolationf + isolationa:isolationb:isolationc:iso +
    isolationa:isolationb:isolationd:iso + isolationa:isolationc:isolationd:iso +
    isolationa:isolationb:isolatione:iso + isolationb:isolationc:isolatione:iso +
    isolationa:isolationd:isolatione:iso + isolationb:isolationd:isolatione:iso +
    isolationa:isolationb:isolationf:iso + isolationb:isolationd:isolationf:iso +
    isolationa:isolatione:isolationf:iso + isolationb:isolatione:isolationf:iso +
    isolationd:isolatione:isolationf:iso + isolationa:isolationb:isolatione:isolationf:iso,
    family = binomial(link = "logit"), data = train)

test_probs <- predict(object = LHCv.GLM.10.12, newdata = test, type = 'response')
print(head(test_probs))
print(head(test$signal))

agreement_probs <- predict(object = LHCv.GLM.10.12, newdata = check_agreement, type = 'response')
correlation_probs <- predict(object = LHCv.GLM.10.12, newdata = check_correlation, type = 'response')

# Functions from 'WeightedROC' package.
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
}

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_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
}

#ks <- compute_ks(agreement_probs[which(check_agreement$signal == 0),2],    # If probs are in matrix with ids
#                agreement_probs[which(check_agreement$signal == 1),2],
ks <- compute_ks(agreement_probs[which(check_agreement$signal == 0)],
                agreement_probs[which(check_agreement$signal == 1)],
                check_agreement$weight[which(check_agreement$signal==0)],
                check_agreement$weight[which(check_agreement$signal==1)])

# R package stats: Function stats::glm was used to generate the classifier model.
# R package caTools: moving window statistics, ROC, AUC, and trapezoidal rule.
# R package ROCR: visualizing the performance of scoring classifiers.
# R package WeightedROC: calculate weighted area under the ROC curve.

print(train[1,])
print(test[1,])


print("Confusion Matrix of prediction on training set")
table(LHCv.GLM.10.12$data$signal, LHCv.GLM.10.12$fitted.values > .5) 
#    FALSE TRUE
#  0   793  522
#  1   209 4752


###########################################
# Kolmogorov-Smirnov Check Agreement Test #
###########################################

# Read in the check_agreement predictions file and calculate the 
# Kolmogorov-Smirnov Test statistic. 

# The KS-stat ensures that my model is not sensitive to differences 
# between simulation and real event data. 

# The null hypothesis is that the model will be sensitive 
# to features not perfected modeled in simulation data. 
# If ks < .09, the null hypothesis may be rejected. 
# My model produces a ks value of 0.0252367434224073.

print("Check Agreement Test")

#head(agreement_probs[which(check_agreement$signal == 0),2])

#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("Kolmogorov-Smirnov Metric ",ks," ks < 0.09 = ",ks < 0.09))
# [1] "KS METRIC 0.0252367434224073 ks < 0.09 = TRUE"
# [1] "2015-10-04 06:39:33 PDT"

###########################################
# Cramer-von Mises Check Correlation Test #
###########################################

# The CvM-stat ensures that my model is not correlated to tau mass.

# The null hypothesis is that the model will correlate to tau mass. 
# If CvM < .002, the null hypothesis may be rejected. 
# My model produces a CvM value of 0.000950743829553582.

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
}

cvm <- function (subindices,total_events){ 
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){ 
predictions = predictions[order(masses)] 
predictions = order(order(predictions)) 
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) 
} 

cvm <- compute_cvm(predictions = correlation_probs, masses = check_correlation[, 'mass'])
#cvm <- compute_cvm(predictions = correlation_probs[,'prediction'], masses = check_correlation[, 'mass'])
print("Check Correlation Test")

print (paste0("Cramer-von Mieses Metric ",cvm, " cvm < 0.002 = ",cvm < 0.002))
# [1] "CvM metric 0.000950743829553582 cvm < 0.002 = TRUE"

#######################################
# Weighted Area Under the Curve (AUC) #
#######################################

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)){ 
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))){ 
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 
}

# ROCR::prediction generates an S4 class object.
pred.GLM.10.12 <- prediction(predictions = LHCv.GLM.10.12$fitted.values,
					labels = LHCv.GLM.10.12$y)	
					
auc_train <- as.numeric(performance(pred.GLM.10.12, "auc")@y.values)	# S4-class slot @y.values					

print("AUC for prediction on training set")
auc_train#
# [1] 0.9067686

print("Weighted AUC for prediction on training set")
roc_auc_truncated(labels = LHCv.GLM.10.12$data$signal, predictions = LHCv.GLM.10.12$fitted.values)
# [1] 0.9699447

# caTools::trapz is a function for  the trapezoidal rule and is used