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

skip <-   1
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=10.0)
  
for (name in c(".sci_prior.",".stat_logl.",".stat_prior.")) {
  for (i in (first:last)) {
    if (i < 10) {
     filename <- paste(prefix,name,"00",i,".dat",sep="")
    } else {
     filename <- paste(prefix,name,"0",i,".dat",sep="")
    }
    tmp <- scan(filename)
    rows <- tmp[1]
    cols <- tmp[2]
    xi <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
    if (i == first) {
      x <- xi
    } else {
      x <- cbind(x,xi)
    }
  }
  if (name == ".sci_prior.") {
    pi <- x
  } else {
    pi = pi + x
  }
}

idx <- seq(1,ncol(pi),skip)
pi <- pi[,idx]

for (i in (first:last)) {
  if (i < 10) {
   filename <- paste(prefix,".sci_rej.","00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".sci_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
  }
  sum <- sum + rej[,1]
  count <- count + 1;
}

rej <- sum/count;

for (i in (first:last)) {
  if (i < 10) {
   filename <- paste(prefix,".sci_parm.","00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".sci_parm.","0",i,".dat",sep="")
  }
  tmp <- scan(filename)
  rows <- tmp[1]
  cols <- tmp[2]
  thetai <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
  if (i == first) {
    theta <- thetai
    print(dim(theta))
  } else {
    theta <- cbind(theta,thetai)
    print(dim(theta))
  }
}

idx <- seq(1,ncol(theta),skip)
theta <-theta[,idx]

for (i in (first:last)) {
  if (i < 10) {
   filename <- paste(prefix,".sci_func.","00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".sci_func.","0",i,".dat",sep="")
  }
  tmp <- scan(filename)
  rows <- tmp[1]
  cols <- tmp[2]
  statsi <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
  if (i == first) {
    stats <- statsi
    print(dim(stats))
  } else {
    stats <- cbind(stats,statsi)
    print(dim(stats))
  }
}

idx <- seq(1,ncol(stats),skip)
stats <-stats[,idx]

if (is.vector(theta)) {
  rows <- 1
} else {
  rows <- nrow(theta)
}

filename <- paste(prefix,".sci_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)
if (is.vector(theta)) {
  plot(idx,theta,type="n",ylab="")
      lines(idx,theta,lty="solid")
} else {
  for (j in 1:rows) {
    plot(idx,theta[j,],type="n",ylab="")
      lines(idx,theta[j,],lty="solid")
  }
}
plot(idx,pi,type="n",ylab="")
  lines(idx,pi,lty="solid")
dev.off()

if (!is.vector(theta)) {
  ttheta <- t(theta)
  ttheta <- as.data.frame(ttheta)
  
  filename <- paste(prefix,".sci_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(ttheta)
  dev.off()
  
  filename <- paste(prefix,".sci_.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(ttheta[,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_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(theta[j,],adjust=4)
    plot(den$x,den$y,type="n",ylab="",yaxt="n")
      lines(den$x,den$y,lty="solid")
  }
  dev.off()
  
  corr <- cor(ttheta)
  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])
}
