prefix <- "sv"
first <-  0
last <-   9
skip <-  10
alag <- 100

library(MASS)
library(ts)

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,".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),skip)
stats <-stats[,idx]

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),skip)
stats <-stats[,idx]

rows <- nrow(stats)

filename <- paste(prefix,".stats.chain.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) {
  plot(idx,stats[j,],type="n",ylab="")
    lines(idx,stats[j,],lty="solid")
}
dev.off()

tstats <- t(stats)
tstats <- as.data.frame(tstats)

filename <- paste(prefix,".stats.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(tstats)
dev.off()

filename <- paste(prefix,".stats.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(tstats[,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="")
    lines(aidx,y,lty="solid")
}
dev.off()

filename <- paste(prefix,".stats.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) {
  den <- density(stats[j,],adjust=4)
  plot(den$x,den$y,type="n",ylab="",yaxt="n")
    lines(den$x,den$y,lty="solid")
}
dev.off()

corr <- cor(tstats)
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)