> # R script to create a dual-histogram ps file from a predict-2nd bychain files.
> 
> # Usage: R --vanilla --quiet --args 'fname1' 'fname2' < bychain2dualhist.r >& R.bychain2dualhist.log
> 
> # input  file:  fname1.by-chain
> # input  file:  fname2.by-chain
> # output file:  fname1.fname2.by-chain.dualhist.pdf
> 
> infile.bychain.suffix      <- ".by-chain"
> outfile.hist.suffix    <- ".by-chain.dualhist"
> 
> # find the fname arguments after --args
> my.args <- commandArgs()
> my.argnum <- which(my.args == "--args")
> fname1 <- NULL
> fname1 <- my.args[my.argnum+1]
> fname2 <- NULL
> fname2 <- my.args[my.argnum+2]
> 
> if (!length(fname1)) {
+     print("Error: expecting --args 'fname1' 'fname2'") 
+     quit() 
+ }
> 
> if (!length(fname2)) {
+     print("Error: expecting --args 'fname1' 'fname2'") 
+     quit() 
+ }
> 
> infile1.prefix <- fname1
> infile2.prefix <- fname2
> outfile.prefix <- paste (fname1, fname2, sep=".")
> 
> infile1.bychain   <- paste (infile1.prefix, infile.bychain.suffix,  sep="")
> infile2.bychain   <- paste (infile2.prefix, infile.bychain.suffix,  sep="")
> outfile.hist.ps   <- paste (outfile.prefix, outfile.hist.suffix,".ps", sep="")
> outfile.hist.pdf  <- paste (outfile.prefix, outfile.hist.suffix,".pdf",sep="")
> 
> 
> bychain1.data <- read.table( infile1.bychain, header=FALSE, as.is=TRUE)
> bychain2.data <- read.table( infile2.bychain, header=FALSE, as.is=TRUE)
> 
> # expected format of by-chain file:
> ## chainID             Bits_saved IDBits  Q7      SOV  object   count   SOV_E   SOV_B   SOV_G  ...
> #1ollA                 1.1881 -------  0.6543  0.6722  2.1785     188  0.8434  1.0000  0.1667  ...
> #1jnrA                 0.8839 -------  0.6199  0.5452  1.7765     642  0.6928  0.0000  0.1071  ...
> 
> bits.saved1   <- bychain1.data$V2
> objective1    <- bychain1.data$V6
> 
> bits.saved2   <- bychain2.data$V2
> objective2    <- bychain2.data$V6
> 
> # sanity check
> bits.saved1[1:5]
[1] 0.5790 0.8464 0.6850 0.8858 1.7558
> objective1 [1:5]
[1] 1.3245 1.6771 1.4314 1.8086 3.0051
> 
> bits.saved2[1:5]
[1] 0.1741 0.5392 0.7445 0.7687 1.6769
> objective2 [1:5]
[1] 0.8493 1.3152 1.5400 1.5883 2.8134
> 
> # #############################################################
> #
> # PLOTTED OUTPUT BEGINS HERE
> #
> 
> # #############
> # histogram plots
> 
> # function for plotting histograms with labels containing mean,sd
> my.plothist <- function (data1, data2, data.nm)
+ {   # arguments
+     # data1:    data to be plotted
+     # data2:    data to be plotted
+     # data.nm   character string of name of data
+ 
+     # get integer min and max of data to compute breaks
+     # as.integer() truncates towards zero, need to round down for negative min
+     if (min(data1,data2,na.rm=TRUE) <  0) data.min <- as.integer(min(data1,data2,na.rm=TRUE)) -1
+     if (min(data1,data2,na.rm=TRUE) >= 0) data.min <- as.integer(min(data1,data2,na.rm=TRUE))
+     data.max <- as.integer(max(data1,data2,na.rm=TRUE)) + 1
+ 
+     data.brk <- seq(data.min,data.max,0.1)
+ 
+     data1.mean <- mean(data1,na.rm=TRUE)
+     data1.sd   <- sd  (data1,na.rm=TRUE)
+     data1.sem  <- data1.sd/sqrt((length(data1)-1))
+ 
+     data2.mean <- mean(data2,na.rm=TRUE)
+     data2.sd   <- sd  (data2,na.rm=TRUE)
+     data2.sem  <- data2.sd/sqrt((length(data2)-1))
+ 
+     data.main <- paste ("Histogram of ",data.nm, sep="")
+ 
+     data1.xlab <- paste(fname1,
+                         " mean=",
+                         formatC(data1.mean,format="f",digits=2),
+                         " std.dev=",
+                         formatC(data1.sd,  format="f",digits=2),
+                         " std.err.mean=",
+                         formatC(data1.sem, format="f",digits=3),
+                         sep="")
+ 
+     data2.xlab <- paste(fname2,
+                         " mean=",
+                         formatC(data2.mean,format="f",digits=2),
+                         " std.dev=",
+                         formatC(data2.sd,  format="f",digits=2),
+                         " std.err.mean=",
+                         formatC(data2.sem, format="f",digits=3),
+                         sep="")
+ 
+     data.xlab <- paste(data1.xlab, data2.xlab, sep="\n")
+ 
+     data1.hist <- hist(data1, breaks = data.brk, freq=FALSE, plot=FALSE)
+     data2.hist <- hist(data2, breaks = data.brk, freq=FALSE, plot=FALSE)
+           
+     ymax  <- max(data1.hist$density, data2.hist$density)
+ 
+     data1.lty <- 2
+     data1.col <- "red"
+ 
+     data2.lty <- 3
+     data2.col <- "blue"
+ 
+     plot  ( data1.hist$mids,  data1.hist$density,
+             main  = data.main,
+ 	    xlab  = data.xlab,
+ 	    ylab  = "density",
+             ylim  = c(0, ymax),
+             type    = "l", lty=data1.lty, col=data1.col,
+           )
+ 
+ 
+    points  ( data2.hist$mids,  data2.hist$density,
+              type    = "l", lty=data2.lty, col=data2.col,
+            )
+ 
+    legend  ( x=1.8, y=1.0,
+              legend = c( fname1,          fname2        ),
+              lty    = c( data1.lty,       data2.lty     ),
+              col    = c( data1.col,       data2.col     ),
+            )
+ 
+ }  # end of function
>                                 
> # set output device and filename
> #pdf(outfile.hist.pdf)
> postscript(outfile.hist.ps,
+     horizontal=FALSE, onefile=FALSE,
+     paper="special", width=6, height=4)
> 
> my.plothist (bits.saved1, bits.saved2,"bits.saved")
> #my.plothist (objective1, objective2,  "objective")
> 
> 
> # restore output device to default
> dev.off()
null device 
          1 
> 
> 
> quit()