###shivplot
###Owen Woody and Robert Nadon, 2006
###Direct questions/concerns about this file to owoody@uwaterloo.ca

####Table of Contents:  (search for the @tag to skip directly to a section of the document)
#@shivargs:   argument descriptions for shivplot function.
#@shivplot:   shivplot code, style of figures 4-5 from manuscript.
#@multiR:     multivariate shivplot, R-language support.
#@multiS:     multivariate shivplot, S-language support.
#@textFigs:   selected code to produce figures from the text.
#@examples:   code producing the example figures
#@support:    the help text, code to detect and accomodate differences between R and S-PLUS

#@shivargs
#shivplot argument list:
#data,                  #either a matrix with rows representing subjects and columns representing observations, or a list composed of a number of such matrices
#locationFCN=mean,          #The function to use to obtain an estimate of location across observations.  The default measure (and the one used throughout this document) is the mean.
#plotNames,                 #Names to attribute to each individual shivplot.  The names will be located above the plots, in the plot frame.  The number of plotNames must match the length of the list used in the 'data' argument.  If these subnames are not desired, enter the string "" for each name.
#plotColors=c(1),           #What colors to use for each created plot.  Expects a vector with a number of entries equal to the length of the list in the 'data' argument.  S-PLUS users must specify a vector of integers for colors, whereas R users may either use numbers, calls to the function "colors()", or strings representing color names.
#outlierColor=2,            #The color to use to mark observations that qualify as "outliers" under default boxplot specifications.
#scaleTop=T,                #When multiple shivplots are drawn in tandem (denoted through this manuscript as "comparative shivplots"), the user has the option to enforce a common scale for the "top" half of the y-axis.  This argument expects a Boolean variable toggling this feature.
#scaleBot=T,                #As above, with bottom half of y-axis
#FUN1=sdVsLocPlugin,            #The function to use to compute the line on the top half of the shivplot.  By default, this function relates standard deviations across replicates to the selected measure of centre.  These functions must follow a special (but simple) 'plugin' format, see the provided code for details.
#FUN2=densityPlugin,            #The 'plugin' function to use to compute the line in the lower half of the shivplot.
#xLab="Mean across replicates",    #The label for the shared x-axis for a single or comparative boxplot.  Similar to xlab from standard R/S-PLUS plotting.  Be sure to change this to reflect the choice of location function (locationFCN).
#boxWidth=.3,               ##
#whiskWidth=.5,             ##Controls the length of the ticks drawn to denote the "box" region from a boxplot, and the "whisker" region denoting the limit past which observations are deemed to be outliers.  Negative values cause the ticks to be reflected about the x-axis.
#outlierCharacter=-1,          #The plotting character to use to mark outliers.  The default (-1) is the system default.

#@shivplot
###This function produces shivplots. 
shivplot<-function(
data,
locationFCN=mean,
plotNames=NULL,
plotColors=1,
outlierColor=2,
scaleTop=T,
scaleBot=T,
FUN1=sdVsLocPlugin,
FUN2=densityPlugin,
xLab="Mean intensity across replicate arrays",
boxWidth=.3,
whiskWidth=.5,
outlierCharacter=-1,
...)

{
    #coerce data frames into matrices to avoid unexpected behaviour.
    if (is.data.frame(data)) {
        print("encountered a data frame, attempting to convert it to matrix");
        data<-as.matrix(data)
    }
    #Determine the number of plots.  Assume there is 1, and if "data" is a list alter this assumption according to the number of entries in the list.
    numPlots=1
    if(is.list(data))
    {
        numPlots=length(data)
    }
    else
    {
        temp<-list()
        temp[[1]]<-data
        data<-temp
    }

    #If there are more plots than plot names, improvise additional names.
    if(length(plotNames) < numPlots) {
        n = numPlots - length(plotNames)
        for(i in 1:n)
        {
            plotNames = c(plotNames , paste("Unnamed", i, sep =" "))
        }
    }
    #If there are more plots than suggested colors, improvise additional colors.
    #This color selection is not smart, and may reuse colors provided by user.
    if(length(plotColors) < numPlots) {
        n = numPlots - length(plotColors)
        plotColors = c(plotColors , 2:(n+1))
    }

    plotData<-list()
    for(i in 1:numPlots)
    {
        #Compute all the required information to produce the graphs.  This is done prior to graphing (as opposed to during) so that graphs can be scaled relative to one another if desired.
        plotData[[i]]<-shiv.prepGraphData(data[[i]], FUN1, FUN2, locationFCN)
    }

    ##Using this plot data, determine the largest values for use in scaling.
    globals<-shiv.determineGlobals(numPlots, plotData)

    #Defining the left edge, right edge, and a unit of measurement relative to the plot width.
    shift=0.05*(globals$maxLoc-globals$minLoc)
    leftBorder= globals$minLoc - shift
    rightBorder= globals$maxLoc + shift

    ###For readability, store all graphing parameter data in a list.
    graphParams<-list(
    leftBorder=leftBorder,
    rightBorder=rightBorder,
    shift=shift,
    scaleTop=scaleTop,
    scaleBot=scaleBot,
    outlierColor=outlierColor,
    outlierCharacter=outlierCharacter,
    boxWidth=boxWidth,
    whiskWidth=whiskWidth
    )

    #Construct the plot device, for subsequent addition of the shivplots.
    shiv.plotSurface(numPlots, plotData, globals, plotNames, plotColors, xLab, graphParams,...)


}

##Calculates all the required plot information for a given data set.  Invoked once for each entry in the 'data' list.  Maximum and minimum values are recorded for comparison later in shiv.determineGlobals (for scaling)
shiv.prepGraphData<-function(datax, FUN1, FUN2, LOCATION.FCN)
{

    #epsilon = 0.000001  ###Used to avoid plotting an outlier at the upper end of the whisker.
    epsilon = 0
    ###general information, information regarding the boxplot ticks, information for the upper half of the plot, and the lower half are all stored in separate lists.
    myList<-list()
    boxpData<-list()
    topData<-list()
    botData<-list()

    #computing location across replicates, recording max, min, average location
    data.locs <- apply(datax, 1, LOCATION.FCN)
    myList$maxLoc <- max(data.locs)
    myList$minLoc <- min(data.locs)
    myList$midLoc <- median(data.locs)

    ##Calculate quantiles in style of S-PLUS/R boxplot
    boxpData$five <- sort(boxplot(data.locs, plot=F)$stats)###NOTE:  S-PLUS and R return $stats in opposite orders.  Applied 'sort' function to be consistent.

    ###for reading convenience.  Note that the extremes are selected to the be the largest values that do not exceed a distance of outThresh*iqr from the quartiles.
    lowerExtreme = boxpData$five[1]
    upperExtreme = boxpData$five[5]

    ###The following line of code attempts to emulate the outlier detection rule used in the R/S-PLUS boxplot function.
    ###outliers are stored, and plotted later.
    boxpData$outs <- c(data.locs[data.locs > (upperExtreme+epsilon)], data.locs[data.locs < (lowerExtreme-epsilon)])  


    topData<-FUN1(datax, LOCATION.FCN)  #Compute the line fit for the top half of the shivplot.
    botData<-FUN2(datax, LOCATION.FCN)#Compute the line fit for the bottom half of the shivplot.
    myList$top<-topData
    myList$bot<-botData
    myList$boxpData<-boxpData
    return(myList)
}


