prefix <- "hab"
first <-   0
last <-    9
stride <- 10
alag <-  100

vars <- c(1,2,3,4,5,6,7)
names <- c("g","sig","sig_w","rho","phi","delta","gamma")

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

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;

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

rownames(rho) <- names[vars]
plotnames <- 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",".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=plotnames[j])
    lines(idx,rho[j,],lty="solid")
}
plot(idx,pi,type="n",ylab="",main="log likelihood + log prior")
  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
  aidx <- 0:alag
  plot(c(0,alag),c(0,1),type="n",ylab="",main=plotnames[j])
    lines(aidx,y,lty="solid")
  acf50[vars[j]] <- y[50]
}
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) {
   hist.scott(rho[j,],main=paste("Histogram of",plotnames[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)

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

for (j in 1:rows) {
  print(c(min(rho[j,]),max(rho[j,])),digits=5)
}

print("rho rejection rates")
for (i in (1:length(rej))) print(rej[i])

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

print("rho means to five digits")
for (j in 1:rows) {
  print(mean(rho[j,]),digits=5)
}

print("rho sdev to five digits")
for (i in (1:rows)) {
  print(sqrt(var[i,i]),digits=5)
}

