<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Install important packages for boosted decision trees
Step1:  Importing packages and defining functions and variables
Step2:  Defining autoencoder model, Training and evaluation functions
Step3:  Mount google drive to access data
Get the data from the link [https
Step4:  Check the files exist and make a copy of the autoencoder model for backup
Step5:  CD to the main data directory and load the trained model
Step6:  Train another round if required
(You are strongly advised not to do this but to just use the trained model...)
Step7:  Evaluation using the trained model
Step8:  Read the important data
Here, we read the autoencoder loss, jet mass and nsubjettiness (mass and nsubjettiness were evaluated using the C programs main.cc, main.hh and all.hh found in [https
Step9:  Plotting and checking
Plot $\epsilon$ (autoencoder loss)
Step10:  Plot $m_J$ (jet mass)
Step11:  Plot jet $\tau_1$ (nsubjettiness)
Step12:  Plot jet $\tau_2$ (nsubjettiness)
Step13:  Plot jet $\tau_3$ (nsubjettiness)
Step14:  Plot jet $\tau_4$ (nsubjettiness)
Step15:  Plot ROC using only $\epsilon$
X axis is $\epsilon_t$
Y axis is $\frac{1}{\epsilon_{\text{QCD}}}$
Step16:  Combining variables
We now repeat the ROC exercise using boosted decision trees to combine the variables.
Step17:  Decision trees using only autoencoder loss
Step18:  Plot the ROC from the above model
Step19:  Train BDT using all variables
Step20:  Plot ROC using above model
Step21:  Not using the autoencoder loss
Step22:  Plot ROC for above model
<ASSISTANT_TASK:>
Python Code:
!pip install lightgbm
!pip install shap
Explanation: Install important packages for boosted decision trees
End of explanation
%tensorflow_version 1.x
import lzma
from google.colab import drive
import numpy as np
import tensorflow as tf
import keras
from keras import backend as K
from keras.layers import Input, Dense
from keras.models import Model
import matplotlib.pyplot as plt
import lightgbm as lgb#t
import shap
import sklearn
from sklearn import svm
from sklearn import preprocessing
from sklearn import datasets
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import AdaBoostRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.metrics import precision_recall_curve
# from sklearn.metrics import plot_precision_recall_curve
from sklearn.metrics import average_precision_score
from sklearn.metrics import roc_curve, auc
from sklearn.preprocessing import label_binarize
from sklearn.multiclass import OneVsRestClassifier
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix
from sklearn.decomposition import PCA
from sklearn.svm import SVC
#from scipy import interp
from sklearn.metrics import roc_auc_score
def READ_XZ (filename):
    file = lzma.LZMAFile(filename)
    type_bytes = file.read(-1)
    type_array = np.frombuffer(type_bytes,dtype='float32')                                                
    return type_array
def Count(array,val):
  count = 0.0
  for e in range(array.shape[0]):
    if array[e]>val :
      count=count+1.0
  return count / array.shape[0]
width=40
batch_size=200
ModelName = "Model_40_24_8_24_40_40"
config = tf.ConfigProto( device_count = {'GPU': 1 , 'CPU': 2} ) 
sess = tf.Session(config=config)
keras.backend.set_session(sess)
K.tensorflow_backend._get_available_gpus()
Explanation: Importing packages and defining functions and variables
End of explanation
# this is our input placeholder
input_img = Input(shape=(width*width,))
# "encoded" is the encoded representation of the input
Layer1 = Dense(24*24, activation='relu')(input_img)
Layer2 = Dense(8*8, activation='relu')(Layer1)
Layer3 = Dense(24*24, activation='relu')(Layer2)
Layer4 = Dense(40*40, activation='relu')(Layer3)
Out = Dense(40*40, activation='softmax')(Layer4)
# this model maps an input to its reconstruction
autoencoder = Model(input_img, Out)
autoencoder.compile(optimizer='adam', loss='mean_squared_error')
def NAME(eventtype,purpose,i,obs) :
  return "./"+eventtype+"/"+purpose+"/"+obs+"."+str(i)+".bin.xz"
#
def EvalOnFile (InFileName,OutFileName):
  data = READ_XZ (InFileName)
  x_train = data.reshape(-1,width*width)
  x_out = autoencoder.predict(x_train,200,use_multiprocessing=True)
  diff = x_train - x_out
  lrnorm = np.ones((diff.shape[0]))
  for e in range(diff.shape[0]):
    lrnorm[e] = np.linalg.norm(diff[e])
  lrnorm.tofile(OutFileName)
  print(lrnorm.shape)
BATCH_SIZE=512
def TrainOnFile (filename,testfilename,totalepochs):
  data = READ_XZ (filename)
  x_train = data.reshape(-1,width*width)
  datatest = READ_XZ (testfilename)
  x_test = datatest.reshape(-1,width*width)
  autoencoder.fit(
      x_train, x_train, epochs=totalepochs,
      batch_size=BATCH_SIZE, shuffle=True,
      validation_data=(x_test, x_test)
  )
  autoencoder.save(ModelName)
Explanation: Defining autoencoder model, Training and evaluation functions
End of explanation
drive.mount('/gdrive')
%cd /gdrive
Explanation: Mount google drive to access data
Get the data from the link [https://drive.google.com/drive/folders/1_voPoiETqfWmCmBeUCKiF5oXqzbr-ZWt?usp=sharing] (this link is public and anyone can download the image data and trained models)
End of explanation
%cd /gdrive/My Drive/S2
!ls
!cp ./Model_40_24_8_24_40_40 ../Model_40_24_8_24_40_40.bak
!ls ../
# !tar -xf S2.tar
Explanation: Check the files exist and make a copy of the autoencoder model for backup
End of explanation
%cd /gdrive/My Drive/S2
autoencoder = keras.models.load_model(ModelName)
Explanation: CD to the main data directory and load the trained model
End of explanation
%cd /gdrive/My Drive/S2
#autoencoder = keras.models.load_model(ModelName)
#!ls ./TOP/TRAIN/*.*.bin.xz
#
for e in range(20):
  for i in range(7):
    TrainOnFile(NAME("QCD","TRAIN",i%7,"out"),NAME("QCD","TEST",i%3,"out"),10)
  #
  for i in range(3):
    TrainOnFile(NAME("QCD","VAL",i%7,"out"),NAME("QCD","TEST",i%3,"out"),10)
  #
#
Explanation: Train another round if required
(You are strongly advised not to do this but to just use the trained model...)
End of explanation
for i in range(7) :
  EvalOnFile(NAME("QCD","TRAIN",i,"out"),NAME("QCD","TRAIN",i,"loss"))
  EvalOnFile(NAME("TOP","TRAIN",i,"out"),NAME("TOP","TRAIN",i,"loss"))
#
for i in range(3) :
  EvalOnFile(NAME("QCD","TEST",i,"out"),NAME("QCD","TEST",i,"loss"))
  EvalOnFile(NAME("TOP","TEST",i,"out"),NAME("TOP","TEST",i,"loss"))
  EvalOnFile(NAME("QCD","VAL",i,"out"),NAME("QCD","VAL",i,"loss"))
  EvalOnFile(NAME("TOP","VAL",i,"out"),NAME("TOP","VAL",i,"loss"))
#
Explanation: Evaluation using the trained model
End of explanation
def ReadLossMassNsub(eventtype,sampletype,i):
  loss = np.fromfile(NAME(eventtype,sampletype,i,"loss"), dtype=float)
  mass = READ_XZ(NAME(eventtype,sampletype,i,"mass"))
  nsub = READ_XZ(NAME(eventtype,sampletype,i,"nsub")).reshape(-1,5)
  #print(nsub.shape)
  out = np.ones((mass.shape[0],7))
  for i in range(mass.shape[0]):
    out[i][0] = loss[i]
    out[i][1] = mass[i]
    out[i][2] = nsub[i][0]
    out[i][3] = nsub[i][1]
    out[i][4] = nsub[i][2]
    out[i][5] = nsub[i][3]
    out[i][6] = nsub[i][4]
  #
  return out
#
vars_qcd_train = ReadLossMassNsub("QCD","TRAIN",0)
vars_qcd_train = np.append (vars_qcd_train,ReadLossMassNsub("QCD","TRAIN",1),0)
vars_qcd_train = np.append (vars_qcd_train,ReadLossMassNsub("QCD","TRAIN",2),0)
vars_qcd_train = np.append (vars_qcd_train,ReadLossMassNsub("QCD","TRAIN",3),0)
vars_qcd_train = np.append (vars_qcd_train,ReadLossMassNsub("QCD","TRAIN",4),0)
vars_qcd_train = np.append (vars_qcd_train,ReadLossMassNsub("QCD","TRAIN",5),0)
vars_qcd_train = np.append (vars_qcd_train,ReadLossMassNsub("QCD","TRAIN",6),0)
vars_qcd_test = ReadLossMassNsub("QCD","TEST",0)
vars_qcd_test = np.append (vars_qcd_test,ReadLossMassNsub("QCD","TEST",1),0)
vars_qcd_test = np.append (vars_qcd_test,ReadLossMassNsub("QCD","TEST",2),0)
vars_qcd_val = ReadLossMassNsub("QCD","VAL",0)
vars_qcd_val = np.append (vars_qcd_val,ReadLossMassNsub("QCD","VAL",1),0)
vars_qcd_val = np.append (vars_qcd_val,ReadLossMassNsub("QCD","VAL",2),0)
vars_top_train = ReadLossMassNsub("TOP","TRAIN",0)
vars_top_train = np.append (vars_top_train,ReadLossMassNsub("TOP","TRAIN",1),0)
vars_top_train = np.append (vars_top_train,ReadLossMassNsub("TOP","TRAIN",2),0)
vars_top_train = np.append (vars_top_train,ReadLossMassNsub("TOP","TRAIN",3),0)
vars_top_train = np.append (vars_top_train,ReadLossMassNsub("TOP","TRAIN",4),0)
vars_top_train = np.append (vars_top_train,ReadLossMassNsub("TOP","TRAIN",5),0)
vars_top_train = np.append (vars_top_train,ReadLossMassNsub("TOP","TRAIN",6),0)
vars_top_test = ReadLossMassNsub("TOP","TEST",0)
vars_top_test = np.append (vars_top_test,ReadLossMassNsub("TOP","TEST",1),0)
vars_top_test = np.append (vars_top_test,ReadLossMassNsub("TOP","TEST",2),0)
vars_top_val = ReadLossMassNsub("TOP","VAL",0)
vars_top_val = np.append (vars_top_val,ReadLossMassNsub("TOP","VAL",1),0)
vars_top_val = np.append (vars_top_val,ReadLossMassNsub("TOP","VAL",2),0)
Explanation: Read the important data
Here, we read the autoencoder loss, jet mass and nsubjettiness (mass and nsubjettiness were evaluated using the C programs main.cc, main.hh and all.hh found in [https://github.com/aravindhv10/CPP_Wrappers/tree/master/AntiQCD4/ML4JETS] )
End of explanation
plt.hist(vars_qcd_test[:,0],100,(0.0,0.4),density=True,histtype='step')
plt.hist(vars_top_test[:,0],100,(0.0,0.4),density=True,histtype='step')
plt.show()
Explanation: Plotting and checking
Plot $\epsilon$ (autoencoder loss)
End of explanation
plt.hist(vars_qcd_test[:,1],100,(0.0,1000),density=True,histtype='step')
plt.hist(vars_top_test[:,1],100,(0.0,1000),density=True,histtype='step')
plt.show()
Explanation: Plot $m_J$ (jet mass)
End of explanation
plt.hist(vars_qcd_test[:,2],100,(0.0,100),density=True,histtype='step')
plt.hist(vars_top_test[:,2],100,(0.0,100),density=True,histtype='step')
plt.show()
Explanation: Plot jet $\tau_1$ (nsubjettiness)
End of explanation
plt.hist(vars_qcd_test[:,3],100,(0.0,100),density=True,histtype='step')
plt.hist(vars_top_test[:,3],100,(0.0,100),density=True,histtype='step')
plt.show()
Explanation: Plot jet $\tau_2$ (nsubjettiness)
End of explanation
plt.hist(vars_qcd_test[:,4],100,(0.0,100),density=True,histtype='step')
plt.hist(vars_top_test[:,4],100,(0.0,100),density=True,histtype='step')
plt.show()
Explanation: Plot jet $\tau_3$ (nsubjettiness)
End of explanation
plt.hist(vars_qcd_test[:,5],100,(0.0,100),density=True,histtype='step')
plt.hist(vars_top_test[:,5],100,(0.0,100),density=True,histtype='step')
plt.show()
Explanation: Plot jet $\tau_4$ (nsubjettiness)
End of explanation
dx = (0.4 - 0.0) / 100.0
qcdeff = np.ones((100))
topeff = np.ones((100))
for i in range(100):
  xval = i*dx
  qcdeff[i]=1.0/(Count(vars_qcd_test[:,0],xval)+0.0000000001)
  topeff[i]=Count(vars_top_test[:,0],xval)
plt.yscale('log')
plt.plot(topeff,qcdeff)
Explanation: Plot ROC using only $\epsilon$
X axis is $\epsilon_t$
Y axis is $\frac{1}{\epsilon_{\text{QCD}}}$
End of explanation
import sklearn
def prepare (qcd_vars,top_vars) :
  out_x = np.append(qcd_vars,top_vars,0)
  out_y = np.append(np.zeros((qcd_vars.shape[0]),dtype='float32'),np.ones((top_vars.shape[0]),dtype='float32'),0)
  return sklearn.utils.shuffle ( out_x , out_y , random_state=0 )
train_x, train_y = prepare(vars_qcd_train,vars_top_train)
test_x, test_y = prepare(vars_qcd_test,vars_top_test)
val_x, val_y = prepare(vars_qcd_val,vars_top_val)
param = { 'objective':'binary' , 'metric':'auc,binary_logloss,binary_error' }
plt.hist(train_x[:,0],100,(0.0,0.4),density=True,histtype='step')
plt.hist(test_x[:,0],100,(0.0,0.4),density=True,histtype='step')
plt.show()
plt.hist(train_x[:,1],100,(0.0,1000),density=True,histtype='step')
plt.hist(test_x[:,1],100,(0.0,1000),density=True,histtype='step')
plt.show()
Explanation: Combining variables
We now repeat the ROC exercise using boosted decision trees to combine the variables.
End of explanation
num_round = 100
#train_data = lgb.Dataset( train_x[:,0:0] , label=train_y )
#val_data   = lgb.Dataset( val_x[:,0:0]   , label=val_y   )
train_data = lgb.Dataset( train_x[:,0].reshape((-1,1)) , label=train_y )
val_data   = lgb.Dataset( val_x[:,0].reshape((-1,1))   , label=val_y   )
bst = lgb.train(param, train_data, num_round, valid_sets=val_data)
pred_qcd_test = bst.predict(vars_qcd_test[:,0].reshape((-1,1)))
pred_top_test = bst.predict(vars_top_test[:,0].reshape((-1,1)))
Explanation: Decision trees using only autoencoder loss
End of explanation
epsilon = 0.0000001
num = 1000
dx = ( 1.0 + (epsilon*2) ) / num
qcdeff_loss = np.ones((num))
topeff_loss = np.ones((num))
for i in range(num):
  xval = (i*dx) - epsilon
  qcdeff_loss[i]=1.0/(Count(pred_qcd_test,xval)+epsilon)
  topeff_loss[i]=Count(pred_top_test,xval)
plt.yscale('log')
plt.plot(topeff_loss,qcdeff_loss)
Explanation: Plot the ROC from the above model
End of explanation
num_round = 100
train_data = lgb.Dataset( train_x[:,0:6] , label=train_y )
val_data   = lgb.Dataset( val_x[:,0:6]   , label=val_y   )
bst = lgb.train(param, train_data, num_round, valid_sets=val_data)
pred_qcd_test = bst.predict(vars_qcd_test[:,0:6])
pred_top_test = bst.predict(vars_top_test[:,0:6])
Explanation: Train BDT using all variables
End of explanation
epsilon = 0.0000001
num = 1000
dx = ( 1.0 + (epsilon*2) ) / num
qcdeff_all = np.ones((num))
topeff_all = np.ones((num))
for i in range(num):
  xval = (i*dx) - epsilon
  qcdeff_all[i]=1.0/(Count(pred_qcd_test,xval)+epsilon)
  topeff_all[i]=Count(pred_top_test,xval)
plt.yscale('log')
plt.plot(topeff_all,qcdeff_all)
Explanation: Plot ROC using above model
End of explanation
num_round = 100
train_data = lgb.Dataset( train_x[:,1:6] , label=train_y )
val_data   = lgb.Dataset( val_x[:,1:6]   , label=val_y   )
bst = lgb.train(param, train_data, num_round, valid_sets=val_data)
pred_qcd_test = bst.predict(vars_qcd_test[:,1:6])
pred_top_test = bst.predict(vars_top_test[:,1:6])
Explanation: Not using the autoencoder loss
End of explanation
epsilon = 0.0000001
num = 1000
dx = ( 1.0 + (epsilon*2) ) / num
qcdeff_noloss = np.ones((num))
topeff_noloss = np.ones((num))
for i in range(num):
  xval = (i*dx) - epsilon
  qcdeff_noloss[i]=1.0/(Count(pred_qcd_test,xval)+epsilon)
  topeff_noloss[i]=Count(pred_top_test,xval)
plt.yscale('log')
plt.plot(topeff_noloss,qcdeff_noloss)
np.savetxt("topeff_loss",topeff_loss)
np.savetxt("qcdeff_loss",qcdeff_loss)
np.savetxt("topeff_all",topeff_all)
np.savetxt("qcdeff_all",qcdeff_all)
np.savetxt("topeff_noloss",topeff_noloss)
np.savetxt("qcdeff_noloss",qcdeff_noloss)
Explanation: Plot ROC for above model
End of explanation
<END_TASK>