#Determines maximum and minimum values, used for scaling measurements.
shiv.determineGlobals<-function(numPlots, plotData)
{

    maxes<-NULL
    mins<-NULL
    middles<-NULL
    maxTops<-NULL
    maxBots<-NULL
    globals<-list()
    ##Finding global max/min for each of the location, top half, and bottom half.
    for(i in 1:numPlots)
    {
        maxes<-c(maxes, plotData[[i]]$maxLoc)
        mins<-c(mins, plotData[[i]]$minLoc)
        middles<-c(middles, plotData[[i]]$midLoc)
        maxTops<-c(maxTops, max(plotData[[i]]$top$y))
        maxBots<-c(maxBots, max(plotData[[i]]$bot$y))
    }
    globals$maxLoc <- max(maxes)
    globals$minLoc <- min(mins)
    globals$centreLoc <- median(middles)#Used to position tile labels (plotNames)
    globals$maxTop <- max(maxTops)
    globals$maxBot <- max(maxBots)
    return(globals)
}

##Prepares a plot surface, invokes shiv.plotIt to add shivplots to surface
shiv.plotSurface<-function(numPlots, plotData, globals, plotNames, plotColors, xLab, graphParams,...)
{
    ###Constructing the plot surface

    ###Constants for beauty/structure, may be changed as desired.
    BOUNDARY.CORRECTION=0.04###removes the upper and lower margin space automatically added by R/S
    NUMBER.OF.LINES.PER.PLOT=4
    TITLE.OFFSET=0.25###how far below dividing line to put tile name (plotName)

    ##Constructing the plot surface.  ylimits set to fit the required number of shivplots, xlimits set to reflect largest and smallest values of location measure across replicates.
    plot(0, 0,
    ylim = c(0+BOUNDARY.CORRECTION*(NUMBER.OF.LINES.PER.PLOT*numPlots), NUMBER.OF.LINES.PER.PLOT*numPlots-BOUNDARY.CORRECTION*(NUMBER.OF.LINES.PER.PLOT*numPlots)),
    xlim = c(globals$minLoc, globals$maxLoc),
    type = "n",

    xlab=xLab,
    ylab="", yaxt="n", xaxs = "r", ...)


    for(i in 1:numPlots)
    {
        j = (numPlots:1)[i]
        ##Construct each plot individually on the prepared surface.
        shiv.plotIt(plotData[[j]], globals, floor=0+(i-1)*4, plotColors[j],graphParams)
        if(i!=numPlots)
        abline(h=i*NUMBER.OF.LINES.PER.PLOT, lwd=3)###Add dividing line
        text(globals$centre, i*NUMBER.OF.LINES.PER.PLOT-TITLE.OFFSET, plotNames[j])###Add plot titles.
    }
}


###Function that adds shivplots to prepared plot surface.
shiv.plotItR<-function(data, globals, floor, plotColor, graphParams)

{
    #Constants
    ##Formats plot surface to have no y-axis, x-axis as normal
    DASHED.LINE = 2  ##Value for a dashed line
    SOLID.LINE = 1 ##Value for a solid line

    ###These values specify the margins within which each shivplot is contained.  The y-axis is divided into blocks, and each shiv plot is constructed within a 4(default) unit block.
    figFloor=floor
    centre = figFloor+2
    TOPLINE.COLOR = plotColor  ###by default, both the upper and lower halves of the shiv plot are drawn using the same color.  This can be changed here if desired.
    BOTLINE.COLOR = plotColor


    ###Recovering the five values required for boxplots
    lowEx = data$boxpData$five[1]
    lowHinge = data$boxpData$five[2]
    dat.median = data$boxpData$five[3]
    upHinge = data$boxpData$five[4]
    upEx = data$boxpData$five[5]

    ##Add bottom line, accounting for whether scaling is desired

    if(graphParams$scaleBot)
    {
        lines(data$bot$x, centre-(data$bot$y)/globals$maxBot, lwd=2, col=BOTLINE.COLOR, err=-1)
        botPeak<-min(centre-(data$bot$y)/globals$maxBot)
    }
    else
    {
        lines(data$bot$x, centre-(data$bot$y)/max(data$bot$y), lwd=2, col=BOTLINE.COLOR, err=-1)
        botPeak<-min(centre-(data$bot$y)/max(data$bot$y))
    }
    ##Add top line, accounting for whether scaling is desired
    if(graphParams$scaleTop)
    {
        lines(data$top$x, centre+(data$top$y)/globals$maxTop, col=TOPLINE.COLOR, lwd=2, lty=SOLID.LINE)
        topPeak<-max(centre+(data$top$y)/globals$maxTop)
    }
    else
    {

        lines(data$top$x, centre+(data$top$y)/max(data$top$y), col=TOPLINE.COLOR, lwd=2, lty=SOLID.LINE)
        topPeak<-max(centre+(data$top$y)/max(data$top$y))
    }

    ###Add boxplot ticks
    segments(lowHinge, centre, lowHinge, centre+graphParams$boxWidth, lty=SOLID.LINE)               ##Box wall, lower quartile
    segments(upHinge    , centre, upHinge, centre+graphParams$boxWidth, lty=SOLID.LINE)         ##Box wall, upper quartile
    segments(lowHinge, centre, lowEx, centre, lty = DASHED.LINE)        ##Lower whisker
    segments(upHinge, centre, upEx, centre, lty = DASHED.LINE)          ##Upper whisker
    segments(upEx, centre, upEx, centre+graphParams$whiskWidth, lty = DASHED.LINE)  #Upper Whisker Wall
    segments(lowEx, centre, lowEx, centre+graphParams$whiskWidth, lty = DASHED.LINE)    #Lower Whisker Wall

    ##Dot the median
    points(dat.median,centre)

    ##Add outlying values
    points(data$boxpData$outs, rep(centre, length(data$boxpData$outs)), pch = graphParams$outlierCharacter, col = graphParams$outlierColor)



    ###add max and min ticks along plot surface y-axis.
    segments(graphParams$rightBorder-graphParams$shift/2.5, botPeak, graphParams$rightBorder, botPeak, lty = SOLID.LINE, xpd=F)
    segments(graphParams$rightBorder-graphParams$shift/2.5, topPeak, graphParams$rightBorder, topPeak, lty = SOLID.LINE, xpd=F)
    segments(graphParams$leftBorder+graphParams$shift/2, centre, graphParams$leftBorder, centre, lty = 1, xpd=F)
    segments(graphParams$leftBorder+graphParams$shift/2, centre+1.25, graphParams$leftBorder, centre+1.25, lty = SOLID.LINE, xpd=F)
    segments(graphParams$leftBorder+graphParams$shift/2, centre-1.25, graphParams$leftBorder, centre-1.25, lty = SOLID.LINE, xpd=F)


    ###R and S-PLUS are highly divergent in their handling of text in the margins.  For this reason, the following section has been provided in two distinct components, each designed to cater to their specified language.
    ###R-Version
    RIGHT.EDGE.TEXT.SIZE=0.6
    LEFT.EDGE.TEXT.SIZE=0.8
    YLAB.TEXT.SIZE=0.65
    mtext(as.character(round(max(data$bot$y),2)), 4, at=botPeak, cex=RIGHT.EDGE.TEXT.SIZE)###Mark the largest value obtained by the bottom half of this shivplot on the y-axis
    mtext(as.character(round(max(data$top$y), 2)), 4, at=topPeak, cex=RIGHT.EDGE.TEXT.SIZE)###Mark the largest value obtained by the top half of this shivplot on the y-axis
    mtext(as.character(0), 2, at=centre, cex=0.8)#Add zero to the y-axis, left side

     ###Mark a tick for a common maximum scale if global scaling desired.  Otherwise, a value a fixed distance from each local maximum is marked.
     if(graphParams$scaleTop)
     {
        mtext(as.character(round(globals$maxTop*1.25, 2)), 2, at=centre+1.25, cex=LEFT.EDGE.TEXT.SIZE)
     }
     else
     {
         mtext(as.character(round(max(data$top$y)*1.25, 2)), 2, at=centre+1.25, cex=LEFT.EDGE.TEXT.SIZE)
     }
     if(graphParams$scaleBot)
     {
         mtext(as.character(round(globals$maxBot*1.25, 2)), 2, at=centre-1.25, cex=LEFT.EDGE.TEXT.SIZE)
     }
     else
     {
         mtext(as.character(round(max(data$bot$y)*1.25, 2)), 2, at=centre-1.25, cex=LEFT.EDGE.TEXT.SIZE)
     }
     ###Add labels to the top and bottom halves of the y-axis.
     mtext(data$bot$ylab, 2, at=centre-1, line=1, cex=YLAB.TEXT.SIZE, crt=90)
     mtext(data$top$ylab, 2, at=centre+1, line=1, cex=YLAB.TEXT.SIZE, crt=90)
}








