prefix <- "tt"
first <- 0
last <-  0

skip <-  1
num  <- 10

alag <- 100
  
library(MASS)
library(ts)
  
ps.options(horizontal=FALSE,onefile=FALSE)
ps.options(pagecentre=TRUE)
ps.options(paper="letter")
ps.options(width=7.0,height=6.0)

for (i in (first:last)) {
  if (i < 10) {
    rootname <- paste(prefix,".sci_sub_logl.","00",i,sep="")
  } else {
    rootname <- paste(prefix,".sci_sub_logl.","0",i,sep="")
  }

  filename <- paste(rootname,".dat",sep="")

  tmp  <- scan(filename)
  lsub <- tmp[1]
  cols <- tmp[2]
  
  logli <- tmp[3:(num*lsub+2)] 
  idxi <- seq(1,num*lsub,skip)
  
  top <- max(logli)
  bot <- min(logli)
  
  filename <- paste(rootname,".eps",sep="")
  postscript(file=filename)
  
  plot(idxi,logli[idxi],ylab="",xlab="")
  for (j in 0:num) {
    lines(c(lsub*j,lsub*j),c(bot,top))
  }
  
  dev.off()

  if (i == first) {
    logl <- logli[idxi]
    idx  <- idxi
  } else {
    logl <- c(logl,logli[idxi])
    idx  <- c(idx,idxi+(i-first)*num*lsub)
  }
  
  if (i < 10) {
   filename <- paste(prefix,".sci_sub_parm.","00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".sci_sub_parm.","0",i,".dat",sep="")
  }

  tmp <- scan(filename)
  rows <- tmp[1]
  cols <- tmp[2]
  etai <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
  if (i == first) {
    eta <- etai[,idxi]
  } else {
    eta <- cbind(eta,etai[,idxi])
  }

  if (i < 10) {
   filename <- paste(prefix,".sci_sub_rej.","00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".sci_sub_rej.","0",i,".dat",sep="")
  }

  tmp <- scan(filename)
  rows <- tmp[1]
  cols <- tmp[2]
  rej <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
  if (i == first) {
    sum <- mat.or.vec(rows,1);
    count <- 0
  }
  for (j in (1:cols)) {
    sum <- sum + rej[,j]
    count <- count + 1;
  }
}

rej <- sum/count;

ps.options(horizontal=FALSE,onefile=FALSE)
ps.options(pagecentre=TRUE)
ps.options(paper="letter")
ps.options(width=7.0,height=10.0)

rows <- nrow(eta)

filename <- paste(prefix,".sci_sub_parm.eps",sep="")
postscript(file=filename)
par(mfrow=c(rows+1,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right)
for (j in 1:rows) {
  plot(idx,eta[j,],type="n",ylab="")
    lines(idx,eta[j,],lty="solid")
    top <- max(eta[j,])
    bot <- min(eta[j,])
    #for (k in 0:(num*(last-first+1))) {
    #  lines(c(lsub*k,lsub*k),c(bot,top),lty="dotted")
    #}
}
plot(idx,logl,type="n",ylab="")
  lines(idx,logl,lty="solid")
  top <- max(logl)
  bot <- min(logl)
  #for (k in 0:(num*(last-first+1))) {
  #  lines(c(lsub*k,lsub*k),c(bot,top),lty="dotted")
  #}

dev.off()


teta <- t(eta)
teta <- as.data.frame(teta)

filename <- paste(prefix,".sci_sub_pairs.eps",sep="")
postscript(file=filename)
par(mfrow=c(1,1),mar=c(4.5,4.0,1.5,1)+0.1,pty="s",xaxt="n",yaxt="n")
pairs(teta)
dev.off()

filename <- paste(prefix,".sci_sub_acf.eps",sep="")
postscript(file=filename)
par(mfrow=c(rows,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right)
for (j in 1:rows) {
  aut <- acf(teta[,j],lag.max=alag,type="correlation",plot=F,demean=T)
  y <- aut$acf
  aidx <- 0:(length(y)-1)
  plot(c(0,alag),c(0,1),type="n",ylab="")
    lines(aidx,y,lty="solid")
}
dev.off()

filename <- paste(prefix,".sci_sub_den.eps",sep="")
postscript(file=filename)
par(mfrow=c(rows,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right)
for (j in 1:rows) {
  den <- density(eta[j,],adjust=4)
  plot(den$x,den$y,type="n",ylab="",yaxt="n")
    lines(den$x,den$y,lty="solid")
}
dev.off()

corr <- cor(teta)
print(corr)
for (j in (1:rows)) {
  for (i in (1:rows)) {
    if (abs(corr[i,j]) < 0.6) {corr[i,j] <- 0.0}
  }
}
print(corr)

for (i in (1:length(rej))) print(rej[i])


