---
title: "...or Doping?"
author: "Vicens Gaitan"
date: "Sunday, October 18, 2015"
output: html_document
---

```{r,warning=FALSE, message=FALSE}
library(xgboost)
library(Matrix)
library(caret)
library(ggplot2)
library(readr)
#source("C:/local/_TOOLS/kaggle/tmmm/metrics.R")

path="C:/local/_TOOLS/kaggle/tmmm/"
```
```{r,echo=FALSE}
library(pracma)
#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)
  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)
  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
}

ptau<-function(d){
  # muon mass in MeV/c^2
  mmu = 105.6583715
 
  # calculate tau energy
  tau_e = sqrt(d$p0_p ** 2 + mmu**2) + sqrt(d$p1_p ** 2 + mmu**2) +sqrt(d$p2_p ** 2 + mmu**2)
  # calculate pz of tau candidate
  tau_pz = d$p0_pt * sinh(d$p0_eta) + d$p1_pt * sinh(d$p1_eta)+d$p2_pt * sinh(d$p2_eta)
  # calculate momentum of tau candidate
  tau_p =sqrt(d$pt ** 2 + tau_pz ** 2)
  # calculate eta of tau candidate
  tau_eta = asinh(tau_pz / d$pt)
  # calculate mass of tau candidate
  tau_m2=tau_e ** 2 - tau_p ** 2
  tau_m2[tau_m2<0]=0
  tau_m=sqrt(tau_m2)
  #M = tau_p*LifeTime*c/FlightDistance
  #c Speed of light
  c= 299.792458
  tau_m=tau_p*d$LifeTime*c/d$FlightDistance
  return(list(tau_e,tau_pz,tau_p,tau_eta,tau_m))
}  
ftrain<-function(data,y_d,params,folds,mask2,nrounds,max,pre=F,scan=F){
  # do the folds
  for(i in 1: length(folds)){
    # get the sampling
    # and construct train and test matrix out of train data
    smpl = folds[[i]]
    g_train = data[-smpl,]
    g_test = data[smpl,]
    y_train = y_d[-smpl]
    y_test = y_d[smpl]
    
    dtrain0=xgb.DMatrix(as.matrix(g_train[,mask2]), label =y_train)
    dtest0=xgb.DMatrix(as.matrix(g_test[,mask2]), label =y_test)

    watchlist0=list(validation=dtest0,train=dtrain0)
    bst0 = xgb.train(param=params, data = dtrain0,nrounds=nrounds,watchlist0,print.every.n = 100,verbose=0,maximize=max,early.stop.round = 100)
    
    #importance_matrix <- xgb.importance(model = bst0,feature_names=names(data[,mask2]))
    #print(importance_matrix)
    
    best0=bst0$bestInd
    
    ptrain[smpl]<-predict(bst0,dtest0,ntreelimit=best0)
    if(pre){
      pagree[,i]=predict(bst0,as.matrix(agree[,mask2]),ntreelimit=best0)
      pcorr[,i]=predict(bst0,as.matrix(corr[,mask2]),ntreelimit=best0)
      #ptest[,i]=predict(bst0,as.matrix(test[,mask2]),ntreelimit=best0)
    }
  }   
  
  return(list(ptrain,pagree,pcorr))
}

fscan<-function(data,y_d,params,folds,mask2,nrounds,max){
  # do the folds
  i=1
    # get the sampling
    # and construct train and test matrix out of train data
    smpl = folds[[i]]
    g_train = data[-smpl,]
    g_test = data[smpl,]
    y_train = y_d[-smpl]
    y_test = y_d[smpl]
    
    dtrain0=xgb.DMatrix(as.matrix(g_train[,mask2]), label =y_train)
    dtest0=xgb.DMatrix(as.matrix(g_test[,mask2]), label =y_test)

    watchlist0=list(validation=dtest0,train=dtrain0)
    bst0 = xgb.train(param=params, data = dtrain0,nrounds=nrounds,watchlist0,print.every.n = 100,verbose=0,maximize=max,early.stop.round = 100)
    
    #importance_matrix <- xgb.importance(model = bst0,feature_names=names(data[,mask2]))
    #print(importance_matrix)
    
    best0=bst0$bestInd
    out=rep(0,best0/10)
    out=cbind(out,out,out)
    names(out)=c("ntree","AUCw","KS")
    out=as.data.frame(out)
    k=1
    for(j in seq(1,best0,10)){
         pred<-predict(bst0,dtest0,ntreelimit=j)
         pagree[,i]=predict(bst0,as.matrix(agree[,mask2]),ntreelimit=j)
         mask=(g_test$min_ANNmuon>.4 )
         ah=roc_auc_truncated(y_test[mask],pred[mask])
         d0=pagree[agree$signal==0,i]
         d1=pagree[agree$signal==1,i]
         w0=agree$weight[agree$signal==0]
         w1=agree$weight[agree$signal==1]
         kks=compute_ks(d0,d1,w0,w1)
         out[k,1]=j
         out[k,2]=ah
         out[k,3]=kks  
         k=k+1
      }
  
  return(out)
}
evaluation<-function(holdout,ragg,rcorr){
  mask=(data$min_ANNmuon>.4 )
  ah=roc_auc_truncated(y_d[mask],holdout[mask])
  cc=compute_cvm(rcorr,corr$mass)
  d0=ragg[agree$signal==0]
  d1=ragg[agree$signal==1]
  w0=agree$weight[agree$signal==0]
  w1=agree$weight[agree$signal==1]
  kks=compute_ks(d0,d1,w0,w1)
  print(paste0("AUCw=",ah," CVM=",mean(cc)," K-S=",kks))
}
```

