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

pincl  <- c(T,T,T,T,T,T,T)
pbeg   <- c(1)
pend   <- c(7)
ppages <- 1

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

fincl  <- c(T,T,T,T,T,T,T)
fbeg   <- c(1)
fend   <- c(5)
fpages <- 1

fnames <- c("min","q_25","med","q_75","max")

pplots.per.page <- 0;
for (p in 1:ppages) {
  if (p<ppages) {
    pplots.per.page <- max(pplots.per.page,pend[p]-pbeg[p]+1)
  } else {
    pplots.per.page <- max(pplots.per.page,pend[p]-pbeg[p]+2)
  }
}

fplots.per.page <- 0;
for (p in 1:fpages) {
  fplots.per.page <- max(fplots.per.page,fend[p]-fbeg[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,".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,".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)

rownames(rho) <- pnames
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]
  funci <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
  if (i == first) {
    func <- funci
    print(dim(func))
  } else {
    func <- cbind(func,funci)
    print(dim(func))
  }
}

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

rows <- nrow(rho)

for (p in 1:ppages) {
  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(pplots.per.page,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(b,l,t,r)
  for (j in pbeg[p]:pend[p]) {
    if (pincl[j]) {
      plot(idx,rho[j,],type="n",ylab="",main=pnames[j])
        lines(idx,rho[j,],lty="solid")
    }
  }
  if (p<ppages) dev.off()
}
plot(idx,pi,type="n",ylab="",main="log likelihood + log prior")
  lines(idx,pi,lty="solid")
dev.off()

for (p in 1:fpages) {
  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(fplots.per.page,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(b,l,t,r)
  for (j in fbeg[p]:fend[p]) {
    if (fincl[j]) {
      plot(idx,func[j,],type="n",ylab="",main=fnames[j])
        lines(idx,func[j,],lty="solid")
    }
  }
  dev.off()
}

#for (p in 1:ppages) {
#  trho <- t(rho[pbeg[p]:pend[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:ppages) {
  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(pplots.per.page,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(b,l,t,r)
  for (j in pbeg[p]:pend[p]) {
    if (pincl[j]) {
      aut <- acf(trho[,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="",main=pnames[j])
        lines(aidx,y,lty="solid")
      acf50lags = min(alag,50)
      acf50[j] <- y[acf50lags]
    }
  }
  dev.off()
}

for (p in 1:ppages) {
  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(pplots.per.page,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(b,l,t,r)
  for (j in pbeg[p]:pend[p]) {
    if (pincl[j]) {
      hist.scott(rho[j,],main=pnames[j])
      #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()
}

#tfunc <- t(func)
#tfunc <- as.data.frame(tfunc)
#
#for (p in 1:fpages) {
#  if (p < 10) {
#    filename <- paste(prefix,".stats_density.0",p,".eps",sep="")
#  } else {
#    filename <- paste(prefix,".stats_density.",p,".eps",sep="")
#  }
#  postscript(file=filename)
#  par(mfrow=c(pplots.per.page,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(b,l,t,r)
#  for (j in fbeg[p]:fend[p]) {
#    if (fincl[j]) {
#      hist.scott(func[j,],main=fnames[j])
#      #den <- density(func[j,],adjust=4)
#      #plot(den$x,den$y,type="n",ylab="",yaxt="n")
#      #  lines(den$x,den$y,lty="solid")
#    }
#  }
#  dev.off()
#}

print("pi range")
for (j in 1:nrow(pi)) {
  print(c(min(pi[j,]),max(pi[j,])),digits=5)
}

trho <- t(rho[pincl,])
trho <- as.data.frame(trho)
tnames <- pnames[pincl];

if (nrow(rho) > 1) {
  print("rho correlations")
  corr <- cor(trho)
  print(corr)
  for (j in (1:ncol(trho))) {
    for (i in (1:ncol(trho))) {
      if (abs(corr[i,j]) < 0.6) {corr[i,j] <- 0.0}
    }
  }
  print("rho correlations larger than 0.6")
  print(corr)
} else {
  print("variance of rho")
  var <- var(trho)
  print(var)
}

print("rho means and standard deviations")
for (i in (1:nrow(rho))) {
  line <- paste(pnames[i],"  ",mean(rho[i,]),"  ",sd(rho[i,]))
  print(line)
}

print("stats means and standard deviations")
for (i in (1:nrow(func))) {
  line <- paste(fnames[i],"  ",mean(func[i,]),"  ",sd(func[i,]))
  print(line)
}

rnames <- c(pnames,"total")

print("rho rejection rates")
for (i in (1:length(rej))) print(paste(rnames[i],"  ",format(rej[i],digits=5)))

print(paste("rho autocorrelations at lag",acf50lags))
for (i in 1:acfln) print(paste(pnames[i],"  ",format(acf50[i],digits=2)))

print("rho means to five digits, more digits printed above")
for (i in (1:nrow(rho))) {
  print(paste(pnames[i],"  ",format(mean(rho[i,]),digits=5)))
}

print("rho sdevs to five digits, more digits printed above")
for (i in (1:nrow(rho))) {
  print(paste(pnames[i],"  ",format(sd(rho[i,]),digits=5)))
}

print("rho t-stats to five digits")
for (i in (1:nrow(rho))) {
  print(paste(pnames[i],"  ",format(mean(rho[i,])/sd(rho[i,]),digits=5)))
}

print("stats means to five digits, more digits printed above")
for (i in (1:nrow(func))) {
  print(paste(fnames[i],"  ",format(mean(func[i,]),digits=5)))
}

print("stats sdevs to five digits, more digits printed above")
for (i in (1:nrow(func))) {
  print(paste(fnames[i],"  ",format(sd(func[i,]),digits=5)))
}

print("stats t-stats to five digits")
for (i in (1:nrow(func))) {
  print(paste(fnames[i],"  ",format(mean(func[i,])/sd(func[i,]),digits=5)))
}
