# 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]
objective1 [1:5]

bits.saved2[1:5]
objective2 [1:5]

# #############################################################
#
# 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()


quit()