###Function that adds shivplots to prepared plot surface.
shiv.plotItS<-function(data, globals, floor, plotColor, graphParams)

{
    #Constants
    ##Formats plot surface to have no y-axis, x-axis as normal
    DASHED.LINE = 2  ##Value for a dashed line
    SOLID.LINE = 1 ##Value for a solid line

    ###These values specify the margins within which each shivplot is contained.  The y-axis is divided into blocks, and each shiv plot is constructed within a 4(default) unit block.
    figFloor=floor
    centre = figFloor+2
    TOPLINE.COLOR = plotColor  ###by default, both the upper and lower halves of the shiv plot are drawn using the same color.  This can be changed here if desired.
    BOTLINE.COLOR = plotColor


    ###Recovering the five values required for boxplots
    lowEx = data$boxpData$five[1]
    lowHinge = data$boxpData$five[2]
    dat.median = data$boxpData$five[3]
    upHinge = data$boxpData$five[4]
    upEx = data$boxpData$five[5]

    ##Add bottom line, accounting for whether scaling is desired

    if(graphParams$scaleBot)
    {
        lines(data$bot$x, centre-(data$bot$y)/globals$maxBot, lwd=2, col=BOTLINE.COLOR, err=-1)
        botPeak<-min(centre-(data$bot$y)/globals$maxBot)
    }
    else
    {
        lines(data$bot$x, centre-(data$bot$y)/max(data$bot$y), lwd=2, col=BOTLINE.COLOR, err=-1)
        botPeak<-min(centre-(data$bot$y)/max(data$bot$y))
    }
    ##Add top line, accounting for whether scaling is desired
    if(graphParams$scaleTop)
    {
        lines(data$top$x, centre+(data$top$y)/globals$maxTop, col=TOPLINE.COLOR, lwd=2, lty=SOLID.LINE)
        topPeak<-max(centre+(data$top$y)/globals$maxTop)
    }
    else
    {

        lines(data$top$x, centre+(data$top$y)/max(data$top$y), col=TOPLINE.COLOR, lwd=2, lty=SOLID.LINE)
        topPeak<-max(centre+(data$top$y)/max(data$top$y))
    }

    ###Add boxplot ticks
    segments(lowHinge, centre, lowHinge, centre+graphParams$boxWidth, lty=SOLID.LINE)               ##Box wall, lower quartile
    segments(upHinge    , centre, upHinge, centre+graphParams$boxWidth, lty=SOLID.LINE)         ##Box wall, upper quartile
    segments(lowHinge, centre, lowEx, centre, lty = DASHED.LINE)        ##Lower whisker
    segments(upHinge, centre, upEx, centre, lty = DASHED.LINE)          ##Upper whisker
    segments(upEx, centre, upEx, centre+graphParams$whiskWidth, lty = DASHED.LINE)  #Upper Whisker Wall
    segments(lowEx, centre, lowEx, centre+graphParams$whiskWidth, lty = DASHED.LINE)    #Lower Whisker Wall

    ##Dot the median
    points(dat.median,centre)

    ##Add outlying values
    points(data$boxpData$outs, rep(centre, length(data$boxpData$outs)), pch = graphParams$outlierCharacter, col = graphParams$outlierColor)



    ###add max and min ticks along plot surface y-axis.
    segments(graphParams$rightBorder-graphParams$shift/2.5, botPeak, graphParams$rightBorder, botPeak, lty = SOLID.LINE, xpd=F)
    segments(graphParams$rightBorder-graphParams$shift/2.5, topPeak, graphParams$rightBorder, topPeak, lty = SOLID.LINE, xpd=F)
    segments(graphParams$leftBorder+graphParams$shift/2, centre, graphParams$leftBorder, centre, lty = 1, xpd=F)
    segments(graphParams$leftBorder+graphParams$shift/2, centre+1.25, graphParams$leftBorder, centre+1.25, lty = SOLID.LINE, xpd=F)
    segments(graphParams$leftBorder+graphParams$shift/2, centre-1.25, graphParams$leftBorder, centre-1.25, lty = SOLID.LINE, xpd=F)

    ###R and S-PLUS are highly divergent in their handling of text in the margins.  For this reason, the following section has been provided in two distinct components, each designed to cater to their specified language.
    ###S-Version

    RIGHT.EDGE.TEXT.SIZE=0.52
    LEFT.EDGE.TEXT.SIZE=0.6
    YLAB.TEXT.SIZE=0.6
    text(graphParams$rightBorder-.1*graphParams$shift, botPeak, round(max(data$bot$y),2), adj=0, cex=RIGHT.EDGE.TEXT.SIZE)###Mark the largest value obtained by the bottom half of this shivplot on the y-axis
    text(graphParams$rightBorder-.1*graphParams$shift, topPeak, round(max(data$top$y), 2), adj=0, cex=RIGHT.EDGE.TEXT.SIZE)###Mark the largest value obtained by the top half of this shivplot on the y-axis
    text(graphParams$leftBorder-.1*graphParams$shift, centre, 0, adj=1, cex=0.6)#Add zero to the y-axis, left side


    ###Mark a tick for a common maximum scale if global scaling desired.  Otherwise, a value a fixed distance from each local maximum is marked.
    if(graphParams$scaleTop)
    {
        text(graphParams$leftBorder-.1*graphParams$shift, centre+1.25, round(globals$maxTop*1.25, 2), adj=1, cex=LEFT.EDGE.TEXT.SIZE)
    }
    else
    {
        text(graphParams$leftBorder-.1*graphParams$shift, centre+1.25, round(max(data$top$y)*1.25, 2), adj=1, cex=LEFT.EDGE.TEXT.SIZE)
    }
    if(graphParams$scaleBot)
    {
        text(graphParams$leftBorder-.1*graphParams$shift, centre-1.25, round(globals$maxBot*1.25, 2), adj=1, cex=LEFT.EDGE.TEXT.SIZE)
    }
    else
    {
        text(graphParams$leftBorder-.1*graphParams$shift, centre-1.25, round(max(data$bot$y)*1.25, 2), adj=1, cex=LEFT.EDGE.TEXT.SIZE)
    }
    ###Add labels to the top and bottom halves of the y-axis.
    mtext(data$bot$ylab, 2, at=centre-1, line=4, cex=YLAB.TEXT.SIZE, srt=90)
    mtext(data$top$ylab, 2, at=centre+1, line=4, cex=YLAB.TEXT.SIZE, srt=90)    

}

