prefix <- "an"
suffix <- "3f"
first <-   0
last <-    9
stride <- 10
alag <-  100
vars <- c(1,2,3,4,5,6)
names <- c("mu_c","rho_c","sigma_c","kappa_c","mu_pi","sigma_pi",
           "gamma_pi","beta","p_a")

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

acfln <- rows
acf50 <- mat.or.vec(acfln,1)

rho <- rho[vars,]
rownames(rho) <- names[vars]
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.",suffix,".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="",main=names[j])
    lines(idx,rho[j,],lty="solid")
}
plot(idx,pi,type="n",ylab="",main="likelihood")
  lines(idx,pi,lty="solid")
dev.off()

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

if (nrow(rho) > 1) {
  filename <- paste(prefix,".rho.pairs.",suffix,".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.",suffix,".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
  aidx <- 0:alag
  plot(c(0,alag),c(0,1),type="n",ylab="",main=names[j])
    lines(aidx,y,lty="solid")
  acf50[vars[j]] <- y[50]
}
dev.off()

filename <- paste(prefix,".rho.density.",suffix,".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) {
   hist.scott(rho[j,],main=paste("Histogram of",names[j]))
#  den <- density(rho[j,],adjust=1)
#  plot(den$x,den$y,type="n",ylab="",yaxt="n")
#    lines(den$x,den$y,lty="solid")
}
dev.off()

var <- var(trho)
print(var)
for (i in (1:rows)) {
  print(sqrt(var[i,i]))
}

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

print(paste("stride",stride))
print("acf50")
for (i in 1:acfln) print(acf50[i],digits=2)