---
title: "Why is so hard to 'learn' the tau mass from kinematics?"
author: "Vicens Gaitan"
date: "Wednesday, October 14, 2015"
output: html_document
---

```{r}
library(xgboost)
library(Matrix)
library(caret)
library(ggplot2)
library(readr)
```
```{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)
}

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
}

```
The final outcome of this competition has been puzzling me in the last days:

Fact 1: The tau invariant mass can be reconstructed with high accuracy from kinematic variables (p0,pt,eta). No need for phi angles as the decay has cylindrical symmetry.

Fact 2: The reconstructed tau mass separates nicely signal and background because the background spectrum has a "hole" for decays coming from a true tau (why????)

So: In theory should  be possible to train a classifier to separate signal and background in an efficient way just Using kinematic variables ( and obviously not passing CVM test...) 

Then: Why rich classifiers as gmb or random forest trained on kinematic variables without explicit mass reconstruction doesn't perform as well as the ones using the mass.

Let's start exploring the two facts and their consequences: 

Firstly: reconstruct the mass using exact relativistic kinematic relations:
```{r}
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))
}  
```

Let's tray on the training data 

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

This analytically calculated mass agree  very well with the existing mass variable in the train dataset
```{r}
plot(train$mm,train$mass,pch='.',col=rgb(0,0,0,.1),xlim=c(1600,1960),ylim=c(1600,1960))
```

Then reconstruct also the mass for the agreement samples and plot it jointly (discriminating signal from background)
```{r}
agree <- read_csv("../input/check_agreement.csv")

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

ggplot(data=train, aes(mm)) +
  geom_histogram(data = subset(train,signal == 1),breaks=seq(1600, 2100, by =5), fill = "red", alpha = 0.3) +
  geom_histogram(data = subset(train,signal == 0),breaks=seq(1600, 2100, by =5),fill = "blue", alpha = 0.3) +
  geom_histogram(data = agree,breaks=seq(1600, 2100, by =5),fill = "green", alpha = 0.3)
```

This plot shows the "Big Trap": mass separates simultaneously  signal vs background, and train vs agreement , but prevents to pass the CVM test

The same plot  for correlation (that fits nicely the train data spectrum) 
```{r}
corr <- read_csv("../input/check_correlation.csv")
m=ptau(corr)
corr$mm=as.data.frame(m[5])[,1]
ggplot(data=train, aes(mm)) + geom_histogram(data = corr,breaks=seq(1600, 2100, by =5),fill = "black", alpha = 0.3) 
```

And finally for the test sample
```{r}
test <- read_csv("../input/test.csv")
m=ptau(test)
test$mm=as.data.frame(m[5])[,1]
ggplot(data=test, aes(mm)) + geom_histogram(data = test,breaks=seq(1600, 2100, by =5),fill = "yellow", alpha = 0.9) 
rm(test)#just to free some memory

```

Now, forgot for the moment the CVM and KS test, and concentrate in the degree of predictability using reconstructed mass.

Lets take a look to the background "hole". From my point of view, there is no reason in introducing this artifact in the training sample, making the problem totally unrealistic as we expect a continuous background superimposed to the tau signal decay
```{r}
ggplot(data=train, aes(mass)) +
  geom_histogram(data = subset(train,signal == 1),breaks=seq(1600, 1960, by =5), fill = "red", alpha = 0.3) +
  geom_histogram(data = subset(train,signal == 0),breaks=seq(1600, 1960, by =5),fill = "blue", alpha = 0.3) 
```

Anyway we'll try to exploit this variable to discriminate signal from background.

First a convenience function to train a classifier using xgboost, returning holdout predictions for the full training sample. For sake of speed, only 2-fold training will be done.

```{r}
ftrain<-function(data,y_d,params,folds,mask2,nrounds,max,pre=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))
}

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

The first model will use the mass appearing in the train sample "mass"

```{r,warning=FALSE, message=FALSE}
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"=.1
)

data=train
y_d=train$signal
mask2=c("mass")

res=ftrain(data,y_d,param,folds,mask2,300,T)
holdout=as.data.frame(res[1])[,1]
mask=(data$min_ANNmuon>.4 )
ah=roc_auc_truncated(y_d[mask],holdout[mask])
```
```{r}
print(format(ah,digits=10,nsmall=20))
```

Excellent... ;)

Now , with the reconstructed mass "mm"

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

plot(train$mm,train$mass,pch='.',col=rgb(train$signal,0,0,.1),xlim=c(1600,1960),ylim=c(1600,1960))
mask2=c("mm")

res=ftrain(data,y_d,param,folds,mask2,300,T)
p=as.data.frame(res[1])[,1]
mask=(data$min_ANNmuon>.4 )
ah=roc_auc_truncated(y_d[mask],p[mask],)
print(format(ah,digits=10,nsmall=20))

```