###A couple example plugins.  These are the plugins described in the manuscript.
sdVsLocPluginR<-function(dataxx, LOCATION.FCN)
{
    ylab="Standard Deviation"###The label to use on the y-axis
    ylablong="Variability vs Intensity trend across arrays"
    x<-apply(dataxx, 1, LOCATION.FCN)
    y<-apply(dataxx, 1, sd)
    loe<-loess(y~x)
    y<-loe$fitt[order(x)]
    x<-sort(x)
    return(list(x=x, y=y, ylab=ylab, ylablong= ylablong))###return value must be a list consisting of x-values, y-values, and a label.
}



sdVsLocPluginS<-function(dataxx, LOCATION.FCN)
{
    ylab="Standard Deviation"###The label to use on the y-axis
    ylablong="Standard deviation across reps"
    x<-apply(dataxx, 1, LOCATION.FCN)
    y<-rowStdevs(dataxx)
    loe<-loess(y~x)
    y<-loe$fitt[order(x)]
    x<-sort(x)
    return(list(x=x, y=y, ylab=ylab, ylablong= ylablong))###return value must be a list consisting of x-values, y-values, and a label.
}

densityPlugin<-function(dataxx, LOCATION.FCN)
{
    ylab="Density"
    ylablong="Density of intensity means across arrays"
    x<-apply(dataxx, 1, LOCATION.FCN)
    temp<-density(x)
    return(list(x=temp$x, y=temp$y, ylab=ylab, ylablong= ylablong))
}

