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

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,".sci_parm.","00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".sci_parm.","0",i,".dat",sep="")
  }
  tmp <- scan(filename)
  rows <- tmp[1]
  cols <- tmp[2]
  thetai <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
  if (i == first) {
    theta <- thetai
    print(dim(theta))
  } else {
    theta <- cbind(theta,thetai)
    print(dim(theta))
  }
}

rows <- nrow(theta)
cols <- ncol(theta)

idx0 <- (1:(cols-1))
idx1 <- (2:cols)

tmp <- theta[,idx0]-theta[,idx1]

for (i in (1:rows)) {
  thetai <- tmp[i,]
  print(hist(thetai,plot=FALSE)$counts)
}

for (i in (first:last)) {
  if (i < 10) {
   filename <- paste(prefix,".sci_rej.","00",i,".dat",sep="")
  } else {
   filename <- paste(prefix,".sci_rej.","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,cols);
    count <- 0
  }
  sum <- sum + rej
  count <- count + 1;
}

print(sum)

ave <- sum/count;

print(ave)

