prefix <- "sv"
first <- 0
last <-  9

nodes <- c("001","002","003","004","005","006","007","008","009","010",
           "011","012","013","014","015")

skip <-  25
alag <- 100

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 (name in c(".sci_prior.",".stat_logl.",".stat_prior.")) {
  for (node in nodes) {
    for (i in (first:last)) {
      if (i < 10) {
       filename <- paste(prefix,name,node,".00",i,sep="")
      } else {
       filename <- paste(prefix,"name",node,".0",i,sep="")
      }
      tmp <- scan(filename)
      rows <- tmp[1]
      cols <- tmp[2]
      xi <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
      if (node == nodes[1]) {
        x <- xi
      } else {
        x <- cbind(x,xi)
      }
    }
  }
  if (name == ".sci_prior.") {
    pi <- x
  } else {
    pi = pi + x
  }
}

idx <- seq(1,ncol(pi),skip)
pi <- pi[,idx]

for (node in nodes) {
  for (i in (first:last)) {
    if (i < 10) {
     filename <- paste(prefix,".sci_rej.",node,".00",i,sep="")
    } else {
     filename <- paste(prefix,".sci_rej.",node,".0",i,sep="")
    }
    tmp <- scan(filename)
    rows <- tmp[1]
    cols <- tmp[2]
    rej <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
    if (node == nodes[1]) {
      sum <- mat.or.vec(rows,1);
      count <- 0
    }
    sum <- sum + rej[,1]
    count <- count + 1;
  }
}

rej <- sum/count;

for (node in nodes) {
  for (i in (first:last)) {
    if (i < 10) {
     filename <- paste(prefix,".sci_parm.",node,".00",i,sep="")
    } else {
     filename <- paste(prefix,".sci_parm.",node,".0",i,sep="")
    }
    tmp <- scan(filename)
    rows <- tmp[1]
    cols <- tmp[2]
    thetai <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
    if (node == nodes[1]) {
      theta <- thetai
      print(dim(theta))
    } else {
      theta <- cbind(theta,thetai)
      print(dim(theta))
    }
  }
}

idx <- seq(1,ncol(theta),skip)
theta <-theta[,idx]

for (node in nodes) {
  for (i in (first:last)) {
    if (i < 10) {
     filename <- paste(prefix,".sci_func.",node,".00",i,sep="")
    } else {
     filename <- paste(prefix,".sci_func.",node,".0",i,sep="")
    }
    tmp <- scan(filename)
    rows <- tmp[1]
    cols <- tmp[2]
    statsi <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
    if (node == nodes[1]) {
      stats <- statsi
      print(dim(stats))
    } else {
      stats <- cbind(stats,statsi)
      print(dim(stats))
    }
  }
}

idx <- seq(1,ncol(stats),skip)
stats <-stats[,idx]

rows <- nrow(theta)

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

ttheta <- t(theta)
ttheta <- as.data.frame(ttheta)

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

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

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

corr <- cor(ttheta)
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 (i in (1:length(rej))) print(rej[i])