```{r,warning=FALSE, message=FALSE}
train<- read_csv("../input/training.csv")
m=ptau(train)
train$mm=as.data.frame(m[5])[,1]
agree <- read_csv("../input/check_agreement.csv")

m=ptau(agree)
agree$mm=as.data.frame(m[5])[,1]

corr <- read_csv("../input/check_correlation.csv")
m=ptau(corr)
corr$mm=as.data.frame(m[5])[,1]

y_d=train$signal
ptrain=y_d
#ptest=as.data.frame(test[,1])
pagree=as.data.frame(agree[,1])
pcorr=as.data.frame(corr[,1])

cv_fold_count=2
set.seed(2603)

folds = createFolds(as.factor(y_d),cv_fold_count);

param <- list(
  "objective" = "binary:logistic",
  "eval_metric"="auc",
  "max_depth"=5,
  "eta"=.05
)


data=train

```

So... what can be done, useful from a physical point of view?

  a) Not to use reconstructed mass, nor features allowing easy mass reconstruction
  b) Try to not use variable regions for which the Monte Carlo simulation doesn't agree with real data
  
The first one is fulfilled in we don't explicit add the reconstructed mass, because the classifier has hard time to built it from the data without external guidance ( tau mass is a tiny effect in multi GeV energy and momentum values)

If we want to use b) we have to break the rules and use a data sample including data and MC, different form the training set: the agreement set.

The idea is to "dope" (in the semiconductor meaning) the training set with a small number of Monte Carlo events from the agreement sample, but labeled as background. This disallow the classifier to pick features discriminating data and Monte Carlo. There are two parameters that regularize the learning: The number of "doping" events and the complexity of the classifier (for instance number of trees)


Let's select the candidates for doping

```{r,warning=FALSE, message=FALSE}
train2=agree[agree$signal==1,]
train2$production=0
train2$min_ANNmuon=0
train2$mass=0
train2$weight<-NULL
train2$signal=0
```

And define a function to see the effect of adding n new events:
```{r,warning=FALSE, message=FALSE}
scan_Regul<-function(ns){
  data=rbind(train,train2[1:ns,])
  y_d=data$signal
  folds = createFolds(as.factor(y_d),cv_fold_count);
  mask2=c(2:47)
  res1=fscan(data,y_d,param,folds,mask2,200,T)
  names(res1)=c("ntree","AUCw","KS")
  res1=res1[2:dim(res1)[1],]
  ggp=ggplot( data =res1 , aes(AUCw, KS ,color=KS )) +geom_point(size = 5)+ scale_colour_gradientn(colours=c("blue", "red"))+ggtitle(paste0("Number of events " ,ns))
 return( list(res1,ggp))
}
```


Let's see the model without doping

```{r,warning=FALSE, message=FALSE}
res=scan_Regul(0)
res[1]
res[2]
```

The plot shows for each number of trees , the obtained K-S value, in vertical, and the AUCw in horizontal. 

For 2000 events, we obtain

```{r,warning=FALSE, message=FALSE}
res=scan_Regul(2000)
res[1]
res[2]
```

As can be see, there is an intermediate value optimizing the AUCw subject to the constrain in K-S.

Optimizing simultaneously for Nvents and Ntrees , shows a global optimum around N=3460


```{r,warning=FALSE, message=FALSE}
res=scan_Regul(3460)
res[1]
res[2]
```

Tunning the xgboost parameters AUCw=0.989x can be obtained 
This model is not correlated with mass and pass the agreement test

```{r,warning=FALSE, message=FALSE}
  ns=3460

  data=rbind(train,train2[1:ns,])
  y_d=data$signal
  folds = createFolds(as.factor(y_d),cv_fold_count);
  mask2=c(2:47)
  res=ftrain(data,y_d,param,folds,mask2,500,T,T)
  holdout0=as.data.frame(res[1])[,1]
  ragree0=rowMeans(as.data.frame(res[2]))
  rcorr0=rowMeans(as.data.frame(res[3]))
  evaluation(holdout0,ragree0,rcorr0)


par(mfrow=c(3,1))
plot(train$mm,holdout0[1:dim(train)[1]],pch='.',col=rgb(train$signal,0,0,.3),xlim=c(1700,2100))
plot(agree$mm,ragree0,pch='.',col=rgb(agree$signal,0,0,.2),xlim=c(1700,2100))
plot(corr$mm,rcorr0,pch='.',col=rgb(0,0,0,1),xlim=c(1700,2100))
 

par(mfrow=c(1,2))
mask=agree$signal==0
plot(ragree0[mask],agree$SPDhits[mask],pch='.',col=rgb(agree$signal[mask],0,0,.4))
mask=agree$signal==1
plot(ragree0[mask],agree$SPDhits[mask],pch='.',col=rgb(agree$signal[mask],0,0,.4))
par(mfrow=c(1,1)) 

```