Not so bad...

Now, as "mm"" is just a deterministic function of ("p0_p","p1_p","p2_p","p0_pt","p1_pt","p2_pt","p0_eta","p1_eta","p2_eta","pt",,"LifeTime","FlightDistance"), let's try a classifier on this variables ( must perform equally well... I hope...)


```{r,warning=FALSE, message=FALSE}
mask2=c("p0_p","p1_p","p2_p","p0_pt","p1_pt","p2_pt","p0_eta","p1_eta","p2_eta","pt","LifeTime","FlightDistance")

res=ftrain(data,y_d,param,folds,mask2,300,T)

p=as.data.frame(res[1])[,1]
mask=(data$min_ANNmuon>.4 )
ah=roc_auc_truncated(y_d[mask],p[mask],)
print(format(ah,digits=10,nsmall=20))
```

but... WTF... wait a moment... let's try to learn "mm" from  ("p0_p","p1_p"...) using a regressor
```{r,warning=FALSE, message=FALSE}
mask2=c("p0_p","p1_p","p2_p","p0_pt","p1_pt","p2_pt","p0_eta","p1_eta","p2_eta","pt","LifeTime","FlightDistance")

param_m <- list(
  "objective" = "reg:linear",
  "max_depth"=5,
  "eta"=.1
)
y_d=train$mm

res=ftrain(data,y_d,param_m,folds,mask2,300,F)
p=as.data.frame(res[1])[,1]
train$mass_1=p
plot(train$mm,train$mass_1,pch='.',col=rgb(train$signal,0,0,.1),xlim=c(1600,1960),ylim=c(1600,1960))

ggplot(data=train, aes(mass_1)) +
  geom_histogram(data = subset(train,signal == 1),breaks=seq(1600, 1960, by =5), fill = "red", alpha = 0.3) +
  geom_histogram(data = subset(train,signal == 0),breaks=seq(1600, 1960, by =5),fill = "blue", alpha = 0.3) 

```


No way... Let's try a simpler thing. mm is just function of 2 variables: total tau energy and total tau momentum. Maybe is possible to learn this simpler relationship. Lets define the dimensionless tau_e and tau_p for tau energy and momentum respectively. Now mm is just 
sqrt(tau_e^2-tau_p^2)
```{r,warning=FALSE, message=FALSE}
mtrain=ptau(train)

m_tau=1776.82#MeV/c^2 

train$tau_e=(as.data.frame(mtrain[1])[,1]/m_tau)
train$tau_p=(as.data.frame(mtrain[3])[,1]/m_tau)

```

And using these variables...
```{r,warning=FALSE, message=FALSE}

mask2=c("tau_e","tau_p")

param_m <- list(
  "objective" = "reg:linear",
  "max_depth"=5,
  "eta"=.1
)
data=train
y_d=train$mm

res=ftrain(data,y_d,param_m,folds,mask2,300,F)
p=as.data.frame(res[1])[,1]

train$mass_3=p

plot(train$mm,train$mass_3,pch='.',col=rgb(train$signal,0,0,.1),xlim=c(1600,1960),ylim=c(1600,1960))

ggplot(data=train, aes(mass_3)) +
  geom_histogram(data = subset(train,signal == 1),breaks=seq(1600, 1960, by =5), fill = "red", alpha = 0.3) +
  geom_histogram(data = subset(train,signal == 0),breaks=seq(1600, 1960, by =5),fill = "blue", alpha = 0.3) 

```

Let's take a look to the variables...

```{r}
ggplot(train, aes(x=tau_p, y=tau_e, color=factor(signal)),alpha=.01) +geom_point()
print(paste0("Correlation=",format(cor(train$tau_p,train$tau_e),digits=10,nsmall=20)))
print(paste0("Mean E=",mean(train$tau_e^2)))
print(paste0("Mean P=",mean(train$tau_p^2)))
```

... to discover that P and E are fully correlated, and the mass is just an effect of 1 part in 2594, hard to discover if you don't know where to look (in this case, not you, but xgboost ) 

In most cases, having 2 variables as correlated as those, suggest to drop one of them, but in this case, let's try a common path , but in a different way.
```{r}
pca=prcomp(train[,c("tau_e","tau_p")])
plot(pca)
pca

rm(pca)
````

PCA analysis tells us that one component fully dominates the data, but in this case, the small one is the interesting, the one correlated to the target. Small component aren't always noise...
```{r}
c1=0.7070256
c2=0.7071880

