require(glmnet)
require(sqldf)
require(ROCR)
require(ggplot2)

fulldata <- read.csv("propTableBPAvgBinaryClass.csv")
new <- fulldata
new$LoS <- ifelse(fulldata$LoS=="long",0,1)
new$ProThr <- fulldata$ProThrCHT227465 + fulldata$ProThrCHT227469

predictedPos = vector()
predictedNeg = vector()

gtPos = vector()
gtNeg = vector()

tn = tp = fp = fn = 0

for(i in 1:5) {
   i
   fileTRpos = paste0("folds/f",i,"/train_pos.R")
   fileTRneg = paste0("folds/f",i,"/train_neg.R")
   fileTepos = paste0("folds/f",i,"/test_pos.R")
   fileTeneg = paste0("folds/f",i,"/test_neg.R")

   tr_pos <- read.csv(fileTRpos)
   tr_neg <- read.csv(fileTRneg)
   te_pos <- read.csv(fileTepos)
   te_neg <- read.csv(fileTeneg)

   jo <- rbind(tr_pos,tr_neg)

   train <- sqldf("select * from new INNER JOIN jo ON new.Icu == jo.ICU")

   testPos <- sqldf("select * from new INNER JOIN te_pos ON new.Icu == te_pos.ICU")

   testNeg <- sqldf("select * from new INNER JOIN te_neg ON new.Icu == te_neg.ICU")

   gtPos <- c(gtPos,testPos$LoS)
   gtNeg <- c(gtNeg,testNeg$LoS)

   model <- glm(formula = LoS ~ Age + ProThr + 
                    Stool + ABPs, family = binomial(link = "logit"), data = train)

#   predPos <- predict(model,newdata=testPos,index=29,terms=pos, type="response")
#   predNeg <- predict(model,newdata=testNeg,index=29,terms=neg, type="response")
   predPos <- predict(model,newdata=testPos,index=29, type="response")
   predNeg <- predict(model,newdata=testNeg,index=29, type="response")

   predictedPos <- c(predictedPos,predPos)
   predictedNeg <- c(predictedNeg,predNeg)
   
   ## for(j in 1:length(predPos)) {
   ##     if (predPos[j] <= 10) { 
   ##         fn = fn + 1
   ##     }
   ##     else {
   ##         tp = tp + 1
   ##     }
   ## }
   ## for(j in 1:length(predNeg)) {
   ##     if (predNeg[j] <= 10) {
   ##         tn = tn + 1
   ##     }
   ##     else {
   ##         fp = fp + 1
   ##     }
   ## }
}

allpreds <- c(predictedPos,predictedNeg)
gts <- c(gtPos,gtNeg)
preds <- prediction(allpreds,gts)
perf <- performance(preds, measure = "tpr", x.measure = "fpr")
# p <- plot(perf, col=rainbow(10),ylim=c(0:1),xlim=c(0:1)) + abline(0,1)
auc <- performance(preds,"auc")
# now converting S4 class to vector
auc <- unlist(slot(auc, "y.values"))
# adding min and max ROC AUC to the center of the plot
minauc<-min(round(auc, digits = 2))
maxauc<-max(round(auc, digits = 2))
#minauct <- paste(c("min(AUC)  = "),minauc,sep="")
minauct <- paste(c("AUC = "),minauc,sep="")
maxauct <- paste(c("max(AUC) = "),maxauc,sep="")
#legend(0.3,0.6,c(minauct,maxauct,"\n"),border="white",cex=1.7,box.col = "white")
# legend(0.0,1.0,c(minauct,"\n"),border="white",box.col = "white")

tpr <- unlist(slot(perf,"y.values"))
fpr <- unlist(slot(perf,"x.values"))

detach("ROCR",unload=T)

bin = 0.1
diag = data.frame(x = seq(0, 1, by = bin), y = seq(0, 1, by = bin))

p <- qplot(x=fpr, y=tpr, xlab="1-Specificity", ylab="Sensitivity") + geom_line(linetype="solid", colour = "lightgray",size = 1)
p <- p + geom_line(data = diag, aes(x = x, y = y), size=0.3, colour = "red")

#) + geom_point(colour = "green") + geom_line(colour = "green") + geom_line(data = diag, aes(x = x, y = y), colour = "red")

# plot ILP points IMV with varying noise
# noise = 30
sf = data.frame(x=0.07, y=0.3)
p <- p + geom_point(data=sf,aes(x,y),colour="red",size=4)
# noise = 50
sf = data.frame(x=0.073, y=0.285)
p <- p + geom_point(data=sf,aes(x,y),colour="green",size=4)
# noise = 10
sf = data.frame(x=0.08, y=0.266)
p <- p + geom_point(data=sf,aes(x,y),colour="blue",size=4)

# plot ILP points for propTableBPAvgBinaryClass
# noise = 0
sf = data.frame(x=0.03,y=0.078)
p <- p + geom_point(data=sf,aes(x,y),colour="red",shape=24,fill="red",size=4) + scale_shape(solid=TRUE) # triangle shaped
# noise = 10
sf = data.frame(x=0.03,y=0.097)
p <- p + geom_point(data=sf,aes(x,y),colour="green",shape=24,fill="green",size=4) + scale_shape(solid=TRUE)
# noise = 500
sf = data.frame(x=0.034,y=0.084)
p <- p + geom_point(data=sf,aes(x,y),colour="blue",shape=24,fill="blue",size=4) + scale_shape(solid=TRUE)