prefix <- "sv97"
first <-   0
last <-    9
stride <- 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=10.0)
  
  
for (i in (first:last)) {
  if (i < 10) {
   filename <- paste(prefix,".reject.00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".reject.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,".pi.00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".pi.0",i,".dat",sep="")
  }
  tmp <- scan(filename)
  rows <- tmp[1]
  cols <- tmp[2]
  pii <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
  if (i == first) {
    pi <- pii
  } else {
    pi <- cbind(pi,pii)
  }
}

idx <- seq(1,ncol(pi),stride)
pi <- pi[1,idx]
if (is.vector(pi)) pi <- matrix(pi,nrow=1,ncol=length(pi))

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

idx <- seq(1,ncol(rho),stride)
rho <-rho[,idx]
if (is.vector(rho)) rho <- matrix(rho,nrow=1,ncol=length(rho))

for (i in (first:last)) {
  if (i < 10) {
   filename <- paste(prefix,".stats.00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".stats.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),stride)
stats <-stats[,idx]
if (is.vector(stats)) stats <- matrix(stats,nrow=1,ncol=length(stats))

rows <- nrow(rho)

filename <- paste(prefix,".rho.chain.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,rho[j,],type="n",ylab="")
    lines(idx,rho[j,],lty="solid")
}
plot(idx,pi,type="n",ylab="")
  lines(idx,pi,lty="solid")
dev.off()

trho <- t(rho)
trho <- as.data.frame(trho)

if (nrow(rho) > 1) {
  filename <- paste(prefix,".rho.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(trho)
  dev.off()
}

filename <- paste(prefix,".rho.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(trho[,j],lag.max=alag,type="correlation",plot=F,demean=T)
  y <- aut$acf
  lo <- min(0,min(y))
  aidx <- 0:alag
  plot(c(0,alag),c(lo,1),type="n",ylab="")
    lines(aidx,y,lty="solid")
}
dev.off()

filename <- paste(prefix,".rho.density.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(rho[j,],adjust=4)
  plot(den$x,den$y,type="n",ylab="",yaxt="n")
    lines(den$x,den$y,lty="solid")
}
dev.off()

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