prefix <- "sv"

first  <-  0
last   <-  5
stride <- 20

alag <-  100

beg   <- c(1, 6, 11)
end   <- c(5,10, 13)
pages <- 3

pnames <- c("a0","a1","a2","b0","b1","b2","c0","c1",
            "r21","r22","r31","r32","r33")

sbeg   <- c(1, 6, 11)
send   <- c(5,10, 13)
spages <- 3

snames <- c("y1min","y2min","y1max","y2max","y1mean","y2mean",
            "y1sdev","y2sdev","y1y2cor","vmin","vmax","vmean","vsdev")

plots.per.page <- 0;
for (p in 1:pages) {
  if (p<pages) {
    plots.per.page <- max(plots.per.page,end[p]-beg[p]+1)
  } else {
    plots.per.page <- max(plots.per.page,end[p]-beg[p]+2)
  }
}

splots.per.page <- 0;
for (p in 1:spages) {
  splots.per.page <- max(splots.per.page,send[p]-sbeg[p]+1)
}

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 (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))

rownames(rho) <-pnames

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))

rownames(stats) <- snames

rows <- nrow(rho)

for (p in 1:pages) {
  if (p < 10) {
    filename <- paste(prefix,".rho.chain.0",p,".eps",sep="")
  } else {
    filename <- paste(prefix,".rho.chain.",p,".eps",sep="")
  }
  postscript(file=filename)
  par(mfrow=c(plots.per.page,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(b,l,t,r)
  for (j in beg[p]:end[p]) {
    plot(idx,rho[j,],type="n",ylab="",main=pnames[j])
     lines(idx,rho[j,],lty="solid")
  }
  if (p<pages) dev.off()
}
plot(idx,pi,type="n",ylab="",main="posterior")
  lines(idx,pi,lty="solid")
dev.off()

for (p in 1:pages) {
  trho <- t(rho[beg[p]:end[p],])
  trho <- as.data.frame(trho)
  if (nrow(rho) > 1) {
    if (p < 10) {
      filename <- paste(prefix,".rho.pairs.0",p,".eps",sep="")
    } else {
      filename <- paste(prefix,".rho.pairs.",p,".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()
  }
}

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

for (p in 1:pages) {
  if (p < 10) {
    filename <- paste(prefix,".rho.acf.0",p,".eps",sep="")
  } else {
    filename <- paste(prefix,".rho.acf.",p,".eps",sep="")
  }
  postscript(file=filename)
  par(mfrow=c(plots.per.page,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(b,l,t,r)    
  for (j in beg[p]:end[p]) {
    aut <- acf(trho[,j],lag.max=alag,type="correlation",plot=F,demean=T)
    y <- aut$acf
    aidx <- 0:alag
    plot(c(0,alag),c(0,1),type="n",ylab="",main=pnames[j])
      lines(aidx,y,lty="solid")
  }
  dev.off()
}
  
for (p in 1:pages) {
  if (p < 10) {
    filename <- paste(prefix,".rho.density.0",p,".eps",sep="")
  } else {
    filename <- paste(prefix,".rho.density.",p,".eps",sep="")
  }
  postscript(file=filename)
  par(mfrow=c(plots.per.page,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(b,l,t,r)    
  for (j in beg[p]:end[p]) {
    den <- density(rho[j,],adjust=4)
    plot(den$x,den$y,type="n",ylab="",yaxt="n",main=pnames[j])
      lines(den$x,den$y,lty="solid")
  }
  dev.off()
}

for (p in 1:spages) {
  if (p < 10) {
    filename <- paste(prefix,".stats.chain.0",p,".eps",sep="")
  } else {
    filename <- paste(prefix,".stats.chain.",p,".eps",sep="")
  }
  postscript(file=filename)
  par(mfrow=c(splots.per.page,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(b,l,t,r)
  for (j in sbeg[p]:send[p]) {
    plot(idx,stats[j,],type="n",ylab="",main=snames[j])
      lines(idx,stats[j,],lty="solid")
  }
  dev.off()
}

#if (nrow(rho) > 1) {
#  corr <- cor(trho)
#  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)
#}