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

pincl  <- c(T,T,T,T,T, T,T,T,T,T, T,T,T,T,T, T,T,T,T,T, T,T)
pbeg   <- c(1, 5,  9, 11, 15, 18, 20)
pend   <- c(4, 8, 10, 14, 17, 19, 22)
ppages <- 7

pnames <-
  c("a0[1]","a0[2]","a0[3]","a0[4]","a0[5]","a0[6]","a0[7]","a0[8]",
  "b0[1]","b0[2]","B(1,1)","B(2,1)","B(1,2)","B(2,2)",
  "R0[1]","R0[2]","R0[3]", "P(1,1)","P(2,1)",
  "cmf_parms[1]","cmf_parms[2]","cmf_parms[3]")

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

fnames <- 
  c("m[1]","m[2]","m[3]","m[4]")

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_chain.00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".pi_chain.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,".parms_chain.","00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".parms_chain.","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)
acf2  <- acf50

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_chain.","00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".stats_chain.","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,".parms.chain.0",p,".eps",sep="")
  } else {
    filename <- paste(prefix,".parms.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,".parms.pairs.0",p,".eps",sep="")
#    } else {
#      filename <- paste(prefix,".parms.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,".parms.acf.0",p,".eps",sep="")
  } else {
    filename <- paste(prefix,".parms.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]
      acf2lags = min(alag,2)
      acf2[j] <- y[acf2lags]
    }
  }
  dev.off()
}

for (p in 1:ppages) {
  if (p < 10) {
    filename <- paste(prefix,".parms.density.0",p,".eps",sep="")
  } else {
    filename <- paste(prefix,".parms.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 autocorrelations at lag 2")
for (i in 1:acfln) print(paste(pnames[i],"  ",format(acf2[i],digits=8)))
#for (i in 1:acfln) {
#  print(paste(pnames[i],"  ",format(acf2[i],digits=8,scientific=TRUE)))
#}

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

probs <- c(0,0.01,0.05,0.25,0.5,0.75,0.95,0.99,1)

print("stats quantiles")
for (i in (1:nrow(func))) {
  #print(paste(fnames[i],"  ",quantile(func[i,],probs),digits=5))
  print(quantile(func[i,],probs))
}