#@textFigs
#####################################################################################
#################################Figures for Manuscript##############################
###This function was only intended to be called from R.
exportAll<-function(outputD="D://", actualData=FALSE)
{
    outputD<-addslash(outputD)
    #ps.options(horizontal=FALSE, paper="special", width=7.187499, height=7.177082, pagecentre=T)
    ps.options(horizontal=FALSE, width = 7, height = 7)
    ###All plots export
    shiv.fig1<-function()
    {
        set.seed(123)
        norm<-rnorm(1000, 10, 1)
        bimode<-c(rnorm(500, 5, 1), rnorm(500, 15, 1))
        temp<-cbind(norm)
        temp<-rbind(temp, cbind(bimode))
        temp<-rbind(temp, cbind(RMAmeans))
        dimnames(temp)<-NULL
        types<-c(rep("a) normal", 1000), rep("b) bimodal", 1000), rep("c) RMA Microarray Data", length(RMAmeans)))
        temp<-data.frame(temp=rev(temp), types=rev(types))

        par(mfrow=c(3,3))
        boxplot(norm, horizontal=T, main = "\na) Boxplot", border="chocolate", xlab="magnitude")
        mtext("Normal Distribution Data", 3, 2, cex=1.2)
        boxplot(bimode, horizontal=T, main = "\nb) Boxplot", border="green", xlab="magnitude")
        mtext("Bimodal Distribution Data", 3, 2, cex=1.2)
        boxplot(RMAmeans, horizontal=T, main = "\nc) Boxplot", border="purple", xlab="RMA\n(mean intensity across replicate arrays)")
        mtext("RMA Microarray Data", 3, 2, cex=1.2)

        hist(norm, xlab="magnitude", ylab="frequency", density=-1, main="", border="chocolate")
        title(main="d) Histogram")
        hist(bimode, xlab="magnitude", ylab="frequency", density=-1, main="", border="green")
        title(main="e) Histogram")
        hist(RMAmeans, xlab="RMA\n(mean intensity across replicate arrays)", ylab="frequency", density=-1, main="", border="purple")
        title(main="f) Histogram")

        plot(density(norm), type="l", xlab="magnitude", ylab="density", main="", col="chocolate", lwd=2)
        title(main="g) Density Plot")
        plot(density(bimode), type="l", xlab="magnitude", ylab="density", main="", col="green", lwd=2)
        title(main="h) Density Plot")
        plot(density(RMAmeans), type="l", xlab="RMA\n(mean intensity across replicate arrays)", ylab="density", main="", col="purple", lwd=2)
        title(main="i) Density Plot")
    }

    shiv.fig2<-function()
    {
        par(mfrow=c(2,2))
        temp<-boxplot(MAS5means, horizontal=T, border="black", main="", xlab="MAS 5.0\n(mean intensity across replicate arrays)")
	points(temp$out, rep(1, length(temp$out)), col="gold")
        title("a)  Boxplot")
	  denstemp<-density(MAS5means)
        plot(density(MAS5means), type="n", main="", xlab="MAS 5.0\n(mean intensity across replicate arrays)", ylab="Density", ylim=c(-max(denstemp$y)*2, max(denstemp$y)*2))
        #lines(density(MAS5means), col="red", lwd=2)
	  lines(denstemp$x, denstemp$y, col="red", lwd=2)
        title("b)  Density Plot")
        plot(MAS5means, MAS5sds, main="", xlab="MAS 5.0\n(mean intensity across replicate arrays)", ylab="Standard Deviation")
        lines(sort(MAS5means), loess(MAS5sds~MAS5means)$fitt[order(MAS5means)], col="red", lwd=2)
        title("c)  Standard Deviation vs Mean")
        shivplot(MAS5data[,1:3], mean, c("MAS 5.0"), plotColors=c("red"),  outlierColor= "gold", T, T, main="d) Shivplot")
    }

    shiv.fig3<-function()
    {
        par(mfrow=c(4,3))
        boxplot(RMAmeans, horizontal=TRUE, border="purple", xlab="RMA\n(mean intensity across replicate arrays)", main="\n   Boxplot")
        mtext("RMA", 3, 2, cex=1.2)
        boxplot(MAS5means, horizontal=TRUE, border="red", xlab="MAS 5.0\n(mean intensity across replicate arrays)", main="\nBoxplot")
        mtext("MAS 5.0", 3, 2, cex=1.2)
        boxplot(dChipmeans, horizontal=TRUE, border="blue", xlab="dChip(v2004)\n(mean intensity across replicate arrays)", main="\nBoxplot")
        mtext("dChip(v2004)", 3, 2, cex=1.2)

        plot(density(RMAmeans), type="l", col="purple", lwd=2, xlab="RMA \n(mean intensity across replicate arrays)", ylab="Density", main="Density Plot")
        plot(density(MAS5means), type="l", col="red", lwd=2, xlab="MAS 5.0 \n(mean intensity across replicate arrays)", ylab="Density", main="Density Plot")
        plot(density(dChipmeans), type="l", col="blue", lwd=2, xlab="dChip(v2004)\n(mean intensity across replicate arrays)", ylab="Density", main="Density Plot")

        plot(RMAmeans, RMAsds, pch=".", xlab="RMA\n(mean intensity across replicate arrays)", ylab="Standard Deviation", main="Standard Deviation\n vs Mean")
        lines(sort(RMAmeans), loe1$fitt[order(RMAmeans)], lwd=2, col="purple")

        plot(MAS5means, MAS5sds, pch=".", xlab="MAS 5.0\n(mean intensity across replicate arrays)", ylab="Standard Deviation", main="Standard Deviation\n vs Mean")
        lines(sort(MAS5means), loe2$fitt[order(MAS5means)], lwd=2, col="red")

        plot(dChipmeans, dChipsds, pch=".", xlab="dChip(v2004)\n(mean intensity across replicate arrays)", ylab="Standard Deviation", main="Standard Deviation\n vs Mean")
        lines(sort(dChipmeans), loe3$fitt[order(dChipmeans)], lwd=2, col="blue")

        plot(RMAmeans, RMAsds, pch=".", type="n", xlab="RMA\n(mean intensity across replicate arrays)", ylab="Standard Deviation",main="Standard Deviation\n vs Mean, loess fit only")
        lines(sort(RMAmeans), loe1$fitt[order(RMAmeans)], lwd=2, col="purple")

        plot(MAS5means, MAS5sds, pch=".", type="n", xlab="MAS 5.0\n(mean intensity across replicate arrays)", ylab="Standard Deviation",main="Standard Deviation\n vs Mean, loess fit only")
        lines(sort(MAS5means), loe2$fitt[order(MAS5means)], lwd=2, col="red")

        plot(dChipmeans, dChipsds, pch=".", type="n", xlab="dChip(v2004)\n(mean intensity across replicate arrays)", ylab="Standard Deviation",main="Standard Deviation\n vs Mean, loess fit only")
        lines(sort(dChipmeans), loe3$fitt[order(dChipmeans)], lwd=2, col="blue")

    }

    shiv.fig4<-function()
    {
        par(mfrow=c(1,1))
        shivplot(data=dataBundle, locationFCN=mean, plotNames=c("RMA", "MAS 5.0", "dChip(v2004)"), plotColors=c("purple", "red", "blue"), outlierColor="orange", scaleTop=T, scaleBot=T, main="Comparative Shivplot")
    }

    shiv.fig5<-function()
    {
        par(mfrow=c(3,1))
        shivplot(data=dataBundle[[1]], locationFCN=mean, plotNames=c("RMA", "MAS 5.0", "dChipv2004")[1], plotColors="purple", outlierColor="orange", scaleTop=T, scaleBot=T, main="RMA\nSingle Shivplot")
        shivplot(data=dataBundle[[2]], locationFCN=mean, plotNames=c("RMA", "MAS 5.0", "dChipv2004")[2], plotColors="red", outlierColor="orange",scaleTop=T, scaleBot=T, main="MAS 5.0\nSingle Shivplot")
        shivplot(data=dataBundle[[3]], locationFCN=mean, plotNames=c("RMA", "MAS 5.0", "dChip (v2004)")[3], plotColors="blue", outlierColor="orange",scaleTop=T, scaleBot=T, main="dChip (v2004)\nSingle Shivplot")
    }


    dumpSite=outputD

    if(actualData)
    {
        cat("Trying to find cel files and dchip expression file in provided directory...\n")
    	  RMAdata<-as.matrix(read.table(paste(outputD, "RMAdata.csv", sep=""), sep=",", header=T, row.names=1))[,1:3]
	  MAS5data<-as.matrix(read.table(paste(outputD, "MASdata.csv", sep=""), sep=",", header=T, row.names=1))[,1:3]
        dChipdata<-as.matrix(read.table(paste(outputD, "dchipdata.csv", sep=""), sep=",", header=T, row.names=1))
    }

    else
    {
        cat("\nSimulating data, content of plots will differ from manuscript\n")
        MAS5data<-matrix(rnorm(22300*3, 10, 2), ncol=3)
        RMAdata<-matrix(runif(22300*3, 7.5, 15), ncol=3)
        dChipdata<-matrix(rgamma(22300*3, 10, 1), ncol=3)
    }

    RMAmeans<-rowMeans(RMAdata[,1:3])
    MAS5means<-rowMeans(MAS5data[,1:3])
    dChipmeans<-rowMeans(dChipdata[,1:3])

    RMAsds<-apply(RMAdata[,1:3], 1, sd)
    MAS5sds<-apply(MAS5data[,1:3], 1, sd)
    dChipsds<-apply(dChipdata[,1:3], 1, sd)


    loe1<- loess(RMAsds~RMAmeans)
    loe2<- loess(MAS5sds~MAS5means)
    loe3<- loess(dChipsds~dChipmeans)
    dataBundle<-list(rma=RMAdata[,1:3], mas=MAS5data[,1:3], dChip=dChipdata[,1:3])

    plotFuncs<-list(
    shiv.fig1,
    shiv.fig2,
    shiv.fig3,
    shiv.fig4,
    shiv.fig5
    )


    for(i in 1:4)
    {
        #postscript(paste(dumpSite, "Figure", i, ".ps", sep=""))
	  pdf(paste(dumpSite, "Figure", i, ".pdf", sep=""), width=7, height=7)
        plotFuncs[[i]]()
        dev.off()
    }

    ps.options(horizontal=FALSE, width = 5, height = 8.5)
    #postscript(paste(dumpSite, "Figure5.ps", sep=""), )
    pdf(paste(dumpSite, "Figure5.pdf", sep=""), width=5, height=8.5)
    shiv.fig5()
    dev.off()
    ps.options(horizontal=FALSE, width = 7, height = 7)
    cat(paste("Output has been saved to ", outputD, "...\n", sep=""))
}


#@examples:
###########################Sample plots###############################


