library(MASS)
library(ltm)
library(biclust)
source("./r/pca.lda.r")
source("./r/rasch.lda.r")
source("./r/my.binarize.r")

load("./data/leukemia.rda")

my.train <- as.matrix(read.table("./subsets/my.train.txt",h=F))

r50 <- as.matrix(read.table("./subsets/r50.txt",h=F))
pca.components.r50 <- read.table("./components/leukemia_pca_components_r50.txt",h=F)$V1
rasch.components.r50 <- read.table("./components/leukemia_rasch_components_r50.txt",h=F)$V1
r100 <- as.matrix(read.table("./subsets/r100.txt",h=F))
pca.components.r100 <- read.table("./components/leukemia_pca_components_r100.txt",h=F)$V1
rasch.components.r100 <- read.table("./components/leukemia_rasch_components_r100.txt",h=F)$V1
r200 <- as.matrix(read.table("./subsets/r200.txt",h=F))
pca.components.r200 <- read.table("./components/leukemia_pca_components_r200.txt",h=F)$V1
rasch.components.r200 <- read.table("./components/leukemia_rasch_components_r200.txt",h=F)$V1

pca.posterior.r50 <- NULL
pca.class.r50 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- pca.lda(Xtrain=lx.original[my.train[i,],r50[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],r50[i,]],ncomp=pca.components.r50[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    pca.posterior.r50 <- c(pca.posterior.r50,predict)
    pca.class.r50 <- c(pca.class.r50,class)
}
pca.roc.r50 <- roc(as.vector(pca.posterior.r50),as.factor(pca.class.r50))

pca.posterior.r100 <- NULL
pca.class.r100 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- pca.lda(Xtrain=lx.original[my.train[i,],r100[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],r100[i,]],ncomp=pca.components.r100[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    pca.posterior.r100 <- c(pca.posterior.r100,predict)
    pca.class.r100 <- c(pca.class.r100,class)
}
pca.roc.r100 <- roc(as.vector(pca.posterior.r100),as.factor(pca.class.r100))

pca.posterior.r200 <- NULL
pca.class.r200 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- pca.lda(Xtrain=lx.original[my.train[i,],r200[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],r200[i,]],ncomp=pca.components.r200[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    pca.posterior.r200 <- c(pca.posterior.r200,predict)
    pca.class.r200 <- c(pca.class.r200,class)   
}
pca.roc.r200 <- roc(as.vector(pca.posterior.r200),as.factor(pca.class.r200))

rasch.posterior.r50 <- NULL
rasch.class.r50 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- rasch.lda.predict(Xtrain=lx.original[my.train[i,],r50[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],r50[i,]],ncomp=rasch.components.r50[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    rasch.posterior.r50 <- c(rasch.posterior.r50,predict)
    rasch.class.r50 <- c(rasch.class.r50,class)
    cat("Complete: ",i,"\n")
}
rasch.roc.r50 <- roc(as.vector(rasch.posterior.r50),as.factor(rasch.class.r50))


rasch.posterior.r100 <- NULL
rasch.class.r100 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- rasch.lda.predict(Xtrain=lx.original[my.train[i,],r100[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],r100[i,]],ncomp=rasch.components.r100[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    rasch.posterior.r100 <- c(rasch.posterior.r100,predict)
    rasch.class.r100 <- c(rasch.class.r100,class)
    cat("Complete: ",i,"\n")
}
rasch.roc.r100 <- roc(as.vector(rasch.posterior.r100),as.factor(rasch.class.r100))

rasch.posterior.r200 <- NULL
rasch.class.r200 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- rasch.lda.predict(Xtrain=lx.original[my.train[i,],r200[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],r200[i,]],ncomp=rasch.components.r200[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    rasch.posterior.r200 <- c(rasch.posterior.r200,predict)
    rasch.class.r200 <- c(rasch.class.r200,class)
    cat("Complete: ",i,"\n")
}
rasch.roc.r200 <- roc(as.vector(rasch.posterior.r200),as.factor(rasch.class.r200))


p50 <- as.matrix(read.table("./subsets/p50.txt",h=F))
pca.components.p50 <- read.table("./components/leukemia_pca_components_p50.txt",h=F)$V1
rasch.components.p50 <- read.table("./components/leukemia_rasch_components_p50.txt",h=F)$V1
p100 <- as.matrix(read.table("./subsets/p100.txt",h=F))
pca.components.p100 <- read.table("./components/leukemia_pca_components_p100.txt",h=F)$V1
rasch.components.p100 <- read.table("./components/leukemia_rasch_components_p100.txt",h=F)$V1
p200 <- as.matrix(read.table("./subsets/p200.txt",h=F))
pca.components.p200 <- read.table("./components/leukemia_pca_components_p200.txt",h=F)$V1
rasch.components.p200 <- read.table("./components/leukemia_rasch_components_p200.txt",h=F)$V1



pca.posterior.p50 <- NULL
pca.class.p50 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- pca.lda(Xtrain=lx.original[my.train[i,],p50[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],p50[i,]],ncomp=pca.components.p50[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    pca.posterior.p50 <- c(pca.posterior.p50,predict)
    pca.class.p50 <- c(pca.class.p50,class)
}
pca.roc.p50 <- roc(as.vector(pca.posterior.p50),as.factor(pca.class.p50))

pca.posterior.p100 <- NULL
pca.class.p100 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- pca.lda(Xtrain=lx.original[my.train[i,],p100[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],p100[i,]],ncomp=pca.components.p100[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    pca.posterior.p100 <- c(pca.posterior.p100,predict)
    pca.class.p100 <- c(pca.class.p100,class)
}
pca.roc.p100 <- roc(as.vector(pca.posterior.p100),as.factor(pca.class.p100))

pca.posterior.p200 <- NULL
pca.class.p200 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- pca.lda(Xtrain=lx.original[my.train[i,],p200[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],p200[i,]],ncomp=pca.components.p200[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    pca.posterior.p200 <- c(pca.posterior.p200,predict)
    pca.class.p200 <- c(pca.class.p200,class)
}
pca.roc.p200 <- roc(as.vector(pca.posterior.p200),as.factor(pca.class.p200))
     
rasch.posterior.p50 <- NULL
rasch.class.p50 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- rasch.lda.predict(Xtrain=lx.original[my.train[i,],p50[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],p50[i,]],ncomp=rasch.components.p50[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    rasch.posterior.p50 <- c(rasch.posterior.p50,predict)
    rasch.class.p50 <- c(rasch.class.p50,class)
    cat("Complete: ",i,"\n")
}
rasch.roc.p50 <- roc(as.vector(rasch.posterior.p50),as.factor(rasch.class.p50))

rasch.posterior.p100 <- NULL
rasch.class.p100 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- rasch.lda.predict(Xtrain=lx.original[my.train[i,],p100[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],p100[i,]],ncomp=rasch.components.p100[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    rasch.posterior.p100 <- c(rasch.posterior.p100,predict)
    rasch.class.p100 <- c(rasch.class.p100,class)
    cat("Complete: ",i,"\n")
}
rasch.roc.p100 <- roc(as.vector(rasch.posterior.p100),as.factor(rasch.class.p100))

rasch.posterior.p200 <- NULL
rasch.class.p200 <- NULL
for(i in 1:nrow(my.train)) {
    predict <- rasch.lda.predict(Xtrain=lx.original[my.train[i,],p200[i,]],Ytrain=ly.original[my.train[i,]],Xtest=lx.original[-(my.train)[i,],p200[i,]],ncomp=rasch.components.p200[i])$posterior[,1]
    class <- ly.original[-(my.train)[i,]]
    rasch.posterior.p200 <- c(rasch.posterior.p200,predict)
    rasch.class.p200 <- c(rasch.class.p200,class)
    cat("Complete: ",i,"\n")
}
rasch.roc.p200 <- roc(as.vector(rasch.posterior.p200),as.factor(rasch.class.p200))


pdf("roc1_bw_temp.pdf",width=20,height=10,paper="special")
par(cex.axis=1.7,cex.lab=1.7,cex.main=1.7,pty="s",mfrow=c(1,2))
plot(1-pca.roc.r50[,3],pca.roc.r50[,2],type="l",cex=1.7,main="(a)",xlab=bquote(1-"Specificity"),ylab="Sensitivity",col="gray60",lty=1,lwd=3)
lines(1-pca.roc.r100[,3],pca.roc.r100[,2],col="gray60",lty=2,lwd=3)
lines(1-pca.roc.r200[,3],pca.roc.r200[,2],col="gray60",lty=3,lwd=3)
lines(1-rasch.roc.r50[,3],rasch.roc.r50[,2],col="black",lty=1,lwd=3)
lines(1-rasch.roc.r100[,3],rasch.roc.r100[,2],col="black",lty=2,lwd=3)
lines(1-rasch.roc.r200[,3],rasch.roc.r200[,2],col="black",lty=3,lwd=3)
legend("bottomright",c(expression(RM~(paste(p,"*",sep="")==50)),expression(RM~(paste(p,"*",sep="")==100)),expression(RM~(paste(p,"*",sep="")==200)),expression(PCA~(paste(p,"*",sep="")==50)),expression(PCA~(paste(p,"*",sep="")==100)),expression(PCA~(paste(p,"*",sep="")==200))),lty=c(1,2,3,1,2,3),lwd=c(3,3,3,3,3,3),col=c("black","black","black","gray60","gray60","gray60"),bty="n",cex=1.7)
plot(1-pca.roc.p50[,3],pca.roc.p50[,2],type="l",cex=1.7,main="(b)",xlab="1-Specificity",ylab="Sensitivity",col="gray60",lty=1,lwd=3)
lines(1-pca.roc.p100[,3],pca.roc.p100[,2],col="gray60",lty=2,lwd=3)
lines(1-pca.roc.p200[,3],pca.roc.p200[,2],col="gray60",lty=3,lwd=3)
lines(1-rasch.roc.p50[,3],rasch.roc.p50[,2],col="black",lty=1,lwd=3)
lines(1-rasch.roc.p100[,3],rasch.roc.p100[,2],col="black",lty=2,lwd=3)
lines(1-rasch.roc.p200[,3],rasch.roc.p200[,2],col="black",lty=3,lwd=3)
legend("bottomright",c(expression(RM~(paste(p,"*",sep="")==50)),expression(RM~(paste(p,"*",sep="")==100)),expression(RM~(paste(p,"*",sep="")==200)),expression(PCA~(paste(p,"*",sep="")==50)),expression(PCA~(paste(p,"*",sep="")==100)),expression(PCA~(paste(p,"*",sep="")==200))),lty=c(1,2,3,1,2,3),lwd=c(3,3,3,3,3,3),col=c("black","black","black","gray60","gray60","gray60"),bty="n",cex=1.7)
dev.off()

pdf("roc1_cl_temp.pdf",width=20,height=10,paper="special")
par(cex.axis=1.7,cex.lab=1.7,cex.main=1.7,pty="s",mfrow=c(1,2))
plot(1-pca.roc.r50[,3],pca.roc.r50[,2],type="l",cex=1.7,main="(a)",xlab=bquote(1-"Specificity"),ylab="Sensitivity",col="blue",lty=1,lwd=3)
lines(1-pca.roc.r100[,3],pca.roc.r100[,2],col="blue",lty=2,lwd=3)
lines(1-pca.roc.r200[,3],pca.roc.r200[,2],col="blue",lty=3,lwd=3)
lines(1-rasch.roc.r50[,3],rasch.roc.r50[,2],col="red",lty=1,lwd=3)
lines(1-rasch.roc.r100[,3],rasch.roc.r100[,2],col="red",lty=2,lwd=3)
lines(1-rasch.roc.r200[,3],rasch.roc.r200[,2],col="red",lty=3,lwd=3)
legend("bottomright",c(expression(RM~(paste(p,"*",sep="")==50)),expression(RM~(paste(p,"*",sep="")==100)),expression(RM~(paste(p,"*",sep="")==200)),expression(PCA~(paste(p,"*",sep="")==50)),expression(PCA~(paste(p,"*",sep="")==100)),expression(PCA~(paste(p,"*",sep="")==200))),lty=c(1,2,3,1,2,3),lwd=c(3,3,3,3,3,3),col=c("red","red","red","blue","blue","blue"),bty="n",cex=1.7)
plot(1-pca.roc.p50[,3],pca.roc.p50[,2],type="l",cex=1.7,main="(b)",xlab="1-Specificity",ylab="Sensitivity",col="blue",lty=1,lwd=3)
lines(1-pca.roc.p100[,3],pca.roc.p100[,2],col="blue",lty=2,lwd=3)
lines(1-pca.roc.p200[,3],pca.roc.p200[,2],col="blue",lty=3,lwd=3)
lines(1-rasch.roc.p50[,3],rasch.roc.p50[,2],col="red",lty=1,lwd=3)
lines(1-rasch.roc.p100[,3],rasch.roc.p100[,2],col="red",lty=2,lwd=3)
lines(1-rasch.roc.p200[,3],rasch.roc.p200[,2],col="red",lty=3,lwd=3)
legend("bottomright",c(expression(RM~(paste(p,"*",sep="")==50)),expression(RM~(paste(p,"*",sep="")==100)),expression(RM~(paste(p,"*",sep="")==200)),expression(PCA~(paste(p,"*",sep="")==50)),expression(PCA~(paste(p,"*",sep="")==100)),expression(PCA~(paste(p,"*",sep="")==200))),lty=c(1,2,3,1,2,3),lwd=c(3,3,3,3,3,3),col=c("red","red","red","blue","blue","blue"),bty="n",cex=1.7)
dev.off()



    






