Deep Learning v2
Gerhard Viljoen
3/15/2019
require(tufte)
#remove all objects from R session
rm(list=ls())
#define a function to concatenate all json objects into an R list object
wrangle <- function(run){
#load required opackages
require(jsonlite)
require(readtext)
#list all files in the run specified as the input parameter to this function
files <- list.files(path=paste0("/Users/gerhard/msc-thesis-data/000",run,"/JS"),
pattern="*json",
full.names=T,
recursive=TRUE)
#use the jsonlite::fromJSON function to convert JSON into list, starting at position 2, since the first JSON object is empty (only contains: "}")
j <- fromJSON(files[2])
#iterate through the rest of the files in the specified run
for(i in 3:length(files)){
#convert each into a list and concatenate them to the preceding list
f <- fromJSON(files[i])
j <- c(j,f)
}
j
}
#run the function defined above on the following runs
dat1 <- wrangle("265309")
dat2 <- wrangle("265377")
dat3 <- wrangle("265378")
#concat all json files
j <- c(dat1,dat2,dat3)
#free up memory
rm(dat1,dat2,dat3)
#extract p (in these earlier files, p is actually PT)
p <- sapply(j, `[[`,"momentum")
#get ADC signal
pads <- sapply(j, `[[`,"layer0")
#set up an empty vector to save the integrated charge deposit per pad
dedX <- numeric(length(pads))
#get particle Id
pid <- sapply(j, `[[`,"pdgCode")
#find out which pads have no data
nulls <- which(sapply(pads, is.null)==T)
#remove those indices from all lists we've created
pads[nulls] <- NULL
dedX <- dedX[-nulls]
p <- p[-nulls]
pid <- pid[-nulls]
#sometimes there are some weird empty lists within lists, not sure what causes them but let's remove them so we can move forward for now
weird.lists <- c()
#in this loop we find the datatype for each element in the list of pads
for(i in pads){
weird.lists <- c(weird.lists,typeof(i))
}
#we then find those entries which are not integer matrices
weird.lists <- which(weird.lists!="integer")
#remove them from all the lists we've created:
pads[weird.lists] <- NULL
dedX <- dedX[-weird.lists]
p <- p[-weird.lists]
pid <- pid[-weird.lists]
#integrated charge is the sum of ADC signal
for(i in 1:length(pads)){
dedX[i] <- sum(pads[[i]])
}
#assign particle ID based on PDG code
pid <- ifelse(pid==211|pid==-211,"pion","electron")
#plotme <- data.frame(cbind(p,dedX,pid))
#to plot the time evolution of the signal we construct an x-axis, which is simply the number of columns in the matrix of the ADC signal for each tracklet
w <- ncol(pads[[1]])
#we will fill each entry in the vector thus constructed, for each tracklet in our dataset
l <- length(pads)
#therefore for the time evolution of the signal, we construct a matrix with an integrated charge deposit signal in each time bin (number of columns), for each tracklet in our dataset (number of rows)
time.evolution <- matrix(ncol = w,nrow=l)
#take the column sums for each entry and be done with it
for(i in 1:nrow(time.evolution)){
time.evolution[i,] <- colSums(pads[[i]])
}
#remove the strange entries spoken about above
#time.evolution <- time.evolution[-c(nulls,weird.lists),]
Momentum Distribution
#display the 5 number summary statistics and quantiles for momenta of electrons and pions in our dataset
require(pander)
pander(summary(p[pid=="electron"]))
pander(summary(p[pid=="pion"]))
pander(quantile(p[pid=="electron"]))
pander(quantile(p[pid=="pion"]))
#plot histograms for momentum distributions for electrons and pions
par(mfrow=c(1,2))
hist(p[pid=="electron"],breaks=1000)
#add a vertical red line to show which particles will be cut from our data if we cut on p<=2
abline(v=2,col="red")
hist(p[pid=="pion"],breaks=1000)
abline(v=2,col="red")
#get the proportion of data we will lose by doing this:
length(p[p<=2])/length(p[p>2])
#rm(i,l,j,nulls,w,weird.lists)
#save.image("/Users/gerhard/msc-thesis-data/concept.rdata")
load("/Users/gerhard/msc-thesis-data/concept.rdata")
Feedforward NN for PID based on Time Evolution of the Signal
plot(x=1:ncol(time.evolution),y=time.evolution[1,],type="l")
mu <- mean(time.evolution)
time.evolution <- time.evolution/mu
k <- nrow(time.evolution)
dP.dt <- matrix(nrow = nrow(time.evolution),ncol=ncol(time.evolution)-1+
ncol(time.evolution)-2+ncol(time.evolution)-2+
ncol(time.evolution)-3)
for(i in 1:k){
dP.dt[i,] <- a <- c(diff(time.evolution[i,]),diff(diff(time.evolution[i,])),
diff(time.evolution[i,],2),diff(diff(time.evolution[i,],2)))
}
x <- cbind(time.evolution,dP.dt)
dedX <- scale(dedX)
x <- cbind(x,dedX)
#rm(pads,time.evolution,dP.dt,dedX,a,i,k,mu)
#save.image("/Users/gerhard/msc-thesis-data/concept2.rdata")
load("/Users/gerhard/msc-thesis-data/concept2.rdata")
library(keras)
use_condaenv('r-tensorflow')
#conda activate py3keras
pid <- ifelse(pid=="pion",0,1)
pid <- to_categorical(pid)
rm(model1)
model1 <- keras_model_sequential()
model1 %>%
layer_dense(units = 256,activation = "relu",input_shape = ncol(x)) %>%
layer_dropout(rate=0.2) %>%
layer_dense(units=256,activation = "relu") %>%
layer_dropout(rate=0.2) %>%
layer_dense(units=2,activation = "softmax")
model1 %>% compile(
optimizer=optimizer_sgd(),
loss="binary_crossentropy",
metrics=c("accuracy")
)
model1
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_1 (Dense) (None, 256) 29184
## ___________________________________________________________________________
## dropout_1 (Dropout) (None, 256) 0
## ___________________________________________________________________________
## dense_2 (Dense) (None, 256) 65792
## ___________________________________________________________________________
## dropout_2 (Dropout) (None, 256) 0
## ___________________________________________________________________________
## dense_3 (Dense) (None, 2) 514
## ===========================================================================
## Total params: 95,490
## Trainable params: 95,490
## Non-trainable params: 0
## ___________________________________________________________________________
model1 %>% fit(x,pid,epochs=20,validation_split=0.2)
model1 %>% save_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel1.h5")
load_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel1.h5")
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_1 (Dense) (None, 256) 29184
## ___________________________________________________________________________
## dropout_1 (Dropout) (None, 256) 0
## ___________________________________________________________________________
## dense_2 (Dense) (None, 256) 65792
## ___________________________________________________________________________
## dropout_2 (Dropout) (None, 256) 0
## ___________________________________________________________________________
## dense_3 (Dense) (None, 2) 514
## ===========================================================================
## Total params: 95,490
## Trainable params: 95,490
## Non-trainable params: 0
## ___________________________________________________________________________
rm(model2)
model2 <- keras_model_sequential()
model2%>%
layer_dense(units = 512,activation = "relu",input_shape = ncol(x)) %>%
layer_dropout(rate=0.2) %>%
layer_dense(units = 256,activation = "relu",input_shape = ncol(x)) %>%
layer_dropout(rate=0.2) %>%
layer_dense(units=128,activation = "relu") %>%
layer_dropout(rate=0.2) %>%
layer_dense(units=2,activation = "softmax")
model2 %>% compile(
optimizer=optimizer_sgd(),
loss="binary_crossentropy",
metrics=c("accuracy")
)
model2
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_4 (Dense) (None, 512) 58368
## ___________________________________________________________________________
## dropout_3 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_5 (Dense) (None, 256) 131328
## ___________________________________________________________________________
## dropout_4 (Dropout) (None, 256) 0
## ___________________________________________________________________________
## dense_6 (Dense) (None, 128) 32896
## ___________________________________________________________________________
## dropout_5 (Dropout) (None, 128) 0
## ___________________________________________________________________________
## dense_7 (Dense) (None, 2) 258
## ===========================================================================
## Total params: 222,850
## Trainable params: 222,850
## Non-trainable params: 0
## ___________________________________________________________________________
model2 %>% fit(x,pid,epochs=20,validation_split=0.2)
model2 %>% save_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel2.h5")
load_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel2.h5")
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_49 (Dense) (None, 512) 58368
## ___________________________________________________________________________
## dense_50 (Dense) (None, 256) 131328
## ___________________________________________________________________________
## dense_51 (Dense) (None, 128) 32896
## ___________________________________________________________________________
## dropout_33 (Dropout) (None, 128) 0
## ___________________________________________________________________________
## dense_52 (Dense) (None, 2) 258
## ===========================================================================
## Total params: 222,850
## Trainable params: 222,850
## Non-trainable params: 0
## ___________________________________________________________________________
p1 <- predict_proba(model1,x)
preds <- cbind(pid[,1],p1[,1])
rm(p1)
p1 <- predict_proba(model2,x)
preds <- cbind(preds,p1[,1])
rm(p1)
load("/Users/gerhard/msc-thesis-data/concept.rdata")
Only time evolution, no lag differentials
rm(model3)
model3 <- keras_model_sequential()
model3 %>%
layer_dense(units = 256,activation = "relu",input_shape = ncol(time.evolution)) %>%
layer_dropout(rate=0.2) %>%
layer_dense(units=256,activation = "relu") %>%
layer_dropout(rate=0.2) %>%
layer_dense(units=2,activation = "softmax")
model3 %>% compile(
optimizer=optimizer_sgd(),
loss="binary_crossentropy",
metrics=c("accuracy")
)
model3
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_8 (Dense) (None, 256) 6400
## ___________________________________________________________________________
## dropout_6 (Dropout) (None, 256) 0
## ___________________________________________________________________________
## dense_9 (Dense) (None, 256) 65792
## ___________________________________________________________________________
## dropout_7 (Dropout) (None, 256) 0
## ___________________________________________________________________________
## dense_10 (Dense) (None, 2) 514
## ===========================================================================
## Total params: 72,706
## Trainable params: 72,706
## Non-trainable params: 0
## ___________________________________________________________________________
pid <- ifelse(pid=="electron",1,0)
pid <- to_categorical(pid)
model3 %>% fit(time.evolution,pid,epochs=20,validation_split=0.2)
model3 %>% save_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel3.h5")
load_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel3.h5")
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_14 (Dense) (None, 256) 6400
## ___________________________________________________________________________
## dropout_10 (Dropout) (None, 256) 0
## ___________________________________________________________________________
## dense_15 (Dense) (None, 256) 65792
## ___________________________________________________________________________
## dropout_11 (Dropout) (None, 256) 0
## ___________________________________________________________________________
## dense_16 (Dense) (None, 2) 514
## ===========================================================================
## Total params: 72,706
## Trainable params: 72,706
## Non-trainable params: 0
## ___________________________________________________________________________
p1 <- predict_proba(model3,time.evolution)
preds <- cbind(preds,p1[,1])
rm(p1)
Back to derived features
rm(model4)
model4 <- keras_model_sequential()
model4 %>%
layer_dense(units = 512,activation = "relu",input_shape = ncol(x),
kernel_regularizer = regularizer_l1_l2(l1 = 0.001,l2=0.001)) %>%
layer_dropout(rate=0.6) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dropout(rate=0.6) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l1_l2(l1 = 0.001,l2=0.001)) %>%
layer_dropout(rate=0.6) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l1_l2(l1 = 0.001,l2=0.001)) %>%
layer_dropout(rate=0.6) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dropout(rate=0.6) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dropout(rate=0.4) %>%
layer_dense(units=2,activation = "softmax")
model4 %>% compile(
optimizer=optimizer_sgd(lr=0.001),
loss="binary_crossentropy",
metrics=c("accuracy")
)
model4
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_11 (Dense) (None, 512) 58368
## ___________________________________________________________________________
## dropout_8 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_12 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_9 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_13 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_10 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_14 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_11 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_15 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_12 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_16 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_13 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_17 (Dense) (None, 2) 1026
## ===========================================================================
## Total params: 1,372,674
## Trainable params: 1,372,674
## Non-trainable params: 0
## ___________________________________________________________________________
rm(j,pads,time.evolution,dedX)
preds <- cbind(preds,p)
preds <- as.data.frame(preds)
names(preds) <- c("actual","p1","p2","p3","momentum")
model4 %>% fit(x,pid,epochs=50,validation_split=0.2, batch_size=128)
model4 %>% save_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel4.h5")
rm(model4.1)
model4.1 <- keras_model_sequential()
model4.1 %>%
layer_dense(units = 512,activation = "relu",input_shape = ncol(x),
kernel_regularizer = regularizer_l1_l2(l1 = 0.001,l2=0.001)) %>%
layer_dropout(rate=0.3) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dropout(rate=0.3) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l1_l2(l1 = 0.001,l2=0.001)) %>%
layer_dropout(rate=0.3) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l1_l2(l1 = 0.001,l2=0.001)) %>%
layer_dropout(rate=0.3) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l1_l2(l1 = 0.001,l2=0.001)) %>%
layer_dropout(rate=0.3) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l1_l2(l1 = 0.001,l2=0.001)) %>%
layer_dropout(rate=0.4) %>%
layer_dense(units=2,activation = "softmax")
model4.1 %>% compile(
optimizer=optimizer_adam(),
loss="binary_crossentropy",
metrics=c("accuracy")
)
model4.1
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_18 (Dense) (None, 512) 58368
## ___________________________________________________________________________
## dropout_14 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_19 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_15 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_20 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_16 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_21 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_17 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_22 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_18 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_23 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_19 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_24 (Dense) (None, 2) 1026
## ===========================================================================
## Total params: 1,372,674
## Trainable params: 1,372,674
## Non-trainable params: 0
## ___________________________________________________________________________
model4.1 %>% fit(x,pid,epochs=50,validation_split=0.3, batch_size=512)
model4.1 %>% save_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel4_1.h5")
load_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel4_1.h5")
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_31 (Dense) (None, 512) 58368
## ___________________________________________________________________________
## dropout_24 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_32 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_25 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_33 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_26 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_34 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_27 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_35 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_28 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_36 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_29 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_37 (Dense) (None, 2) 1026
## ===========================================================================
## Total params: 1,372,674
## Trainable params: 1,372,674
## Non-trainable params: 0
## ___________________________________________________________________________
rm(model4.2)
model4.2 <- keras_model_sequential()
model4.2 %>%
layer_dense(units = 512,activation = "relu",input_shape = ncol(x),
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dense(units=512,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dense(units=256,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dense(units=256,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dense(units=256,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dense(units=2,activation = "softmax")
model4.2 %>% compile(
optimizer=optimizer_adam(),
loss="binary_crossentropy",
metrics=c("accuracy")
)
model4.2
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_25 (Dense) (None, 512) 58368
## ___________________________________________________________________________
## dense_26 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dense_27 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dense_28 (Dense) (None, 256) 131328
## ___________________________________________________________________________
## dense_29 (Dense) (None, 256) 65792
## ___________________________________________________________________________
## dense_30 (Dense) (None, 256) 65792
## ___________________________________________________________________________
## dense_31 (Dense) (None, 2) 514
## ===========================================================================
## Total params: 847,106
## Trainable params: 847,106
## Non-trainable params: 0
## ___________________________________________________________________________
model4.2 %>% fit(x,pid,epochs=50,validation_split=0.3, batch_size=512)
model4.2 %>% save_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel4_2.h5")
load_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel4_2.h5")
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_38 (Dense) (None, 512) 58368
## ___________________________________________________________________________
## dense_39 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dense_40 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dense_41 (Dense) (None, 256) 131328
## ___________________________________________________________________________
## dense_42 (Dense) (None, 256) 65792
## ___________________________________________________________________________
## dense_43 (Dense) (None, 256) 65792
## ___________________________________________________________________________
## dense_44 (Dense) (None, 2) 514
## ===========================================================================
## Total params: 847,106
## Trainable params: 847,106
## Non-trainable params: 0
## ___________________________________________________________________________
rm(model2.1)
model2.1 <- keras_model_sequential()
model2.1%>%
layer_dense(units = 512,activation = "relu",input_shape = ncol(x),
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dense(units = 256,activation = "relu",input_shape = ncol(x),
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dense(units=128,activation = "relu") %>%
layer_dropout(rate=0.3) %>%
layer_dense(units=2,activation = "softmax")
model2.1 %>% compile(
optimizer=optimizer_adamax(),
loss="binary_crossentropy",
metrics=c("accuracy","categorical_accuracy","mae")
)
model2
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_4 (Dense) (None, 512) 58368
## ___________________________________________________________________________
## dropout_3 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_5 (Dense) (None, 256) 131328
## ___________________________________________________________________________
## dropout_4 (Dropout) (None, 256) 0
## ___________________________________________________________________________
## dense_6 (Dense) (None, 128) 32896
## ___________________________________________________________________________
## dropout_5 (Dropout) (None, 128) 0
## ___________________________________________________________________________
## dense_7 (Dense) (None, 2) 258
## ===========================================================================
## Total params: 222,850
## Trainable params: 222,850
## Non-trainable params: 0
## ___________________________________________________________________________
model2.1 %>% fit(x,pid,epochs=50,validation_split=0.3)
model2.1 %>% save_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodel2.h5")
New features:
p1 <- predict_proba(model4.2,x)
preds <- cbind(preds,p1[,1])
rm(p1)
p1 <- predict_proba(model2.1,x)
preds <- cbind(preds,p1[,1])
rm(p1)
preds <- as.data.frame(preds)
names(preds)[c(6,7)] <- c("p4.2","p2.1")
preds <- preds[,c(1:4,6,7,5)]
table(preds$actual)
##
## 0 1
## 20777 174448
20777/174448
## [1] 0.1191014
load("/Users/gerhard/msc-thesis-data/concept.rdata")
threesum1 <- time.evolution[,1]+time.evolution[,2]+time.evolution[,3]
threesum2 <- time.evolution[,4]+time.evolution[,5]+time.evolution[,6]
threesum3 <- time.evolution[,7]+time.evolution[,7]+time.evolution[,9]
threesum4 <- time.evolution[,10]+time.evolution[,11]+time.evolution[,12]
threesum5 <- time.evolution[,13]+time.evolution[,14]+time.evolution[,15]
threesum6 <- time.evolution[,16]+time.evolution[,17]+time.evolution[,18]
threesum7 <- time.evolution[,19]+time.evolution[,20]+time.evolution[,21]
threesum8 <- time.evolution[,22]+time.evolution[,23]+time.evolution[,24]
threesum <- cbind(threesum1,threesum2,
threesum3,threesum4,threesum5,threesum6,threesum7,
threesum8)
rm(pads,threesum1,threesum2,threesum3,threesum4,threesum5,threesum6,threesum7,threesum8,dedX)
rm(p,x,model1,model2,model2.1,model3,model4,model4.1,model4.2)
mu <- mean(threesum)
threesum <- threesum/mu
k <- nrow(threesum)
dP.dt <- matrix(nrow = nrow(threesum),ncol=ncol(threesum)-1+ncol(threesum)-2)
for(i in 1:k){
dP.dt[i,] <- a <- c(diff(threesum[i,]),diff(diff(threesum[i,])))
}
x_2 <- cbind(threesum,dP.dt)
pid <- ifelse(pid=="electron",1,0)
pid <- to_categorical(pid)
rm(model.z)
model.z <- keras_model_sequential()
model.z%>%
layer_dense(units = 512,activation = "relu",input_shape = ncol(x_2)) %>%
layer_dropout(rate=0.35) %>%
layer_dense(units = 512,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dense(units = 256,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dense(units=128,activation = "relu",
kernel_regularizer = regularizer_l2(l = 0.001)) %>%
layer_dense(units=2,activation = "softmax")
model.z %>% compile(
optimizer=optimizer_adamax(),
loss="binary_crossentropy",
metrics=c("accuracy","categorical_accuracy","mae")
)
model.z
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_36 (Dense) (None, 512) 11264
## ___________________________________________________________________________
## dropout_21 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_37 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dense_38 (Dense) (None, 256) 131328
## ___________________________________________________________________________
## dense_39 (Dense) (None, 128) 32896
## ___________________________________________________________________________
## dense_40 (Dense) (None, 2) 258
## ===========================================================================
## Total params: 438,402
## Trainable params: 438,402
## Non-trainable params: 0
## ___________________________________________________________________________
model.z %>% fit(x_2,pid,epochs=25,validation_split=0.3)
model.z %>% save_model_hdf5("/Users/gerhard/Msc-thesis/NEW/Rmodelz.h5")
p1 <- predict_proba(model.z,x_2)
preds <- cbind(preds[,1:6],p1[,1],preds[,7])
names(preds)[c(7,8)] <- c("pz","momentum")
rm(dP.dt,threesum,x_2,a,k,i,model.z,mu)
rm(p1)
rs <- numeric(nrow(time.evolution))
for(i in 1:nrow(time.evolution)){
rs[i] <- sum(time.evolution[i,])
}
zeros <- which(rs==0)
pid <- pid[-zeros,]
preds <- preds[-zeros,]
time.evolution <- time.evolution[-zeros]
rm(i, zeros)
pred.f <- as.matrix(preds[,c(2:7)])
p.f <- numeric(nrow(pred.f))
for(i in 1:length(p.f)){
p.f[i] <- mean(as.numeric(unlist(pred.f[i,])))
}
p.a <- data.frame(cbind(preds$actual,p.f))
p.a$err <- p.a$V1-p.a$p.f
hist(p.a$err)
hist(p.a$p.f,breaks=100)
p.a.low.p <- which(preds$momentum<=2)
preds <- preds[p.a.low.p,]
p.a <- p.a[p.a.low.p,]
final.prediction <- ifelse(p.a$p.f>.5,1,0)
p.a <- data.frame(cbind(p.a,final.prediction))
pions <- which(p.a$V1==1)
electrons <- which(p.a$V1==0)
pions <- p.a[pions,]
electrons <- p.a[electrons,]
pions <- pions[,c(1,4)]
electrons <- electrons[,c(1,4)]
electrons$correct <- electrons[,1]==electrons[,2]
pions$correct <- pions[,1]==pions[,2]
p.top <- length(which(pions$correct==T))
p.bottom <- length(pions$correct)
p.accuracy <- p.top/p.bottom
e.top <- length(which(electrons$correct==T))
e.bottom <- length(electrons$correct)
e.accuracy <- e.top/e.bottom
e.accuracy
## [1] 0.680606
p.accuracy
## [1] 0.368967