shivSampleS<-function(outputD=NULL)
{
    graphsheet(Name="shivsample")
    set.seed(123)
    data=list(
    standardNormal=matrix(rnorm(3000, 5, 2), ncol=3),
    bimodal=rbind(matrix(rnorm(1500, 3, 2),ncol=3), matrix(rnorm(1500, 7, 2), ncol=3)),
    uniform=matrix(runif(3000, 0, 10), ncol=3)
    )
    names=c("Standard Normal", "Bimodal", "Uniform")

    par(mfrow=c(3,3))
    for(i in 1:3){
        boxplot(rowMeans(data[[i]]), xlab="mean intensity across replicate arrays", horizontal=T, main=names[i], boxcol=i+1)
    }

    for (i in 1:3){
        dens<-density(rowMeans(data[[i]]))
        plot(dens, ylab="density", xlab="mean intensity across replicate arrays", main="", type="n")
        lines(dens, col=i+1)
    }

    for(i in 1:3){
        curr<-data[[i]]
        means<-rowMeans(curr)
        sds<-apply(curr, 1, stdev)
        loe <- loess(sds ~ means)
        loefits <- loe$fitt ##Scale values to fit in plot area
        plot(means, sds, ylab="standard deviation across reps", xlab="mean intensity across replicate arrays", type="n")
        points(means, sds, pch=".", col=i+1)
        lines(sort(means), loefits[order(means)], lwd=2, col=5)
    }

    par(mfrow=c(1,3))
    for(i in 1:3)
        shivplot(data=data[[i]], locationFCN=mean, plotNames=c("Standard Normal", "Bimodal", "Uniform")[i], plotColors=c(2,3,4)[i], outlierColor=5, scaleTop=T, scaleBot=T, main="Single Shivplot")
    par(mfrow=c(1,1))


    shivplot(data= data, locationFCN=mean, plotNames=c("Standard Normal", "Bimodal", "Uniform"), plotColors=c(2,3,4), outlierColor=5, scaleTop=T, scaleBot=T, main="Comparative Shivplot")

    graphsheet(Name="multishivsample")
    multishivplot(data= data, plotNames=c("Standard Normal", "Bimodal", "Uniform"), plotColors=c(2,3,4), outlierColor=5, titles="Multi-sample Shivplot")

    if(!is.null(outputD)) {
        guiSave("GraphSheet", Name="shivsample", FileName=paste(outputD, "shivsample.sgr", sep=""))
        guiClose("GraphSheet", name = "shivsample")
        guiSave("GraphSheet", Name="multishivsample", FileName=paste(outputD, "multishivsample.sgr", sep=""))
        guiClose("GraphSheet", name = "multishivsample")
        print(paste("Output saved in ...", outputD))
    }
}

shivSampleR<-function(outputD=NULL)
{
    if(!is.null(outputD)){
        outputD<-addslash(outputD)
    }
    ps.options(horizontal=FALSE, width = 7, height = 7)
    set.seed(123)  #For reproducibility
    data.dir=outputD  #Desired output directory

    #This following declaration illustrates how to construct a list that is suitable to be passed in as the "data" argument for shivplot.
    data=list(
    standardNormal=matrix(rnorm(3000, 5, 2), ncol=3),
    bimodal=rbind(matrix(rnorm(1500, 3, 2),ncol=3), matrix(rnorm(1500, 7, 2), ncol=3)),
    uniform=matrix(runif(3000, 0, 10), ncol=3)
    )

    names=c("Standard Normal", "Bimodal", "Uniform")##names of distributions for plot titles

    if(!is.null(outputD)){
        postscript(paste(data.dir, "example1.ps", sep=""))
    }
    else {
        x11()
    }
    par(mfrow=c(3,3))##First illustration is a 3x3 plot showing all the component subplots of a shivplot on the various distributions.
    for(i in 1:3){
        boxplot(rowMeans(data[[i]]), xlab="mean intensity across replicate arrays", horizontal=T, main=names[i], col=i+1)
    }

    for (i in 1:3){
        plot(density(rowMeans(data[[i]])), ylab="density", xlab="mean intensity across replicate arrays", main="", col=i+1)
    }
    for(i in 1:3)
    {
        curr<-data[[i]]
        means<-rowMeans(curr)
        sds<-apply(curr, 1, sd)
        loe <- loess(sds ~ means)
        loefits <- loe$fitt ##Scale values to fit in plot area
        plot(means, sds, pch=".", ylab="standard deviation across reps", xlab="mean intensity across replicate arrays", col=i+1)
        lines(sort(means), loefits[order(means)], lwd=2, col="violet")
    }
    if(!is.null(outputD)){
        dev.off()##Save as a postscript in the designated output directory
        postscript(paste(data.dir, "example2.ps", sep=""))
    }
    else
    {
        x11()
    }


    par(mfrow=c(1,3))  ##Next illustration is of the three shivplots that represent the multiple subplots shown previously.  In this figure, the shivplots are not scaled relative to one another.
    for(i in 1:3)
        shivplot(data=data[[i]], locationFCN=mean, plotNames=c("Standard Normal", "Bimodal", "Uniform")[i], plotColors=c(2,3,4)[i], outlierColor=5, scaleTop=F, scaleBot=F, main="Single Shivplot")
    if(!is.null(outputD)){
        dev.off()##Save as a postscript in the designated output directory
        postscript(paste(data.dir, "example3.ps", sep=""))
    }
    else {
        x11()
    }

    shivplot(data = data, locationFCN=mean, plotNames=c("Standard Normal", "Bimodal", "Uniform"), plotColors=c(2,3,4), outlierColor=5, scaleTop=T, scaleBot=T, main="Comparative Shivplot")
    if(!is.null(outputD)){
        dev.off()##Save as a postscript in the designated output directory
	  postscript(paste(data.dir, "example4.ps", sep=""))        

    }
    else {
	x11()
    }

    multishivplot(data = data, plotNames=c("Standard Normal", "Bimodal", "Uniform"), plotColors=c(2,3,4), outlierColor=5, titles="Multi-sample Shivplot")
    if(!is.null(outputD)){
        dev.off()##Save as a postscript in the designated output directory
	  cat(paste("Output has been saved to ", outputD, "...\n", sep=""))
    }


}

#@multiS
###################################################################################
#shivplot argument list:
#data,                  #either a matrix with rows representing subjects and columns representing observations, or a list composed of a number of such matrices
#locationFCN=mean,          #The function to use to obtain an estimate of location across observations.  The default measure (and the one used throughout this document) is the mean.
#plotNames,                 #Names to attribute to each individual shivplot.  The names will be located above the plots, in the plot frame.  The number of plotNames must match the length of the list used in the 'data' argument.  If these subnames are not desired, enter the string "" for each name.
#plotColors=c(1),           #What colors to use for each created plot.  Expects a vector with a number of entries equal to the length of the list in the 'data' argument.  S-PLUS users must specify a vector of integers for colors, whereas R users may either use numbers, calls to the function "colors()", or strings representing color names.
#outlierColor=2,            #The color to use to mark observations that qualify as "outliers" under default boxplot specifications.
#FUN1=sdVsLocPlugin,            #The function to use to compute the line on the top half of the shivplot.  By default, this function relates standard deviations across replicates to the selected measure of centre.  These functions must follow a special (but simple) 'plugin' format, see the provided code for details.
#FUN2=densityPlugin,            #The 'plugin' function to use to compute the line in the lower half of the shivplot.
#xLab="mean intensity across replicate arrays",    #The label for the shared x-axis for a single or comparative boxplot.  Similar to xlab from standard R/S-PLUS plotting.  Be sure to change this to reflect the choice of location function (locationFCN).
#boxWidth=.2,               ##
#whiskWidth=.35,             ##Controls the length of the ticks drawn to denote the "box" region from a boxplot, and the "whisker" region denoting the limit past which observations are deemed to be outliers.  Negative values cause the ticks to be reflected about the x-axis.
#outlierCharacter=-1,          #The plotting character to use to mark outliers.  The default (-1) is the system default.