train$pca1=c1*train$tau_e+c2*train$tau_p
train$pca2=-c1*train$tau_e+c2*train$tau_p

ggplot(train, aes(x=pca1, y=pca2, color=factor(signal)),alpha=.01) +geom_point()
```

We augment the other datasets with the new variables

```{r}


m=ptau(agree)
agree$tau_e=(as.data.frame(m[1])[,1]/m_tau)
agree$tau_p=(as.data.frame(m[3])[,1]/m_tau)
agree$pca1=c1*agree$tau_e+c2*agree$tau_p
agree$pca2=-c1*agree$tau_e+c2*agree$tau_p

m=ptau(corr)
corr$tau_e=(as.data.frame(m[1])[,1]/m_tau)
corr$tau_p=(as.data.frame(m[3])[,1]/m_tau)
corr$pca1=c1*corr$tau_e+c2*corr$tau_p
corr$pca2=-c1*corr$tau_e+c2*corr$tau_p

ggplot(agree, aes(x=pca1, y=pca2, color=factor(signal)),alpha=.01) +geom_point()
```

And using this rotated variables  the final mass model becomes

```{r,warning=FALSE, message=FALSE}
mask2=c("pca1","pca2")

param_m <- list(
  "objective" = "reg:linear",
  "max_depth"=7,
  "eta"=.1
)
data=train
y_d=train$mass

res=ftrain(data,y_d,param_m,folds,mask2,500,F,T)
p=as.data.frame(res[1])[,1]
magree=rowMeans(as.data.frame(res[2]))
mcorr=rowMeans(as.data.frame(res[3]))

train$mass_4=p
agree$mass_4=magree
corr$mass_4=mcorr

plot(train$mass,train$mass_4,pch='.',col=rgb(train$signal,0,0,.1),xlim=c(1600,1960),ylim=c(1600,1960))

ggplot(data=train, aes(mass_4)) +
  geom_histogram(data = subset(train,signal == 1),breaks=seq(1600, 2100, by =5), fill = "red", alpha = 0.3) +
  geom_histogram(data = subset(train,signal == 0),breaks=seq(1600, 2100, by =5),fill = "blue", alpha = 0.3) +
  geom_histogram(data = agree,breaks=seq(1600, 2100, by =5),fill = "green", alpha = 0.3)
```

Good. Now return to the classification problem

```{r,warning=FALSE, message=FALSE}
data=train
y_d=train$signal
mask2=c("pca1","pca2")
param <- list(
  "objective" = "binary:logistic",
  "eval_metric"="auc",
  "max_depth"=7,
  "eta"=.9
)

res=ftrain(data,y_d,param,folds,mask2,500,T)
p=as.data.frame(res[1])[,1]

mask=(data$min_ANNmuon>.4 )
ah=roc_auc_truncated(y_d[mask],p[mask])
print(ah)

```

It seems that gradient boosting  is  able to take profit from the mass information in a not supervised way.

And finally, introducing the calculated mass as a feature (stacking correctly the model: notice than the folds hasn't changed from the start) and adding all kinematic variables
```{r,warning=FALSE, message=FALSE}
data=train
y_d=train$signal
mask2=c("mass_4","pca1","pca2","p0_p","p1_p","p2_p","p0_pt","p1_pt","p2_pt","p0_eta","p1_eta","p2_eta","pt","LifeTime","FlightDistance")
param <- list(
  "objective" = "binary:logistic",
  "eval_metric"="auc",
  "max_depth"=7,
  "eta"=.1
)

res=ftrain(data,y_d,param,folds,mask2,500,T,T)
p=as.data.frame(res[1])[,1]
p3agree=rowMeans(as.data.frame(res[2]))
p3corr=rowMeans(as.data.frame(res[3]))

mask=(data$min_ANNmuon>.4 )
ah=roc_auc_truncated(y_d[mask],p[mask],)
print(ah)
```

Fairly good. But... what about agreement and correlation tests?


```{r} 
cc=compute_cvm(p3corr,corr$mass)
print(paste0("CVM=",cc))
  
d0=p3agree[agree$signal==0]
d1=p3agree[agree$signal==1]
w0=agree$weight[agree$signal==0]
w1=agree$weight[agree$signal==1]
kks=compute_ks(d0,d1,w0,w1)
print(paste0("K-S=",kks))
```

Too bad...    In the next chapter we'll see how to tame the mass_4 wildness, and the still more scary Curse of SPDhits...
But some kind of regularization will help.... (to be continued)