prefix <- "sv"
first <- 0
last <-  9

skip <-  10
alag <- 100

library(MASS)

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

for (inode in (1:15)) {

  if (inode < 10) {
    node <- paste("00",inode,sep="")
  } else {
    node <- paste("0",inode,sep="")
  }
    
  for (name in c(".sci_prior.",".stat_logl.",".stat_prior.")) {
    for (i in (first:last)) {
      if (i < 10) {
       filename <- paste(prefix,name,node,".00",i,sep="")
      } else {
       filename <- paste(prefix,"name",node,".0",i,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.",node,".00",i,sep="")
    } else {
     filename <- paste(prefix,".sci_rej.",node,".0",i,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.",node,".00",i,sep="")
    } else {
     filename <- paste(prefix,".sci_parm.",node,".0",i,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.",node,".00",i,sep="")
    } else {
     filename <- paste(prefix,".sci_func.",node,".0",i,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]
  
  rows <- nrow(theta)
  
  filename <- paste(prefix,".sci_parm.",node,".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,theta[j,],type="n",ylab="")
      lines(idx,theta[j,],lty="solid")
  }
  plot(idx,pi,type="n",ylab="")
    lines(idx,pi,lty="solid")
  dev.off()
  
  ttheta <- t(theta)
  ttheta <- as.data.frame(ttheta)
  
  filename <- paste(prefix,".sci_pairs.",node,".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.",node,".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.",node,".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])

}