#multishivplotS<-function(data, ltys, cols, names, titles)
multishivplotS<-function(data, plotColors, lineTypes, plotNames, titles="", FUN1=sdVsLocPlugin, FUN2=densityPlugin, locationFCN=mean, boxWidth=.3, whiskWidth=.5, outlierColor=2, outlierCharacter=-1,xLab="mean intensity across replicate arrays" )
{

	names<-plotNames
	ncomponents = length(data)
	colors<-1:ncomponents
	if(!missing(plotColors))
	{
		colors[1:min(length(plotColors), ncomponents)]<-plotColors
	}
	if(missing(lineTypes))
	{
		lineTypes<-rep(1, ncomponents)
	}
	else{
		if(length(lineTypes)==1)
		{
			lineTypes<-rep(lineTypes, ncomponents)
		}
	}
	if(missing(plotNames))
	{
		plotNames<-rep("", ncomponents)
	}
	locs<-NULL
	sds<-NULL
	topx<-NULL
	topy<-NULL
	botx<-NULL
	boty<-NULL
	fives<-NULL
	toplab<-NULL
	botlab<-NULL

	for(i in 1:ncomponents)
	{
		curr<-data[[i]]
		locs<-cbind(locs, apply(curr, 1, locationFCN))
		#sds<-cbind(sds, rowStDevs(curr)) 
		#top<-lowess(locs[,i], sds[,i])
		top<-FUN1(curr, locationFCN)
		topx<-cbind(topx, top$x)
		topy<-cbind(topy, top$y)
		toplab<-top$ylablong

		bot<-FUN2(curr, locationFCN)
		botx<-cbind(botx, bot$x)
		boty<-cbind(boty, bot$y)
		botlab<-bot$ylablong
		fives<-cbind(fives, sort(boxplot(locs[,i], plot=F)$stats))

	}
	dimnames(locs)[[2]]<-names   #####
	topymax<-max(topy)
	botymax<-max(boty)
	meanmax<-max(locs)
	meanmin<-min(locs)
	#frame()
	#graphsheet()
	plot(c(0,1), c(0,1), type="n", xaxt="n", yaxt="n", ylab="", xlab="", bty="n")
	par(mfrow=c(1,1))
	subplot(fun=
		{
			plot(1:10, 1:10, xlim=c(meanmin, meanmax), ylim=c(0, topymax), type="n", xaxt="n", ylab = "", xlab="")
			mtext(toplab, 2, line=1.9, cex=.75)
			for(i in 1:ncomponents)
			{
				lines(topx[,i], topy[,i], col=colors[i], lty=lineTypes[i], lwd=2)
			}
		},
	x=c(0, 1), y= c(.583, 1))

	title(titles)

	subplot(fun=
		{
			plot(1:10, 1:10, xlim=c(meanmin, meanmax), ylim=c(0.8, (length(data)+.2)), type="n", axes=F, frame.plot=T, ylab="", xlab="")         #####
			box(1)
			axis(2, labels = rep("", length(data)), las=1, at = 1:length(data))
			#mtext(names, side=2, line= 1.5, at=1:length(data))
			for(i in 1:ncomponents)
			{
				newstyle=F###toggle
				if(newstyle)
				{
					points(fives[,i], rep(i, nrow(fives)), pch="I", cex=c(1,1,1,1,1), col= colors[i])
	
					lines(fives[c(2,4),i], rep(i, 2), col=colors[i], lty=lineTypes[i], lwd=2)
				}
				else
				{
					#print(fives[3,i])	
					points(fives[3,i], i, col=colors[i])
					#lines(rep(fives[2, i], 2), c(i, i+0.2), col=colors[i], lty=lineTypes[i])
					#lines(rep(fives[4, i], 2), c(i, i+0.2), col=colors[i], lty=lineTypes[i])
					#lines(rep(fives[5, i], 2), c(i, i+0.35), col=colors[i], lty=lineTypes[i])
					#lines(rep(fives[1, i], 2), c(i, i+0.35), col=colors[i], lty=lineTypes[i])

					solid= 1
					lines(rep(fives[2, i], 2), c(i, i+boxWidth), col=colors[i], lty=solid)
					lines(rep(fives[4, i], 2), c(i, i+boxWidth), col=colors[i], lty=solid)
					lines(rep(fives[5, i], 2), c(i, i+whiskWidth), col=colors[i], lty=solid)
					lines(rep(fives[1, i], 2), c(i, i+whiskWidth), col=colors[i], lty=solid)

					lines(fives[c(1,2), i], c(i, i), col=colors[i], lty=lineTypes[i])
					lines(fives[c(4,5), i], c(i, i), col=colors[i], lty=lineTypes[i])
					points(locs[,i][locs[,i]<fives[1,i]], rep(i, sum(locs[,i]<fives[1,i])), col=outlierColor, pch=outlierCharacter)
					points(locs[,i][locs[,i]>fives[5,i]], rep(i, sum(locs[,i]>fives[5,i])), col=outlierColor, pch=outlierCharacter)


				}
			}
		}, 
		x=c(0, 1), y= c(.416, .583) )
		

	subplot(fun=
		{	
			plot(1, xlim=c(meanmin, meanmax), ylim=c(-botymax,0), axes=F, frame.plot=T, ylab="")
			box(1)
			mtext(xLab, 1, line=1.9, cex=.75)
			mtext(botlab, 2, line=1.9, cex=.75)
			axis(1)
			axis(2, labels = as.character(-pretty(c(-botymax, 0), 5)), at = pretty(c(-botymax, 0), 5), srt=90)
			key(corner = c(1,0), text= list(a = rev(names)),  line = list(lty = rev(lineTypes), col = rev(colors), lwd = 1), col=rev(colors), border=T)
			for(i in 1:ncomponents)
			{
				lines(botx[,i], -boty[,i], col=colors[i], lty= lineTypes[i], lwd=2)
			}
		},
		x=c(0, 1), y= c(0, .416) )
		return()
}

#@multiR
################################################################################
#shivplot argument list:
#data,                  #either a matrix with rows representing subjects and columns representing observations, or a list composed of a number of such matrices
#locationFCN=mean,          #The function to use to obtain an estimate of location across observations.  The default measure (and the one used throughout this document) is the mean.
#plotNames,                 #Names to attribute to each individual shivplot.  The names will be located above the plots, in the plot frame.  The number of plotNames must match the length of the list used in the 'data' argument.  If these subnames are not desired, enter the string "" for each name.
#plotColors=c(1),           #What colors to use for each created plot.  Expects a vector with a number of entries equal to the length of the list in the 'data' argument.  S-PLUS users must specify a vector of integers for colors, whereas R users may either use numbers, calls to the function "colors()", or strings representing color names.
#outlierColor=2,            #The color to use to mark observations that qualify as "outliers" under default boxplot specifications.
#scaleTop=T,                #When multiple shivplots are drawn in tandem (denoted through this manuscript as "comparative shivplots"), the user has the option to enforce a common scale for the "top" half of the y-axis.  This argument expects a Boolean variable toggling this feature.
#scaleBot=T,                #As above, with bottom half of y-axis
#FUN1=sdVsLocPlugin,            #The function to use to compute the line on the top half of the shivplot.  By default, this function relates standard deviations across replicates to the selected measure of centre.  These functions must follow a special (but simple) 'plugin' format, see the provided code for details.
#FUN2=densityPlugin,            #The 'plugin' function to use to compute the line in the lower half of the shivplot.
#xLab="mean intensity across replicate arrays",     #The label for the shared x-axis for a single or comparative boxplot.  Similar to xlab from standard R/S-PLUS plotting.  Be sure to change this to reflect the choice of location function (locationFCN).
#boxWidth=.3,               ##
#whiskWidth=.5,             ##Controls the length of the ticks drawn to denote the "box" region from a boxplot, and the "whisker" region denoting the limit past which observations are deemed to be outliers.  Negative values cause the ticks to be reflected about the x-axis.
#outlierCharacter=-1,          #The plotting character to use to mark outliers.  The default (-1) is the system default.



multishivplotR<-function(data, plotColors, lineTypes, plotNames, titles="", FUN1=sdVsLocPlugin, FUN2=densityPlugin,locationFCN=mean, boxWidth=.3, whiskWidth=.5, outlierColor="gold", outlierCharacter=-1, xLab="mean intensity across replicate arrays")
{
	oldpar<-par(no.readonly=TRUE)
	names<-plotNames
	ncomponents = length(data)
	colors<-1:ncomponents
	colors<-1:ncomponents
	if(!missing(plotColors))
	{
		colors[1:min(length(plotColors), ncomponents)]<-plotColors
	}
	if(missing(lineTypes))
	{
		lineTypes<-rep(1, ncomponents)
	}
	else{
		if(length(lineTypes)==1)
		{
			lineTypes<-rep(lineTypes, ncomponents)
		}
	}
	if(missing(plotNames))
	{
		plotNames<-rep("", ncomponents)
	}
	locs<-NULL
	sds<-NULL
	topx<-NULL
	topy<-NULL
	botx<-NULL
	boty<-NULL
	fives<-NULL
	toplab<-NULL
	botlab<-NULL
	for(i in 1:ncomponents)
	{
		curr<-data[[i]]
		locs<-cbind(locs, apply(curr, 1, locationFCN))
		#sds<-cbind(sds, apply(curr, 1, sd))
		#top<-lowess(locs[,i], sds[,i])
		top<-FUN1(curr, locationFCN)
		topx<-cbind(topx, top$x)
		topy<-cbind(topy, top$y)
		toplab<-top$ylablong
		bot<-FUN2(curr, locationFCN)
		botx<-cbind(botx, bot$x)
		boty<-cbind(boty, bot$y)
		botlab<-bot$ylablong
		#fives<-cbind(fives, fivenum(locs[,i]))
		fives<-cbind(fives, boxplot(locs[,i], plot=F)$stats)

	}
	colnames(locs)<-names   #####
	topymax<-max(topy)
	botymax<-max(boty)
	meanmax<-max(locs)
	meanmin<-min(locs)

	layout(matrix(c(1,2,3, 0), 4, 1, byrow=T), heights = c(5, 1, 5, .2))
	par(mai=c(0, .76875, .3, 0.39375))
	plot(1:10, xlim=c(meanmin, meanmax), ylim=c(0, topymax), type="n", xaxt="n", ylab = toplab)
	for(i in 1:ncomponents)
	{
		lines(topx[,i], topy[,i], col=colors[i], lty=lineTypes[i], lwd=2)
	}
	title(titles)
	par(mai=c(0, .76875, 0, 0.39375))
	plot(1:10, xlim=c(meanmin, meanmax), ylim=c(0.8, (length(data)+.2)), type="n", axes=F, frame.plot=T, ylab="")         #####
	####axis(2, labels = names, las=1, at = 1:length(data))
	for(i in 1:ncomponents)
	{
		newstyle=F###toggle
		if(newstyle)
		{
			points(fives[,i], rep(i, nrow(fives)), pch="I", cex=c(1,1,1,1,1), col= colors[i])

			lines(fives[c(2,4),i], rep(i, 2), col=colors[i], lty=lineTypes[i], lwd=2)
		}
		else
		{
			shift = ncomponents+1-i
			solid=1
			points(fives[3,i], shift, col=colors[i])
			lines(rep(fives[2, i], 2), c(shift, shift+boxWidth), col=colors[i], lty=solid)
			lines(rep(fives[4, i], 2), c(shift, shift+boxWidth), col=colors[i], lty=solid)
			lines(rep(fives[5, i], 2), c(shift, shift+whiskWidth), col=colors[i], lty=solid)
			lines(rep(fives[1, i], 2), c(shift, shift+whiskWidth), col=colors[i], lty=solid)
			lines(fives[c(1,2), i], c(shift, shift), col=colors[i], lty=lineTypes[i])
			lines(fives[c(4,5), i], c(shift, shift), col=colors[i], lty=lineTypes[i])
			points(locs[,i][locs[,i]<fives[1,i]], rep(shift, sum(locs[,i]<fives[1,i])), col=outlierColor, pch=outlierCharacter)
			points(locs[,i][locs[,i]>fives[5,i]], rep(shift, sum(locs[,i]>fives[5,i])), col=outlierColor, pch=outlierCharacter)
		}
	}
	par(mai=c(.3, .76875, 0, 0.39375))
	plot(1, xlim=c(meanmin, meanmax), ylim=c(-botymax,0), axes=F, frame.plot=T, ylab=botlab, xlab = xLab)
	mtext("mean intensity across replicate arrays", 1, line=1.9, cex=.75)
	axis(1)
	axis(2, labels = as.character(-axTicks(2)), at = axTicks(2))
	for(i in 1:ncomponents)
	{
		lines(botx[,i], -boty[,i], col=colors[i], lty= lineTypes[i], lwd=2)
	}
	legend(x="bottomright", legend=names, col=colors, lty=lineTypes, lwd=1)
	par(oldpar)
}

############################################################


#@support
{
if(match("language", names(version), nomatch=F)) {
	shiv.lang = "R"
}
else{
	shiv.lang="S"
}
}

{
if(shiv.lang=="R")
{
	shiv.plotIt<-shiv.plotItR
	multishivplot<-multishivplotR
	stdev<-sd
	sdVsLocPlugin<-sdVsLocPluginR
	shivSample<-shivSampleR
}
else   ##shiv.lang=="S"
{
	
	shiv.plotIt<-shiv.plotItS
	multishivplot<-multishivplotS
	sd<-stdev
	sdVsLocPlugin<-sdVsLocPluginS
	shivSample<-shivSampleS
}
}

addslash<-function(word)
{
    tail<-substring(word, nchar(word), nchar(word))
    if (!is.element(tail, c("\\", "/")))
    {
        return(paste(word, "//", sep=""))
    }
    else
    {
        return(word)
    }
}

####################################################################################

shivHelp<-function() {
cat("Okay, shivplot should be ready to go.
In order to invoke any of the samples illustrated below, just type (or
copy and paste) any of the instructions given at left
at the command prompt.
---------
Sample Commands:
1.\tshivSample(outputD=\"D://\")\t\t#Produce a series of example plots
\t\t\t\t\t\tIf a path is provided, graphical output is saved to files.
---------
2.\tshivplot(
\tdata=matrix(rnorm(300, 0, 1), ncol=3),
\tplotNames=\"Silly Test\",
\tplotColors=13)\t\t\t#Build a shivplot out of a 100x3 random data set
\t\t\t\t\t\t(100 criteria, 3 replicate observations)
---------
3.\texportAll(outputD=\"D://\")\t\tCreate plots from the manuscript
\t\t\t\t\t\tusing simulated data, output to D://
\t\t\t\t\t\t(Only R is supported for the generation
\t\t\t\t\t\tof manuscript plots)
---------
To see this message at any time, type 'shivHelp()' at the command prompt
i.e.
 >  shivHelp()
...
---------
")

}

shivHelp()