.packageName <- "gregmisc"
Args <- function(name, sort.args=FALSE)
{
  a <- formals(get(as.character(substitute(name))))
  if(is.null(a))
    return(NULL)
  arg.labels <- names(a)
  arg.values <- as.character(a)
  char <- sapply(a, is.character)
  arg.values[char] <- paste("\"", arg.values[char], "\"", sep="")

  if(sort.args)
  {
    ord <- order(arg.labels)
    if(any(arg.labels == "..."))
      ord <- c(ord[-1], length(ord))
    arg.labels <- arg.labels[ord]
    arg.values <- arg.values[ord]
  }

  output <- data.frame(value=I(arg.values), row.names=arg.labels)
  print(output, right=FALSE)

  invisible(output)
}
# $Log: CrossTable.R,v $
# Revision 1.4  2004/05/25 02:57:01  warnes
# Updates from Mark Schwartz.
#
# Revision 1.4  2004/03/01 schwartz
#
# - Thanks to Gabor Grothendieck for providing ideas and
#   comments on the following updates
# - Modified function to allow a single vector to be
#   tabulated in 'x'
# - Added 'chisq' argument to disable the generation and
#   printing of the chi-square tests, producing a simple
#   table with no tests. Default value = FALSE.
# - Removed restriction on minimum 2 x 2 tables (ie. can do a
#   1 x n table, but no stats will be generated
# - Added code to create generic row/col names if argument 'x'
#   is a matrix without dimnames
# - Added  'max.width' argument to enable wrapped output for
#   1 x n tables. If n > max.width, output will wrap to next
#   line
# - Cleaned up list return code
# - Separated functions to print m x n table and 1 x n vector
#   to enable flexibility in vector output formats and to faciliate
#   future creation and separation of table generation and print methods.


# Revision 1.3  2003/01/30 21:41:35  warnes
#
# - Removed argument 'correct' and now print separate corrected values
#   for 2 x 2 tables.
# - Added arguments 'prop.r', 'prop.c' and 'prop.t' to toggle printing
#   of row, col and table percentages. Default is TRUE.
# - Added argument 'fisher' to toggle fisher exact test. Default is
#   FALSE.
# - Added McNemar test to statistics and argument 'mcnemar' to toggle
#   test. Default is FALSE.
# - Added code to generate an invisible return list containing table
#   counts, proportions and the results of the appropriate statistical tests.
#
# Revision 1.3  2003/01/26
# Updates from Marc Schwartz:
#
# - Removed argument 'correct' and now print separate corrected values for 2 x 2 tables.
# - Added arguments 'prop.r', 'prop.c' and 'prop.t' to toggle printing of row, col
#   and table percentages. Default is TRUE.
# - Added argument 'fisher' to toggle fisher exact test. Default is FALSE.
# - Added McNemar test to statistics and argument 'mcnemar' to toggle test.
#   Default is FALSE.
# - Added code to generate an invisible return list containing table counts, proportions
#   and the results of the appropriate statistical tests.
#
#
# $Log: CrossTable.R,v $
# Revision 1.4  2004/05/25 02:57:01  warnes
# Updates from Mark Schwartz.
#
# Revision 1.3  2003/01/30 21:41:35  warnes
#
# - Removed argument 'correct' and now print separate corrected values
#   for 2 x 2 tables.
# - Added arguments 'prop.r', 'prop.c' and 'prop.t' to toggle printing
#   of row, col and table percentages. Default is TRUE.
# - Added argument 'fisher' to toggle fisher exact test. Default is
#   FALSE.
# - Added McNemar test to statistics and argument 'mcnemar' to toggle
#   test. Default is FALSE.
# - Added code to generate an invisible return list containing table
#   counts, proportions and the results of the appropriate statistical tests.
#
# Revision 1.2  2002/11/04 14:13:57  warnes
#
# - Moved fisher.test() to after table is printed, so that table is
#   still printed in the event that fisher.test() results in errors.
# ------------------------------------------------------------------------
#


CrossTable <- function (x, y,
                        digits = 3,
                        max.width = 5,
                        expected = FALSE,
                        prop.r = TRUE,
                        prop.c = TRUE,
                        prop.t = TRUE,
                        chisq = FALSE,
                        fisher = FALSE,
                        mcnemar = FALSE)
{

  # Ensure that max.width >= 1
  if (max.width < 1)
    stop("max.width must be >= 1")
  # Set 'x' vector flag
  vector.x <- FALSE
  # Ensure that if (expected), a chisq is done
  if (expected)
    chisq <- TRUE
  
  if (missing(y))
  {
    # is x a vector?
    if (is.null(dim(x)))
    {
      # Remove any unused factor levels
      x <- factor(x)
      t <- t(as.matrix(table(x)))
      vector.x <- TRUE
    }  
    # is x a matrix?
    else if (length(dim(x) == 2))
    {
      if(any(x < 0) || any(is.na(x)))
        stop("all entries of x must be nonnegative and finite")
      
      # Add generic dimnames if required
      # check each dimname separately, in case user has defined one or the other
      if (is.null(rownames(x)))
        rownames(x) <- paste("[", 1:nrow(x), ",]", sep = "")
      if (is.null(colnames(x)))
        colnames(x) <- paste("[,", 1:ncol(x), "]", sep = "")

      t <- x
    }
    else
      stop("x must be either a vector or a 2 dimensional matrix, if y is not given")
  }
  else
  {
    if(length(x) != length(y))
      stop("x and y must have the same length")

    # Create Titles for Table From Vector Names
    RowData <- deparse(substitute(x))
    ColData <- deparse(substitute(y))

    # Remove unused factor levels from vectors
    x <- factor(x)
    y <- factor(y)

    # Generate table
    t <- table(x, y)
  }

  # if t is not at least a 2 x 2, do not do stats
  # even if any set to TRUE. Do not do col/table props
  if (any(dim(t) < 2))
  {  
    prop.c <- prop.r <- chisq <- expected <- fisher <- mcnemar <- FALSE
  }  
  
  # Generate cell proportion of row
  CPR <- prop.table(t, 1)

  # Generate cell proportion of col
  CPC <- prop.table(t, 2)

  # Generate cell proportion of total
  CPT <- prop.table(t)

  # Generate summary counts
  GT <- sum(t)
  RS <- rowSums(t)
  CS <- colSums(t)

  # Column and Row Total Headings
  ColTotal <- "Column Total"
  RowTotal <- "Row Total"

  # Set consistent column widths based upon dimnames and table values
  CWidth <- max(digits + 2, c(nchar(t), nchar(dimnames(t)[[2]]), nchar(RS), nchar(CS), nchar(RowTotal)))
  RWidth <- max(c(nchar(dimnames(t)[[1]]), nchar(ColTotal)))

  # Adjust first column width if Data Titles present
  if (exists("RowData"))
    RWidth <- max(RWidth, nchar(RowData))

  # Create row separators
  RowSep <- paste(rep("-", CWidth + 2), collapse = "")
  RowSep1 <- paste(rep("-", RWidth + 1), collapse = "")
  SpaceSep1 <- paste(rep(" ", RWidth), collapse = "")
  SpaceSep2 <- paste(rep(" ", CWidth), collapse = "")

  # Create formatted Names
  FirstCol <- formatC(dimnames(t)[[1]], width = RWidth, format = "s")
  ColTotal <- formatC(ColTotal, width = RWidth, format = "s")
  RowTotal <- formatC(RowTotal, width = CWidth, format = "s")

  # Perform Chi-Square Tests
  # Needs to be before the table output, in case (expected = TRUE)
  if (chisq)
  {  
    if (all(dim(t) == 2))
      CSTc <- chisq.test(t, correct = TRUE)

    CST <- chisq.test(t, correct = FALSE)
  }

  # Default print function for n x m table
  print.CrossTable.default <- function()
  {
    # Print Column headings
    if (exists("RowData"))
    {
      cat(SpaceSep1, "|", ColData, "\n")
      cat(formatC(RowData, width = RWidth, format = "s"), formatC(dimnames(t)[[2]], width = CWidth, format = "s"),
          RowTotal, sep = " | ", collapse = "\n")
    }
    else
      cat(SpaceSep1, formatC(dimnames(t)[[2]], width = CWidth, format = "s"), RowTotal,
          sep = " | ", collapse = "\n")

    cat(RowSep1, rep(RowSep, ncol(t) + 1), sep = "|", collapse = "\n")

    # Print table cells
    for (i in 1:nrow(t))
    {
      cat(FirstCol[i], formatC(c(t[i, ], RS[i]), width = CWidth), sep = " | ", collapse = "\n")

      if (expected)
        cat(SpaceSep1, formatC(CST$expected[i, ], digits = digits, format = "f", width = CWidth), SpaceSep2,
            sep = " | ", collapse = "\n")

      if (prop.r)
        cat(SpaceSep1, formatC(c(CPR[i, ], RS[i] / GT), width = CWidth, digits = digits, format = "f"),
            sep = " | ", collapse = "\n")
        
      if (prop.c)
        cat(SpaceSep1, formatC(CPC[i, ], width = CWidth, digits = digits, format = "f"), SpaceSep2,
            sep = " | ", collapse = "\n")
      
      if (prop.t)
        cat(SpaceSep1, formatC(CPT[i, ], width = CWidth, digits = digits, format = "f"), SpaceSep2,
            sep = " | ", collapse = "\n")

      cat(RowSep1, rep(RowSep, ncol(t) + 1), sep = "|", collapse = "\n")
    }

    # Print Column Totals
    cat(ColTotal, formatC(c(CS, GT), width = CWidth), sep = " | ", collapse = "\n")
    
    if (prop.c)
      cat(SpaceSep1, formatC(CS / GT, width = CWidth, digits = digits, format = "f"), SpaceSep2,
          sep = " | ", collapse = "\n")

    cat(RowSep1, rep(RowSep, ncol(t) + 1), sep = "|", collapse = "\n")
  }

  # Print function for 1 x n vector
  print.CrossTable.vector <- function()
  {
    if (length(t) > max.width)
    {
      # set breakpoints for output based upon max.width
      final.row <- length(t) %% max.width
      max <- length(t) - final.row
      # Define breakpoint indices for each row
      start <- seq(1, max, max.width)
      end <- start + (max.width - 1)
      # Add final.row if required
      if (final.row > 0)
      {
         start <- c(start, end[length(end)] + 1)
         end <- c(end, end[length(end)] + final.row)
      }
    }
    else
    {
      # Each value printed horizontally in a single row
      start <- 1
      end <- length(t)
    }

    SpaceSep3 <- paste(SpaceSep2, " ", sep = "")
    
    for (i in 1:length(start))
    {
      # print column labels
      cat(SpaceSep2, formatC(dimnames(t)[[2]][start[i]:end[i]], width = CWidth, format = "s"),
          sep = " | ", collapse = "\n")

      cat(SpaceSep3, rep(RowSep, (end[i] - start[i]) + 1), sep = "|", collapse = "\n")
      cat(SpaceSep2, formatC(t[, start[i]:end[i]], width = CWidth), sep = " | ", collapse = "\n")
      cat(SpaceSep2, formatC(CPT[, start[i]:end[i]], width = CWidth, digits = digits, format = "f"),
          sep = " | ", collapse = "\n")
      cat(SpaceSep3, rep(RowSep, (end[i] - start[i]) + 1), sep = "|", collapse = "\n")
      cat("\n\n")
    }
  }

  # Print Cell Layout

  cat(rep("\n", 2))
  cat("   Cell Contents\n")

  cat("|-----------------|\n")
  cat("|               N |\n")
  if (expected)
    cat("|      Expected N |\n")
  if (prop.r)
    cat("|   N / Row Total |\n")
  if (prop.c)
    cat("|   N / Col Total |\n")
  if (prop.t)
    cat("| N / Table Total |\n")
  cat("|-----------------|\n")
  cat(rep("\n", 2))
  cat("Total Observations in Table: ", GT, "\n")
  cat(rep("\n", 2))

  if (!vector.x)
    print.CrossTable.default()
  else
    print.CrossTable.vector()

  # Print Statistics
  if (chisq)
  {
    cat(rep("\n", 2))
    cat("Statistics for All Table Factors\n\n\n")

    cat(CST$method,"\n")
    cat("------------------------------------------------------------\n")
    cat("Chi^2 = ", CST$statistic, "    d.f. = ", CST$parameter, "    p = ", CST$p.value, "\n\n")

    if (all(dim(t) == 2))
    {
      cat(CSTc$method,"\n")
      cat("------------------------------------------------------------\n")
      cat("Chi^2 = ", CSTc$statistic, "    d.f. = ", CSTc$parameter, "    p = ", CSTc$p.value, "\n")
    }
  }

  # Perform McNemar tests
  if (mcnemar)
  {
    McN <- mcnemar.test(t, correct = FALSE)
    cat(rep("\n", 2))
    cat(McN$method,"\n")
    cat("------------------------------------------------------------\n")
    cat("Chi^2 = ", McN$statistic, "    d.f. = ", McN$parameter, "    p = ", McN$p.value, "\n\n")

    if (all(dim(t) == 2))
    {
      McNc <- mcnemar.test(t, correct = TRUE)
      cat(McNc$method,"\n")
      cat("------------------------------------------------------------\n")
      cat("Chi^2 = ", McNc$statistic, "    d.f. = ", McNc$parameter, "    p = ", McNc$p.value, "\n")
    }
  }


  # Perform Fisher Tests
  if (fisher)
  {
    cat(rep("\n", 2))
    FTt <- fisher.test(t, alternative = "two.sided")

    if (all(dim(t) == 2))
    {
      FTl <- fisher.test(t, alternative = "less")
      FTg <- fisher.test(t, alternative = "greater")
    }

    cat("Fisher's Exact Test for Count Data\n")
    cat("------------------------------------------------------------\n")

    if (all(dim(t) == 2))
    {
      cat("Sample estimate odds ratio: ", FTt$estimate, "\n\n")

      cat("Alternative hypothesis: true odds ratio is not equal to 1\n")
      cat("p = ", FTt$p.value, "\n")
      cat("95% confidence interval: ", FTt$conf.int, "\n\n")

      cat("Alternative hypothesis: true odds ratio is less than 1\n")
      cat("p = ", FTl$p.value, "\n")
      cat("95% confidence interval: ", FTl$conf.int, "\n\n")

      cat("Alternative hypothesis: true odds ratio is greater than 1\n")
      cat("p = ", FTg$p.value, "\n")
      cat("95% confidence interval: ", FTg$conf.int, "\n\n")
    }
    else
    {
      cat("Alternative hypothesis: two.sided\n")
      cat("p = ", FTt$p.value, "\n")
    }
  }

  cat(rep("\n", 2))
  
  # Create list of results for invisible()

  CT <- list(t = t, prop.row = CPR, prop.col = CPC, prop.tbl = CPT)
  
  if (any(chisq, fisher, mcnemar))
  {
    if (all(dim(t) == 2))
    {
      if (chisq)
        CT <- c(CT, list(chisq = CST, chisq.corr = CSTc))

      if (fisher)
        CT <- c(CT, list(fisher.ts = FTt, fisher.tl = FTl, fisher.gt = FTg))

      if (mcnemar)
        CT <- c(CT, list(mcnemar = McN, mcnemar.corr = McNc))
    }
    else
    {
      if (chisq)
        CT <- c(CT, list(chisq = CST))

      if (fisher)
        CT <- c(CT, list(fisher.ts = FTt))

      if (mcnemar)
        CT <- c(CT, list(mcnemar = McN))
    }
  }  

  # return list(CT)
  invisible(CT)
}

# $Id: RSCompat.S,v 1.8 2003/04/04 13:58:59 warnes Exp $
#
# $Log: RSCompat.S,v $
# Revision 1.8  2003/04/04 13:58:59  warnes
#
# - Replace 'T' with 'TRUE'
#
# Revision 1.7  2003/03/07 15:48:35  warnes
#
# - Minor changes to code to allow the package to be provided as an
#   S-Plus chapter.
#
# Revision 1.6  2003/01/02 15:42:00  warnes
# - Add nlevels function.
#
# Revision 1.5  2002/03/20 03:44:32  warneg
# - Added definition of is.R function.
#
# - Added boxplot.formula
#
# Revision 1.4  2002/02/05 02:20:07  warneg
#
# - Fix typo that caused code meant to run only under S-Plus to run
#   under R, causing problems.
#
# Revision 1.3  2001/12/19 22:45:44  warneg
# - Added code for %in%.
#
# Revision 1.2  2001/09/18 14:15:44  warneg
#
# Release 0.3.2
#
# Revision 1.1  2001/09/01 19:19:13  warneg
#
# Initial checkin.
#
#
# Code necessary for contrast.lm, boxplot.n to work in S-Plus
 
if(!exists("is.R") || !is.R() )
  {
    is.R <- function() FALSE 
    
    getOption <- function(...) options(...)
    
    if(!exists("parent.frame")) parent.frame <- sys.parent
    
    colnames <- function (x, do.NULL = TRUE, prefix = "col") 
      {
        dn <- dimnames(x)
        if (!is.null(dn[[2]])) 
          dn[[2]]
        else {
          if (do.NULL) 
            NULL
          else paste(prefix, seq(length = NCOL(x)), sep = "")
        }
      }
    
    rownames <- function (x, do.NULL = TRUE, prefix = "row") 
      {
        dn <- dimnames(x)
        if (!is.null(dn[[1]])) 
          dn[[1]]
        else {
          if (do.NULL) 
            NULL
          else paste(prefix, seq(length = NROW(x)), sep = "")
        }
      }
    
    "rownames<-" <- function (x, value) 
      {
        dn <- dimnames(x)
        ndn <- names(dn)
        dn <- list(value, if (!is.null(dn)) dn[[2]])
        names(dn) <- ndn
        dimnames(x) <- dn
        x
      }
    
    "colnames<-" <- function (x, value) 
      {
       dn <- dimnames(x)
       ndn <- names(dn)
       dn <- list(if (!is.null(dn)) dn[[1]], value)
       names(dn) <- ndn
       dimnames(x) <- dn
       x
     }
    
   # from the MASS library by Venables & Ripley 
    ginv <- function (X, tol = sqrt(.Machine$double.eps))
      {
        if (length(dim(X)) > 2 || !(is.numeric(X) || is.complex(X)))
          stop("X must be a numeric or complex matrix")
        if (!is.matrix(X))
          X <- as.matrix(X)
        Xsvd <- svd(X)
        if (is.complex(X))
          Xsvd$u <- Conj(Xsvd$u)
        Positive <- Xsvd$d > max(tol * Xsvd$d[1], 0)
        if (all(Positive)) Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u))
        else if (!any(Positive)) array(0, dim(X)[2:1])
        else Xsvd$v[, Positive] %*% ((1/Xsvd$d[Positive]) * t(Xsvd$u[, Positive]))
      }
    
    
    "format.pval" <- 
      function (pv, digits = max(1, getOption("digits") - 2),
                eps = .Machine$double.eps, 
                na.form = "NA") 
        {
          if ((has.na <- any(ina <- is.na(pv)))) 
            pv <- pv[!ina]
          r <- character(length(is0 <- pv < eps))
          if (any(!is0)) {
            rr <- pv <- pv[!is0]
            expo <- floor(log10(pv))
            fixp <- expo >= -3 | (expo == -4 & digits > 1)
            if (any(fixp)) 
              rr[fixp] <- format(pv[fixp], dig = digits)
            if (any(!fixp)) 
              rr[!fixp] <- format(pv[!fixp], dig = digits)
            r[!is0] <- rr
          }
          if (any(is0)) {
            digits <- max(1, digits - 2)
            if (any(!is0)) {
              nc <- max(nchar(rr))
              if (digits > 1 && digits + 6 > nc) 
                digits <- max(1, nc - 7)
              sep <- if (digits == 1 && nc <= 6) 
                ""
              else " "
            }
            else sep <- if (digits == 1) 
              ""
            else " "
            r[is0] <- paste("<", format(eps, digits = digits), sep = sep)
          }
          if (has.na) {
            rok <- r
            r <- character(length(ina))
            r[!ina] <- rok
            r[ina] <- na.form
          }
          r
        }
    
    "%in%" <- function (x, table)  match(x, table, nomatch = 0) > 0
 
    strwidth   <-  function(...)
      {
        par("cin")[1] / par("fin")[1] * (par("usr")[2] - par("usr")[1])
      }
    
    strheight <-  function(...)
      {
        par("cin")[2] / par("fin")[2] * (par("usr")[4] - par("usr")[3])
      }
    
    boxplot.formula <- function(x, data = sys.parent(), ..., ask = TRUE)
      {
        if(!inherits(x, "formula"))
          x <- as.formula(x)
        
        mf <- model.frame(x, data, na.action = function(z) 	z)
        if(length(names(mf)) > 2) 
          stop("boxplot.formula only accepts models with 1 predictor")
        
        resp <- attr(attr(mf, "terms"), "response")
        class(mf) <- NULL
        y <- mf[[resp]]
        x <- mf[[-resp]]
        xlab <- names(mf)[-resp]
        ylab <- names(mf)[resp]
	
        boxplot(split(y, x), xlab = xlab, ylab = ylab, ...) 
      }

    nlevels <- function(x) length(levels(x))

    NULL
    
  }
# $Id: aggregate.table.R,v 1.3 2002/09/23 13:59:30 warnes Exp $
#
# $Log: aggregate.table.R,v $
# Revision 1.3  2002/09/23 13:59:30  warnes
# - Modified all files to include CVS Id and Log tags.
#
# Revision 1.2  2002/02/21 21:45:01  warneg
#
# - Fixed bug where row and column labels didn't always correspond to the
# contents.  This only occured when a factor was used for by1 or by2 and
# the factors levels weren't in the default sort order.
#
# Revision 1.1  2002/02/20 22:10:08  warneg
#
# New function.
#
#

aggregate.table <- function(x, by1, by2, FUN=mean, ... )
  {
    if(!is.factor(by1)) by1 <- as.factor(by1)
    if(!is.factor(by2)) by2 <- as.factor(by2)
    
    ag <- aggregate(x, by=list(by1,by2), FUN=FUN, ... )
    tab <- matrix( nrow=nlevels(by1), ncol=nlevels(by2) )
    dimnames(tab) <- list(levels(by1),levels(by2))

    for(i in 1:nrow(ag))
      tab[ as.character(ag[i,1]), as.character(ag[i,2]) ] <- ag[i,3]
    tab
  }
# $Id: balloonplot.R,v 1.4 2003/04/04 13:49:27 warnes Exp $
#
# $Log: balloonplot.R,v $
# Revision 1.4  2003/04/04 13:49:27  warnes
#
# - Change occurences of 'T' to 'TRUE'
#
# Revision 1.3  2003/03/08 16:22:47  warnes
#
# - Added parameters for rotation of and amount of space allocated for
#   the row and column labels.
#
# Revision 1.2  2003/01/20 18:35:47  warnes
#
# - Updated balloonplot help page.
#
# Revision 1.1  2003/01/03 21:34:06  warnes
# - Initial checkin of balloonplot functions and documentation.
#
#

balloonplot <- function(x,...)
  UseMethod("balloonplot",x)



balloonplot.table <- function(x, xlab, ylab, zlab, ... )
  {
    obj <- x
    tmp <- as.data.frame(x)
    x <- tmp[,1]
    y <- tmp[,2]
    z <- tmp[,3]
    tableflag <- TRUE
      
    if(missing(xlab)) xlab <- names(dimnames(obj))[1] 
    if(missing(ylab)) ylab <- names(dimnames(obj))[2]
    if(missing(zlab)) zlab <- "Freq"

    balloonplot.default(x, y, z, xlab=xlab, ylab=ylab, zlab=zlab, ...)
  }



balloonplot.default <- function(x,y,z,
                                xlab=deparse(substitute(x)),
                                ylab=deparse(substitute(y)),
                                zlab=deparse(substitute(z)),
                                dotsize=2/max(strwidth(19),strheight(19)),
                                dotchar=19,
                                dotcolor="skyblue",
                                main,
                                label=TRUE,
                                label.digits=2,
                                scale.method=c("volume","diameter"),
                                colsrt=par("srt"),
                                rowsrt=par("srt"),
                                colmar=1,
                                rowmar=1,
                                ... )
{

  scale.method <- match.arg(scale.method)
  
  if(missing(main))
    {
      if(scale.method=="volume")
        main <- paste("Balloon Plot for ", xlab," by ", ylab,
                      ".\nArea is proportional to ", zlab, ".", sep='')
      else
        main <- paste("Balloon Plot for ", xlab," by ", ylab,
                      ".\nDiameter is proportional to ", zlab, ".", sep='')
      }
      

    
  x <- as.factor(x)
  y <- as.factor(y)
  z <- as.numeric(z)

  if( any(z < 0 ) )
    warning("z value(s) below zero detected. These will have zero size.")

  scale <- function(X, min=0, max=16, scale.method)
    {
      #X <- ( X -min(X) )
      if(scale.method=="volume")
        {
          X[X<0] <- 0
          X <- sqrt(X)
        }
      
      X <- min + (X/max(X, na.rm=TRUE) * (max - min) )
      X
    }
  
  plot(x=as.numeric(x),
       y=nlevels(y) - as.numeric(y) + 1,
       cex=scale(z, max=dotsize, scale.method=scale.method),
       pch=dotchar, # plot character
       col=dotcolor, # dot color
       xlab=xlab,
       ylab=ylab,
       xaxt="n", # no x axis lables
       yaxt="n", # no y axis lables
       bty="n",  # no box around the plot
       xlim=c(0.5-rowmar,nlevels(x)+0.5), # extra space on either end of plot
       ylim=c(0.5,nlevels(y)+1.5+colmar-1),  # so dots don't cross into margins,
       ...
     )

  # add text labels
  par(adj=0.5)
  text(x=1:nlevels(x),
       y=nlevels(y)+(colmar/2)+0.5,
       labels=levels(x),
       srt=colsrt)

  text(y=nlevels(y):1,
       x=-rowmar/2+0.5,
       labels=levels(y),
       srt=rowsrt)
  
  # add borders between cells
  abline(v=(0:nlevels(x)+0.5))
  abline(h=(0:nlevels(y)+0.5))
  
  # annotate with actual values
  if(label)
    text(x=as.numeric(x),     # as.numeric give numeric values
         y=nlevels(y) - as.numeric(y) + 1, 
         labels=format(z, digits=label.digits),       # label value
         col="black", # textt color
         )
  
  # put a nice title
  title(main=main)
}

# $Id: bandplot.R,v 1.3 2002/04/09 00:51:28 warneg Exp $
#
# $Log: bandplot.R,v $
# Revision 1.3  2002/04/09 00:51:28  warneg
#
# Checkin for version 0.5.3
#
# Revision 1.2  2002/02/20 20:06:45  warneg
#
# - Generalized to allow specification of the multiples of the standard deviation levels to be plotted (0=mean, 1=1 sd, ..).
# - Now (invisibly) returnes computed smooths.
#
# Revision 1.1  2001/08/31 23:41:58  warneg
#
# Used wrong comment character (% instead of #) in header.  Fixed.
#
#

bandplot  <-  function(x,y,
                       ..., 
                       add=FALSE,
                       sd=c(-2:2),
                       sd.col=c("lightblue","blue","red",
                                 "blue","lightblue"),
                       method="frac", width=1/5,
                       n=50
                       )
  {

    if(length(sd.col)<length(sd)) sd <-rep(sd.col,length=length(sd))
    
    if(!add)
      {
        m <- match.call(expand.dots = FALSE)
        m$width  <- m$add  <- m$sd  <- m$sd.col  <- NULL
        m$method <- m$n <- NULL
        m[[1]] <- as.name("plot")
        mf <- eval(m, parent.frame())
      }

    x  <- as.numeric(x)
    y  <- as.numeric(y)
    
    CL <- function(x,sd)
      if(length(x)<2)
        return( NA )
      else
        mean(x)+sd*sqrt(var(x))

    sdplot <- function(S, COL)
                  {
                    where <- wapply(x,y,CL,sd=S,width=width,method=method,n=n)
                    lines(where,col=COL,...)
                    where
                  }

    stdevs <- list()
    for( i in 1:length(sd) )
      stdevs[[as.character(sd[i])]] <- sdplot( sd[i], sd.col[i] )

    invisible( stdevs )
  }
# $Log: barplot2.R,v $
# Revision 1.7  2003/12/01 15:55:50  warnes
#
# - Follow patches applied to barplot() in base.
#
# Revision 1.6  2003/04/22 15:34:25  warnes
# Update from Marc Schwartz, modified by Gregory Warnes:
#
# -  Modified dim() checks for 'ci.l' and 'ci.u' against 'height'
#    to remove R v1.7.0 if() based error msgs for vector conditions.
#
# Revision 1.5  2003/01/30 21:43:05  warnes
#
# - Added argument 'add' to allow for the addition of a barplot to an
#   existing graphic. Default is FALSE
#
# revision 1.4  2003/01/02 16:09:46 warnes
#
#- Changed assignment statements that used "=" to "<-" to avoid syntax
#  errors in older versions of the S language.
#
# Revision 1.3  2002/11/04 14:21:40  warnes
# Updates from Marc Schwartz:
#
# - Updated underlying code to be based upon the new barplot() in R v1.6.1
# - This now uses the 'axis.lty' and 'border' arguments
# - In R v1.6.0, R Core introduced a new function called axTicks().
#   This is an R equivalent of the C code for CreateAtVector in plot.c.
#   This now enables me to get the axis tick mark positions consistently
#   when the 'height' related axis is either linear or log.  Thus, I can
#   now have consistent tick marks and can plot grid lines in either
#   situation.  If 'plot.grid = TRUE' and 'grid.inc' is specified, then
#   I still use pretty() to determine the tick marks and lines.
# - This code now depends on R 1.6.0 or later.
#
# Revision 1.2  2002/10/11 18:02:15  warnes
#
# - Fixed log scale errors in legend() call
#
# Revision 1.1  2002/09/23 13:38:53  warnes
#
# - Added CrossTable() and barplot2() code and docs contributed by Marc Schwartz.
# - Permit combinations() to be used when r>n provided repeat.allowed=TRUE
# - Bumped up version number
#

barplot2 <- function(height, ...) UseMethod("barplot2")

barplot2.default <-
function(height, width = 1, space = NULL, names.arg = NULL,
       legend.text = NULL, beside = FALSE, horiz = FALSE,
       density = NULL, angle = 45,
       col = heat.colors(NR), prcol = NULL, border = par("fg"),
       main = NULL, sub = NULL, xlab = NULL, ylab = NULL,
       xlim = NULL, ylim = NULL, xpd = TRUE, log = "",
       axes = TRUE, axisnames = TRUE,
       cex.axis = par("cex.axis"), cex.names = par("cex.axis"),
       inside = TRUE, plot = TRUE, axis.lty = 0,
       plot.ci = FALSE, ci.l = NULL, ci.u = NULL,
       ci.color = "black", ci.lty = "solid", ci.lwd = 1,
       plot.grid = FALSE, grid.inc = NULL,
       grid.lty = "dotted", grid.lwd = 1, grid.col = "black",
       add = FALSE, ...)
{
    if (!missing(inside)) .NotYetUsed("inside", error = FALSE)# -> help(.)

    if (missing(space))
      space <- if (is.matrix(height) && beside) c(0, 1) else 0.2

    space <- space * mean(width)

    if (plot && axisnames && missing(names.arg))
      names.arg <- if(is.matrix(height)) colnames(height) else names(height)

    if (is.vector(height))
    {
      height <- cbind(height)
      beside <- TRUE
    }
    else if (is.array(height) && (length(dim(height)) == 1))
    {
      height <- rbind(height)
      beside <- TRUE
    }
    else if (!is.matrix(height))
      stop("`height' must be a vector or a matrix")

    if(is.logical(legend.text))
      legend.text <-
        if(legend.text && is.matrix(height)) rownames(height)

    # Check for log scales
    logx <- FALSE
    logy <- FALSE

    if (log != "")
    {
      if (any(grep("x", log)))
        logx <- TRUE

      if (any(grep("y", log)))
        logy <- TRUE
    }

    # Cannot "hatch" with rect() when log scales used
    if ((logx || logy) && !is.null(density))
      stop("Cannot use shading lines in bars when log scale is used")

    NR <- nrow(height)
    NC <- ncol(height)

    if (beside)
      {
        if (length(space) == 2)
          space <- rep.int(c(space[2], rep.int(space[1], NR - 1)), NC)
        
        width <- rep(width, length = NR * NC)
      }
    else
      {
        width <- rep(width, length = NC)
      }

    delta <- width / 2
    w.r <- cumsum(space + width)
    w.m <- w.r - delta
    w.l <- w.m - delta

    #if graphic will be stacked bars, do not plot ci
    if (!beside && (NR > 1) && plot.ci)
      plot.ci = FALSE

    # error check ci arguments
    if (plot && plot.ci)
    {
      if ((missing(ci.l)) || (missing(ci.u)))
        stop("confidence interval values are missing")

      if (is.vector(ci.l))
        ci.l <- cbind(ci.l)
      else if (is.array(ci.l) && (length(dim(ci.l)) == 1))
        ci.l <- rbind(ci.l)
      else if (!is.matrix(ci.l))
        stop("`ci.l' must be a vector or a matrix")

      if (is.vector(ci.u))
        ci.u <- cbind(ci.u)
      else if (is.array(ci.u) && (length(dim(ci.u)) == 1))
        ci.u <- rbind(ci.u)
      else if (!is.matrix(ci.u))
        stop("`ci.u' must be a vector or a matrix")

      if ( any(dim(height) != dim(ci.u) ) )
        stop("'height' and 'ci.u' must have the same dimensions.")
      else if ( any( dim(height) != dim(ci.l) ) )
        stop("'height' and 'ci.l' must have the same dimensions.")
    }

    if (beside)
      w.m <- matrix(w.m, nc = NC)

    # check height/ci.l if using log scale to prevent log(<=0) error
    # adjust appropriate ranges and bar base values
    if ((logx && horiz) || (logy && !horiz))
    {
      if (min(height) <= 0)
        stop("log scale error: at least one 'height' value <= 0")

      if (plot.ci && (min(ci.l) <= 0))
        stop("log scale error: at least one lower c.i. value <= 0")

      if (logx && !is.null(xlim) && (xlim[1] <= 0))
        stop("log scale error: 'xlim[1]' <= 0")

      if (logy && !is.null(ylim) && (ylim[1] <= 0))
        stop("'log scale error: 'ylim[1]' <= 0")

      # arbitrary adjustment to display some of bar for min(height) or min(ci.l)
      if (plot.ci)
        rectbase <- rangeadj <- (0.9 * min(c(height, ci.l)))
      else
        rectbase <- rangeadj <- (0.9 * min(height))

      # if axis limit is set to < above, adjust bar base value
      # to draw a full bar
      if (logy && !is.null(ylim))
        rectbase <- ylim[1]
      else if (logx && !is.null(xlim))
        rectbase <- xlim[1]

      # if stacked bar, set up base/cumsum levels, adjusting for log scale
      if (!beside)
        height <- rbind(rectbase, apply(height, 2, cumsum))

      # if plot.ci, be sure that appropriate axis limits are set to include range(ci)
      lim <-
        if (plot.ci)
          c(height, ci.l, ci.u)
        else
          height
    }
    else
    {
      # Use original bar base value
      rectbase <- 0

      # if stacked bar, set up base/cumsum levels
      if (!beside)
        height <- rbind(rectbase, apply(height, 2, cumsum))

      # if plot.ci, be sure that appropriate axis limits are set to include range(ci)
      lim <-
        if (plot.ci)
          c(height, ci.l, ci.u)
        else
          height

      # use original range adjustment factor
      rangeadj <- (-0.01 * lim)
    }

    # define xlim and ylim, adjusting for log-scale if needed
    if (horiz)
    {
      if (missing(xlim)) xlim <- range(rangeadj, lim, na.rm=TRUE)
      if (missing(ylim)) ylim <- c(min(w.l), max(w.r))
    }
    else
    {
      if (missing(xlim)) xlim <- c(min(w.l), max(w.r))
      if (missing(ylim)) ylim <- range(rangeadj, lim, na.rm=TRUE)
    }

    if(plot) ##-------- Plotting :
    {
      opar <-
        if (horiz)
          par(xaxs = "i", xpd = xpd)
        else
          par(yaxs = "i", xpd = xpd)

      on.exit(par(opar))

      # If add = FALSE open new plot window
      # else allow for adding new plot to existing window
      if (!add)
      {
        plot.new()
        plot.window(xlim, ylim, log = log, ...)
      }

      # Set plot region coordinates
      usr <- par("usr")

      # adjust par("usr") values if log scale(s) used
      if (logx)
      {
        usr[1] <- 10 ^ usr[1]
        usr[2] <- 10 ^ usr[2]
      }

      if (logy)
      {
        usr[3] <- 10 ^ usr[3]
        usr[4] <- 10 ^ usr[4]
      }

      # if prcol specified, set plot region color
      if (!missing(prcol))
        rect(usr[1], usr[3], usr[2], usr[4], col = prcol)

      # if plot.grid, draw major y-axis lines if vertical or x axis if horizontal
      # R V1.6.0 provided axTicks() as an R equivalent of the C code for
      # CreateAtVector.  Use this to determine default axis tick marks when log
      # scale used to be consistent when no grid is plotted.
      # Otherwise if grid.inc is specified, use pretty()

      if (plot.grid)
      {
        par(xpd = FALSE)

        if (is.null(grid.inc))
        {
          if (horiz)
          {
            grid <- axTicks(1)
            abline(v = grid, lty = grid.lty, lwd = grid.lwd, col = grid.col)
          }
          else
          {
            grid <- axTicks(2)
            abline(h = grid, lty = grid.lty, lwd = grid.lwd, col = grid.col)
          }
        }
        else
        {
          if (horiz)
          {
            grid <- pretty(xlim, n = grid.inc)
            abline(v = grid, lty = grid.lty, lwd = grid.lwd, col = grid.col)
          }
          else
          {
            grid <- pretty(ylim, n = grid.inc)
            abline(h = grid, lty = grid.lty, lwd = grid.lwd, col = grid.col)
          }
        }

         par(xpd = xpd)
      }

      xyrect <- function(x1,y1, x2,y2, horizontal = TRUE, ...)
      {
        if(horizontal)
          rect(x1,y1, x2,y2, ...)
        else
          rect(y1,x1, y2,x2, ...)
      }

      if (beside)
        xyrect(rectbase, w.l, c(height), w.r, horizontal=horiz,
               angle = angle, density = density, col = col, border = border)
      else
      {
        for (i in 1:NC)
          xyrect(height[1:NR, i], w.l[i], height[-1, i], w.r[i],
                 horizontal=horiz,  angle = angle, density = density,
                 col = col, border = border)
      }

      if (plot.ci)
      {
        # CI plot width = barwidth / 2
        ci.width = width / 4

        if (horiz)
        {
          segments(ci.l, w.m, ci.u, w.m, col = ci.color, lty = ci.lty, lwd = ci.lwd)
          segments(ci.l, w.m - ci.width, ci.l, w.m + ci.width, col = ci.color, lty = ci.lty, lwd = ci.lwd)
          segments(ci.u, w.m - ci.width, ci.u, w.m + ci.width, col = ci.color, lty = ci.lty, lwd = ci.lwd)
        }
        else
        {
          segments(w.m, ci.l, w.m, ci.u, col = ci.color, lty = ci.lty, lwd = ci.lwd)
          segments(w.m - ci.width, ci.l, w.m + ci.width, ci.l, col = ci.color, lty = ci.lty, lwd = ci.lwd)
          segments(w.m - ci.width, ci.u, w.m + ci.width, ci.u, col = ci.color, lty = ci.lty, lwd = ci.lwd)
        }
      }

      if (axisnames && !is.null(names.arg)) # specified or from {col}names
      {
        at.l <-
        if (length(names.arg) != length(w.m))
        {
          if (length(names.arg) == NC) # i.e. beside (!)
            colMeans(w.m)
          else
            stop("incorrect number of names")
        }
        else w.m

        axis(if(horiz) 2 else 1, at = at.l, labels = names.arg, lty = axis.lty, cex.axis = cex.names, ...)
      }

      if(!is.null(legend.text))
      {
        legend.col <- rep(col, length = length(legend.text))

        if((horiz & beside) || (!horiz & !beside))
        {
          legend.text <- rev(legend.text)
          legend.col <- rev(legend.col)
          density <- rev(density)
          angle <- rev(angle)
        }

        # adjust legend x and y values if log scaling in use
        if (logx)
          legx <- usr[2] - ((usr[2] - usr[1]) / 10)
        else
          legx <- usr[2] - xinch(0.1)

        if (logy)
          legy <- usr[4] - ((usr[4] - usr[3]) / 10)
        else
          legy <- usr[4] - yinch(0.1)

        legend(legx, legy,
        legend = legend.text, angle = angle, density = density,
        fill = legend.col, xjust = 1, yjust = 1)
      }

      title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)

      # if axis is to be plotted, adjust for grid "at" values
      if (axes)
      {
        if(plot.grid)
          axis(if(horiz) 1 else 2, at = grid, cex.axis = cex.axis, ...)
        else
          axis(if(horiz) 1 else 2, cex.axis = cex.axis, ...)
      }

      invisible(w.m)

    }

    else w.m
}
# $Id: boxplot.n.R,v 1.4 2002/09/23 13:59:30 warnes Exp $
#
# $Log: boxplot.n.R,v $
# Revision 1.4  2002/09/23 13:59:30  warnes
# - Modified all files to include CVS Id and Log tags.
#
#

boxplot.n  <- function( ..., top=FALSE, shrink=1.0, textcolor=NULL )
  {
    boxcall <- match.call()           # get call
    boxcall$top  <- boxcall$shrink  <- boxcall$textcolor  <- NULL
    boxcall[[1]]  <- as.name("boxplot")
    
    if(is.R())
      {
        box <- eval(boxcall, parent.frame())
        mids <- 1:length(box$n)
      }
    else
      {
        mids <- eval(boxcall, parent.frame())
        boxcall$plot <- FALSE
        box <- eval(boxcall, parent.frame())
      }
    
    if(top)
      {
        where  <- par("usr")[4]
        adj  <- c(0.5,1)
      }
    else
      {
        where  <- par("usr")[3]
        adj  <- c(0.5,0)
      }
    tcex <- par("cex")
    par(cex=shrink*tcex)
    
    if(is.R())
      text( x=mids, y=where, labels=paste("n=",box$n,sep=""), adj=adj, 
			      col=textcolor)
    else
      {
        if( is.null(textcolor) ) 
          textcolor <- 1
        space <- ifelse(top, -1, 1) * par("1em")[2] / 2
        
        text( x=mids, y=where + space, labels=paste("n=",box$n,sep=""), adj=adj[1], 
             col=textcolor)
        
      }
    
    par(cex=tcex)
    
    invisible(box)
  }

# $Id: capture.R,v 1.4 2004/01/21 04:31:25 warnes Exp $
#
# $Log: capture.R,v $
# Revision 1.4  2004/01/21 04:31:25  warnes
# - Mark sprint() as depreciated.
# - Replace references to sprint with capture.output()
# - Use match.arg for halign and valign arguments to textplot.default.
# - Fix textplot.character so that a vector of characters is properly
#   displayed. Previouslt, character vectors were plotted on top of each
#   other.
#
# Revision 1.3  2003/11/10 22:11:13  warnes
#
# - Add files contributed by Arni Magnusson
#   <arnima@u.washington.edu>. As well as some of my own.
#
# Revision 1.2  2003/04/04 13:45:21  warnes
#
# - Allow optional arguments to sprint to be passed to print
#
# Revision 1.1  2003/04/02 22:28:32  warnes
#
# - Added file 'capture.R' containing capture() and sprint().
#
#

capture <- function( expression, collapse="\n")
  {
    warning("Depreciated.  Use capture.output(...) from base instead.")

    resultText <- capture.output( expression )

    return( paste( c(resultText, ""), collapse=collapse, sep="" ) )
    # the reason for c(result, "") is so that we get the line
    # terminator on the last line of output.  Otherwise, it just shows
    # up between the lines.
  }


sprint <- function(x,...)
  {
    warning("Depreciated.  Use capture.output(print(...)) from base instead.")
    capture(print(x,...))
  }
# $Id: ci.R,v 1.8 2003/08/07 03:49:54 warnes Exp $
#
# $Log: ci.R,v $
# Revision 1.8  2003/08/07 03:49:54  warnes
# - Fixed incorrect denominator in standard error for mean in ci.default.
#
# Revision 1.7  2002/09/23 13:59:30  warnes
# - Modified all files to include CVS Id and Log tags.
#
#

ci  <-  function(x, confidence=0.95,alpha=1-confidence,...)
  UseMethod("ci")

ci.default <- function(x, confidence=0.95,alpha=1-confidence,na.rm=FALSE,...) {
  est <- mean(x, na.rm=na.rm)
  stderr <-  sd(x, na.rm=na.rm)/sqrt(nobs(x));
  ci.low  <- est + qt(alpha/2, nobs(x)-1) * stderr
  ci.high <- est - qt(alpha/2, nobs(x)-1) * stderr
  retval  <- c(
               Estimate=est,
               "CI lower"=ci.low,
               "CI upper"=ci.high,
               "Std. Error"=stderr,
               )

  retval
}
  
ci.lm  <-  function(x,confidence=0.95,alpha=1-confidence,...)
{
  x  <-  summary(x)
  est  <-  coef(x)[,1] ;
  ci.low  <- est + qt(alpha/2, x$df[2]) * coef(x)[,2] ;
  ci.high <- est - qt(alpha/2, x$df[2]) * coef(x)[,2] ;
  retval  <- cbind(Estimate=est,
                   "CI lower"=ci.low,
                   "CI upper"=ci.high,
                   "Std. Error"= coef(x)[,2],
                   "p-value" = coef(x)[,4])

  retval
}

ci.lme <- function(x,confidence=0.95,alpha=1-confidence,...)
  {
  x  <-  summary(x)
  est  <-  x$tTable[,"Value"] ;
  ci.low  <- est + qt(alpha/2, x$tTable[,"DF"]) * x$tTable[,"Std.Error"] ;
  ci.high <- est - qt(alpha/2, x$tTable[,"DF"]) * x$tTable[,"Std.Error"] ;
  retval  <- cbind(Estimate=est,
                   "CI lower"=ci.low,
                   "CI upper"=ci.high,
                   "Std. Error"= x$tTable[,"Std.Error"],
                   "DF" = x$tTable[,"DF"],
                   "p-value" = x$tTable[,"p-value"])
  rownames(retval)  <-  rownames(x$tTable)
  retval
}
colorpanel <- function(n,low,mid,high)
  {
    isodd <- odd(n)

    if(isodd)
      {
        n <- n+1
      }
    
    # convert to rgb
    low <- col2rgb(low)
    mid <- col2rgb(mid)
    high <- col2rgb(high)

    # determine length of each component
    lower <- floor(n/2)
    upper <- n - lower

    red  <- c(
              seq(low[1,1], mid [1,1], length=lower),
              seq(mid[1,1], high[1,1], length=upper)
              )/255

    green <- c(
               seq(low[3,1], mid [3,1], length=lower),
               seq(mid[3,1], high[3,1], length=upper)
               )/255
    
    blue <- c(
              seq(low[2,1], mid [2,1], length=lower),
              seq(mid[2,1], high[2,1], length=upper)
              )/255

    if(isodd)
      {
        red   <- red  [-(lower+1)]
        green <- green[-(lower+1)]
        blue  <- blue [-(lower+1)]
      }

    rgb(red,blue,green)
  }



# Generate red-to-green colorscale
redgreen <- function(n) colorpanel(n, 'red', 'black', 'green')
greenred <- function(n) colorpanel(n, 'green', 'black', 'red' )
bluered  <- function(n) colorpanel(n, 'blue','white','red')
redblue  <- function(n) colorpanel(n, 'red','white','blue')
# $Id: combinations.R,v 1.6 2004/05/26 13:18:45 warnes Exp $
#

##
## From email by Brian D Ripley <ripley@stats.ox.ac.uk> to r-help
## dated Tue, 14 Dec 1999 11:14:04 +0000 (GMT) in response to
## Alex Ahgarin <datamanagement@email.com>.  Original version was
## named "subsets" and was Written by Bill Venables.  
##

combinations <- function(n, r, v = 1:n, set = TRUE, repeats.allowed=FALSE)
{
  if(mode(n) != "numeric" || length(n) != 1 
     || n < 1 || (n %% 1) != 0) stop("bad value of n") 
  if(mode(r) != "numeric" || length(r) != 1 
     || r < 1 || (r %% 1) != 0) stop("bad value of r") 
  if(!is.atomic(v) || length(v) < n) 
    stop("v is either non-atomic or too short")
  if( (r > n) & repeats.allowed==FALSE)
    stop("r > n and repeats.allowed=FALSE")
  if(set) {
    v <- unique(sort(v))
    if (length(v) < n) stop("too few different elements")
  }
  v0 <- vector(mode(v), 0)
  ## Inner workhorse
  if(repeats.allowed)
    sub <- function(n, r, v)
      { 
        if(r == 0) v0 else
        if(r == 1) matrix(v, n, 1) else
        if(n == 1) matrix(v, 1, r) else
        rbind( cbind(v[1], Recall(n, r-1, v)),
              Recall(n-1, r, v[-1]))
      }
  else
    sub <- function(n, r, v)
      { 
        if(r == 0) v0 else
        if(r == 1) matrix(v, n, 1) else
        if(r == n) matrix(v, 1, n) else
        rbind(cbind(v[1], Recall(n-1, r-1, v[-1])),
              Recall(n-1, r, v[-1]))
      }
  sub(n, r, v[1:n])
}

##
## Original version by Bill Venables and cited by by Matthew
## Wiener (mcw@ln.nimh.nih.gov) in an email to R-help dated
## Tue, 14 Dec 1999 09:11:32 -0500 (EST) in response to
## Alex Ahgarin <datamanagement@email.com>
##
##


permutations <- function(n, r, v = 1:n, set = TRUE, repeats.allowed=FALSE)
{
  if(mode(n) != "numeric" || length(n) != 1 
     || n < 1 || (n %% 1) != 0) stop("bad value of n") 
  if(mode(r) != "numeric" || length(r) != 1 
     || r < 1 || (r %% 1) != 0) stop("bad value of r") 
  if(!is.atomic(v) || length(v) < n) 
    stop("v is either non-atomic or too short")
  if( (r > n) ) #& repeats.allowed==FALSE)
    stop("r > n") # and repeats.allowed=FALSE")
  if(set) {
    v <- unique(sort(v))
    if (length(v) < n) stop("too few different elements")
  }
  v0 <- vector(mode(v), 0)
  ## Inner workhorse
  if(repeats.allowed)
    sub <- function(n, r, v)
      {
        if(r==1) matrix(v,n,1) else
        if(n==1) matrix(v,1,r) else
        {
          inner  <-  Recall(n, r-1, v)
          cbind( rep( v, rep(nrow(inner),n)  ),
                 matrix( t(inner), ncol=ncol(inner), nrow=nrow(inner) * n ,
                        byrow=TRUE )
                )
        }
      }
  else
    sub <- function(n, r, v)
      {
        if(r==1) matrix(v,n,1) else
        if(n==1) matrix(v,1,r) else
        {
        X  <-  NULL
        for(i in 1:n)
          X  <-  rbind( X, cbind( v[i], Recall(n-1, r - 1, v[-i])))
        X
        }
      }

  sub(n, r, v[1:n])
}
# $Id: combine.R,v 1.2 2002/09/23 13:59:30 warnes Exp $
#
# $Log: combine.R,v $
# Revision 1.2  2002/09/23 13:59:30  warnes
# - Modified all files to include CVS Id and Log tags.
#
#

combine  <-  function(..., names=NULL)
  {
    tmp  <-  list(...)
    if(is.null(names)) names  <- names(tmp)
    if(is.null(names)) names  <- sapply( as.list(match.call()), deparse)[-1]

    if( any(
            sapply(tmp, is.matrix)
            |
            sapply(tmp, is.data.frame) ) )
      { 
        len  <-  sapply(tmp, function(x) c(dim(x),1)[1] )
        len[is.null(len)]  <-  1
        data <-  rbind( ... )
      }
    else
      {
        len  <- sapply(tmp,length)
        data  <-  unlist(tmp)
        
      }

    namelist  <- factor(rep(names, len), levels=names)
        
    return( data.frame( data, source=namelist) )
  }
# $Id: contrast.lm.R,v 1.16 2003/01/30 21:53:07 warnes Exp $
#
# $Log: contrast.lm.R,v $
# Revision 1.16  2003/01/30 21:53:07  warnes
#
# - Renamed 'contrast.lm' to 'fit.contrast'.  This new name is more
#   descriptive and makes it easier to create and use methods for other
#   classes, eg lme.
#
# - Enabled fit.contrast for lme object now that Doug Bates has provided
#   the necessary support for contrasts in the nlme package.
#
# - New contrast.lm function which generates a 'depreciated' warning and
#   calls fit.contrast
#
# - Updated help text to match changes.
#
# Revision 1.15  2003/01/02 16:04:42  warnes
#
# - Now will run on objects of class aov (but not aovlist!)
#
# Revision 1.14  2002/10/29 23:00:42  warnes
#
# - Moved make.contrasts to a separate file.
# - Enhanced make contrasts to better label contrast matrix, to give
#   how.many a default value, and to coerce vectors into row matrixes.
# - Added help page for make.contrasts.
# - Added link from contrasts.lm seealso to make.contrasts.
#
# Revision 1.13  2002/10/29 19:55:13  warnes
#
# - Fix bug that prevented contrast.lm() from working on 'aov' objects
# - Add 'aov' examples to documentation for contrast.lm
# - Added note about future support for contrast.lme.
#
# Revision 1.12  2002/09/24 15:55:16  warnes
#
# - Add ability to show confidence intervals when showall=TRUE.
#
# Revision 1.11  2002/08/01 19:37:14  warnes
#
# - Corrected documentation mismatch for ci, ci.default.
#
# - Replaced all occurences of '_' for assignment with '<-'.
#
# - Replaced all occurences of 'T' or 'F' for 'TRUE' and 'FALSE' with
#   the spelled out version.
#
# - Updaded version number and date.
#
# Revision 1.10  2002/04/15 21:28:51  warneg
# - Separated out, then commented out contrast.lme code.  The
#   contrast.lme function will become part of Bates and Pinhiero's NLME
#   library.
#
# Revision 1.9  2002/04/09 00:51:29  warneg
#
# Checkin for version 0.5.3
#
# Revision 1.8  2002/04/05 18:23:17  warneg
#
# - Updated contrast.lm to handle lme objects
# - Modified contrast.lm to compute confidence intervals even when
#   showall is true.
# - Added check and warning if conf.int is outside (0,1).  This will ensuere
#   that conf.int=TRUE does not cause nonsense results.
#
# Revision 1.7  2001/12/10 19:29:21  warneg
# incorrectly put contrast.coeff.R (now estimable.R) here.  Now correctly put contrast.factor.R back here
#
# Revision 1.1  2001/12/07 19:53:49  warneg
# Renamed 'contrast.lm.R' to 'contrast.lm.R' to highlight that this function only works on individual factors.
#
# Revision 1.5  2001/11/13 21:19:22  warneg
# - Fixed error that occured when a factor has 2 levels and only one
#   contrast is specified
#
# Revision 1.4  2001/09/18 14:14:34  warneg
#
# Fixed bug in make.contrasts.  There was leftover code expecting a
# parameter 'x' which is no lonber provided.
#
# Revision 1.3  2001/08/31 23:36:52  warneg
#
# Added make.contrasts() to allow S-Plus so that the remaining
# unspecified df for contrats are filled by contrasts orthogonal to the
# specified ones.  This results in getting the exact same value from
# computing contrasts in a one-way anova as would be obtained by
# directly computing the means and performing the approprial linalg.
#
# Revision 1.2  2001/08/31 20:46:55  warneg
# Previous version did not actually work.  This version now correctly
# computes contrasts.  It will also (in conjunction wiht RSCompat.S)
# work in S-Plus.
#


contrast.lm <- function(...)
  {
    warning("'contrast.lm' is depreciated.",
            " Please use 'fit.contrast' instead.")
    UseMethod("fit.contrast")
  }

# Posted by Ben Bolker to R-News on Fri Dec 15 2000
# http://www.r-project.org/nocvs/mail/r-help/2000/3865.html
#
# Some code (originally contributed by Ian Wilson
# <i.wilson@maths.abdn.ac.uk>

           
#  functions for the "Dirichlet function", the multidimensional
#  generalization of the beta distribution: it's the Bayesian
#  canonical # distribution for the parameter estimates of a
#  multinomial distribution.

# "pdirichlet" and "qdirichlet" (distribution function and quantiles)
# would be more difficult because you'd first have to decide how to
# define the distribution function for a multivariate distribution
# ... I'm sure this could be done but I don't know how



ddirichlet<-function(x,alpha)
## probability density for the Dirichlet function, where x=vector of
## probabilities
## and (alpha-1)=vector of observed samples of each type
## ddirichlet(c(p,1-p),c(x1,x2)) == dbeta(p,x1,x2)
{

  dirichlet1 <- function(x, alpha)
    {
      logD <- sum(lgamma(alpha)) - lgamma(sum(alpha))
      s<-sum((alpha-1)*log(x))
      exp(sum(s)-logD)
      
    }

  # make sure x is a matrix
  if(!is.matrix(x))
    if(is.data.frame(x))
      x <- as.matrix(x)
    else
      x <- t(x)

  if(!is.matrix(alpha))
    alpha <- matrix( alpha, ncol=length(alpha), nrow=nrow(x), byrow=TRUE)

  if( any(dim(x) != dim(alpha)) )
    stop("Mismatch between dimensions of 'x' and 'alpha'.")
  
  pd <- vector(length=nrow(x))
  for(i in 1:nrow(x))
    pd[i] <- dirichlet1(x[i,],alpha[i,])

  # Enforce 0 <= x[i,j] <= 1, sum(x[i,]) = 1
  pd[ apply( x, 1, function(z) any( z <0 | z > 1)) ] <- 0
  pd[ apply( x, 1, function(z) all.equal(sum( z ),1) !=TRUE) ] <- 0 
  pd
}


rdirichlet<-function(n,alpha)
## generate n random deviates from the Dirichlet function with shape
## parameters alpha
{
    l<-length(alpha);
    x<-matrix(rgamma(l*n,alpha),ncol=l,byrow=TRUE);
    sm<-x%*%rep(1,l);
    x/as.vector(sm);
}

elem <- function(object, unit=c("KB","MB","bytes"), digits=0, dimensions=FALSE)
{
  get.element.classname <- function(object, element.name)
  {
    element <- object[[match(element.name, names(object))]]
    classname <- class(element)[1]
    return(classname)
  }

  get.element.dimensions <- function(object, element.name)
  {
    element <- object[[match(element.name, names(object))]]
    if(!is.null(dim(element)))
      dimensions <- paste(dim(element), collapse=" x ")
    else
      dimensions <- length(element)
    return(dimensions)
  }

  get.element.size <- function(object, element.name)
  {
    element <- object[[match(element.name, names(object))]]
    use.zero <- c("classRepresentation", "ClassUnionRepresentation", "grob")
    if(class(element)[1] %in% use.zero)
      size <- 0
    else
      size <- object.size(element)
    return(size)
  }

  unit <- match.arg(unit)
  if(is.null(names(object)))
  {
    element.frame <- data.frame()
  }
  else
  {
    class.vector <- sapply(names(object), get.element.classname, object=object)
    size.vector <- sapply(names(object), get.element.size, object=object)
    if(is.data.frame(object))
    {
      class.vector <- c(class(row.names(object))[1], class.vector)
      size.vector <- c("<row.names>"=object.size(row.names(object)),
                       size.vector)
    }
    denominator <- switch(unit, "KB"=1024, "MB"=1024^2, 1)
    size.vector <- round(size.vector/denominator, digits)
    element.frame <- data.frame(class.vector=class.vector,
                       size.vector=size.vector, row.names=names(size.vector))
    names(element.frame) <- c("Class", unit)
    if(dimensions==TRUE && is.data.frame(object))
    {
      element.frame <- cbind(element.frame, Dim=c(nrow(object),
                         sapply(names(object),get.element.dimensions,
                         object=object)))
    }
    else if(dimensions==TRUE && !is.data.frame(object))
      element.frame <- cbind(element.frame, Dim=sapply(names(object),
                         get.element.dimensions,object=object))
  }

  return(element.frame)
}

env <- function(unit=c("KB","MB","bytes"), digits=0)
{
  get.object.size <- function(object.name, pos)
  {
    object <- get(object.name, pos=pos)
    use.zero <- c("classRepresentation", "ClassUnionRepresentation", "grob")
    if(class(object)[1] %in% use.zero)
      size <- 0
    else
      size <- object.size(object)
    return(size)
  }

  get.environment.size <- function(pos)
  {
    if(search()[pos]=="Autoloads" || length(ls(pos,all=TRUE))==0)
      size <- 0
    else
      size <- sum(sapply(ls(pos,all=TRUE), get.object.size, pos=pos))
    return(size)
  }

  get.environment.nobjects <- function(pos)
  {
    nobjects <- length(ls(pos,all=TRUE))
    return(nobjects)
  }

  unit <- match.arg(unit)
  denominator <- switch(unit, "KB"=1024, "MB"=1024^2, 1)
  size.vector <- sapply(seq(along=search()), get.environment.size)
  size.vector <- round(size.vector/denominator, digits)
  nobjects.vector <- sapply(seq(along=search()), get.environment.nobjects)
  env.frame <- data.frame(search(), nobjects.vector, size.vector)
  names(env.frame) <- c("Environment", "Objects", unit)

  print(env.frame, right=FALSE)
  invisible(env.frame)
}

estimable <- function (obj, cm, beta0, conf.int=NULL, joint.test=FALSE,
                           show.beta0) 
{
  if (!is.matrix(cm) && !is.data.frame(cm)) 
    cm <- matrix(cm, nrow=1)

  if(missing(show.beta0))
    {
      if(!missing(beta0))
        show.beta0=TRUE
      else
        show.beta0=FALSE
    }
       
  
  if (missing(beta0))
    {
      beta0 = rep(0, ifelse(is.null(nrow(cm)), 1, nrow(cm)))

    }
  
  
  if (joint.test==TRUE)
    {
      .wald(obj, cm, beta0)
    }
  else
    {
      if ("lme" %in% class(obj)) {
        stat.name <- "t.stat"
        cf <- summary(obj)$tTable
        rho <- summary(obj)$cor
        vcv <- rho * outer(cf[, 2], cf[, 2])
        tmp <- cm
        tmp[tmp==0] <- NA
        df.all <- t(abs(t(tmp) * obj$fixDF$X))
        df <- apply(df.all, 1, min, na.rm=TRUE)
        problem <- apply(df.all !=df, 1, any, na.rm=TRUE)
        if (any(problem)) 
          warning(paste("Degrees of freedom vary among parameters used to ", 
                        "construct linear contrast(s): ", 
                        paste((1:nrow(tmp))[problem], collapse=", "), 
                        ". Using the smallest df among the set of parameters.", 
                        sep=""))
      }
      else if ("lm" %in% class(obj))
        {
          stat.name <- "t.stat"
          cf <- summary.lm(obj)$coefficients
          vcv <- summary.lm(obj)$cov.unscaled * summary.lm(obj)$sigma^2
          df <- obj$df.residual
          if ("glm" %in% class(obj))
            {
              vcv <- summary(obj)$cov.scaled
              if(family(obj)[1] %in% c("poisson", "binomial"))
                {
                  stat.name <- "X2.stat"            
                  df <- 1
                }
              else
                {              
                  stat.name <- "t.stat"
                  df <- obj$df.residual
                }
            }
        }
      else if ("geese" %in% class(obj))
        {
          stat.name <- "X2.stat"
          cf <- summary(obj)$mean
          vcv <- obj$vbeta
          df <- 1
        }
      else if ("gee" %in% class(obj))
        {
          stat.name <- "X2.stat"
          cf <- summary(obj)$coef
          vcv <- obj$robust.variance
          df <- 1
        }
      else
        {
          stop("obj must be of class 'lm', 'glm', 'aov', 'lme', 'gee', 'geese' or 'nlme'")
        }
      if (is.null(cm)) 
        cm <- diag(dim(cf)[1])
      if (!dim(cm)[2]==dim(cf)[1]) 
        stop(paste("\n Dimension of ", deparse(substitute(cm)), 
                   ": ", paste(dim(cm), collapse="x"),
                   ", not compatible with no of parameters in ", 
                   deparse(substitute(obj)), ": ", dim(cf)[1], sep=""))
      ct <- cm %*% cf[, 1] 
      ct.diff <- cm %*% cf[, 1] - beta0
      
      vc <- sqrt(diag(cm %*% vcv %*% t(cm)))
      if (is.null(rownames(cm))) 
        rn <- paste("(", apply(cm, 1, paste, collapse=" "), 
                    ")", sep="")
      else rn <- rownames(cm)
      switch(stat.name, 
             t.stat={
               prob <- 2 * (1 - pt(abs(ct.diff/vc), df))
             }, 
             X2.stat={
               prob <- 1 - pchisq((ct.diff/vc)^2, df=1)
             })
      
      if (stat.name=="X2.stat")
        {
          retval <- cbind(hyp=beta0, est=ct, stderr=vc,
                          "X^2 value"=(ct.diff/vc)^2, 
                          df=df, prob=1 - pchisq((ct.diff/vc)^2, df=1))
          dimnames(retval) <- list(rn, c("beta0", "Estimate", "Std. Error", 
                                         "X^2 value", "DF", "Pr(>|X^2|)"))
        }
      else if (stat.name=="t.stat")
        {
          retval <- cbind(hyp=beta0, est=ct, stderr=vc, "t value"=ct.diff/vc, 
                          df=df, prob=2 * (1 - pt(abs(ct.diff/vc), df)))
          dimnames(retval) <- list(rn, c("beta0", "Estimate", "Std. Error", 
                                         "t value", "DF", "Pr(>|t|)"))
        }
      
      if (!is.null(conf.int))
        {
          if (conf.int <=0 || conf.int >=1) 
            stop("conf.int should be betweeon 0 and 1. Usual values are 0.95, 0.90")
          alpha <- 1 - conf.int
          switch(stat.name, t.stat={
            quant <- qt(1 - alpha/2, df)
          }, X2.stat={
            quant <- qt(1 - alpha/2, 100)
          })
          nm <- c(colnames(retval), "Lower.CI", "Upper.CI")
          retval <- cbind(retval, lower=ct.diff - vc * quant, upper=ct.diff + 
                          vc * quant)
          colnames(retval) <- nm
        }
      retval <- as.data.frame(retval)
      if(!show.beta0) retval$beta0 <- NULL
      return(retval)
    }
  }

.wald <- function (obj, cm,
                   beta0=rep(0, ifelse(is.null(nrow(cm)), 1, nrow(cm)))) 
{
    if (!is.matrix(cm) && !is.data.frame(cm)) 
        cm <- matrix(cm, nrow=1)
    df <- nrow(cm)
    if ("geese" %in% class(obj))
      {
        cf <- obj$beta
        vcv <- obj$vbeta
      }
    else if ("gee" %in% class(obj))
      {
        cf <- summary(obj)$coef
        vcv <- obj$robust.variance
      }
    else if ("lm" %in% class(obj))
      {
        cf <- summary.lm(obj)$coefficients[, 1]
        vcv <- summary.lm(obj)$cov.unscaled * summary.lm(obj)$sigma^2
        if ("glm" %in% class(obj))
          {
            vcv <- summary(obj)$cov.scaled
          }
      }
    else
      stop("obj must be of class 'lm', 'glm', 'aov', 'gee' or 'geese'")
    u <- (cm %*% cf)-beta0
    vcv.u <- cm %*% vcv %*% t(cm)
    W <- t(u) %*% solve(vcv.u) %*% u
    prob <- 1 - pchisq(W, df=df)
    retval <- as.data.frame(cbind(W, df, prob))
    names(retval) <- c("X2.stat", "DF", "Pr(>|X^2|)")
    print(as.data.frame(retval))
}
# The current implementation of the function svd() in S-Plus and R is
# much slower when operating on a matrix with a large number of
# columns than on the transpose of this matrix, which has a large
# number of rows. R's La.svd does not suffer from this problem.
#
# For R, the simple solution is to use La.svd instead of svd.  For
# S-Plus the solution is to check if the matrix is wider than tall,
# and to SVD on the transpose when this is the case.

if(exists("is.R") && is.R()==TRUE)
  {
    # This fast.prcomp() function is a slight modification of the
    # standard R prcomp function from package mva.  It uses La.svd instead
    # of the standard svd.  Consequently, it can be used with matrices
    # with many rows without a performance penalty.

    fast.prcomp <- function (x, retx = TRUE, center = TRUE, scale. = FALSE,
                        tol = NULL)
    {
        x <- as.matrix(x)
        x <- scale(x, center = center, scale = scale.)
        s <- La.svd(x, nu = 0)
        if (!is.null(tol)) {
            rank <- sum(s$d > (s$d[1] * tol))
            if (rank < ncol(x)) 
                s$vt <- s$vt[, 1:rank, drop = FALSE]
        }
        s$d <- s$d/sqrt(max(1, nrow(x) - 1))
    
        dimnames(s$vt) <- list(paste("PC", seq(len = nrow(s$vt)), sep = ""),
                               colnames(x), )
        r <- list(sdev = s$d, rotation = t(s$vt) )
        if (retx) 
            r$x <- x %*% t(s$vt)
        class(r) <- "prcomp"
        r
    }

    fast.svd <- function( x, nu = min(n, p), nv = min(n, p), ...) 
      {
        x <- as.matrix(x)
        dx <- dim(x)
        n <- dx[1]
        p <- dx[2]

        retval <- La.svd(x, nu=nu, nv=nv,  ... )
        retval$v <- t(retval$vt)
        retval$vt <- NULL
        retval
      }
    
  } else 
  {
      # The fast.svd() function checks if the number of columns is
      # larger than the number of rows.  When this is the case, it
      # transposes the matrix, calles svd, and then flips the returned
      # u and v matrixes. Otherwise it just calls svd.
      # 
      # This permits an SVD to be computed efficiently regardless of whether
      # n >>p or vice versa.

      # The fast.prcomp() function is simply a copy of the standard R
      # prcomp function from package mva which calls fast.svd instead of the
      # standard svd.  Consequently, it can be used with matrices with many
      # rows without a performance penalty.

    
    fast.prcomp <- function (x, retx = TRUE, center = TRUE, scale. = FALSE,
                        tol = NULL)
    {
        x <- as.matrix(x)
        x <- scale(x, center = center, scale = scale.)
        s <- fast.svd(x, nu = 0)
        if (!is.null(tol)) {
            rank <- sum(s$d > (s$d[1] * tol))
            if (rank < ncol(x)) 
                s$v <- s$v[, 1:rank, drop = FALSE]
        }
        s$d <- s$d/sqrt(max(1, nrow(x) - 1))
        dimnames(s$v) <- list(colnames(x), paste("PC", seq(len = ncol(s$v)), 
            sep = ""))
        r <- list(sdev = s$d, rotation = s$v)
        if (retx) 
            r$x <- x %*% s$v
        class(r) <- "prcomp"
        r
    }
    
    
    fast.svd <- function( x, nu = min(n, p), nv = min(n, p), ...) 
      {
        x <- as.matrix(x)
        dx <- dim(x)
        n <- dx[1]
        p <- dx[2]
    
        if(  p <= n )
          return( svd( x, nu, nv, ... ) )
        else
          {
          s <- svd( t(x), nu=nv, nv=nu, ...)
          retval <- list()
          retval$d <- s$d
          retval$u <- s$v
          retval$v <- s$u
          return(retval)
        }
      }

    NULL
    
   }
# $Id: fit.contrast.R,v 1.3 2003/11/17 21:40:15 warnes Exp $
#
# $Log: fit.contrast.R,v $
# Revision 1.3  2003/11/17 21:40:15  warnes
#
# - Fix incorrect handling of glm objects by fit.contrast, as reported
#   by Ulrich Halekoh, Phd <ulrich.halekoh@agrsci.dk>.
#
# - Add regression test code to for this bug.
#
# Revision 1.2  2003/04/22 17:24:29  warnes
#
# - the variable 'df' was used within the lme code section overwriting
#   the argument 'df'.
#
# Revision 1.1  2003/01/30 21:53:06  warnes
#
# - Renamed 'contrast.lm' to 'fit.contrast'.  This new name is more
#   descriptive and makes it easier to create and use methods for other
#   classes, eg lme.
#
# - Enabled fit.contrast for lme object now that Doug Bates has provided
#   the necessary support for contrasts in the nlme package.
#
# - New contrast.lm function which generates a 'depreciated' warning and
#   calls fit.contrast
#
# - Updated help text to match changes.
#
# Revision 1.15  2003/01/02 16:04:42  warnes
#
# - Now will run on objects of class aov (but not aovlist!)
#
# Revision 1.14  2002/10/29 23:00:42  warnes
#
# - Moved make.contrasts to a separate file.
# - Enhanced make contrasts to better label contrast matrix, to give
#   how.many a default value, and to coerce vectors into row matrixes.
# - Added help page for make.contrasts.
# - Added link from contrasts.lm seealso to make.contrasts.
#
# Revision 1.13  2002/10/29 19:55:13  warnes
#
# - Fix bug that prevented contrast.lm() from working on 'aov' objects
# - Add 'aov' examples to documentation for contrast.lm
# - Added note about future support for contrast.lme.
#
# Revision 1.12  2002/09/24 15:55:16  warnes
#
# - Add ability to show confidence intervals when showall=TRUE.
#
# Revision 1.11  2002/08/01 19:37:14  warnes
#
# - Corrected documentation mismatch for ci, ci.default.
#
# - Replaced all occurences of '_' for assignment with '<-'.
#
# - Replaced all occurences of 'T' or 'F' for 'TRUE' and 'FALSE' with
#   the spelled out version.
#
# - Updaded version number and date.
#
# Revision 1.10  2002/04/15 21:28:51  warneg
# - Separated out, then commented out contrast.lme code.  The
#   contrast.lme function will become part of Bates and Pinhiero's NLME
#   library.
#
# Revision 1.9  2002/04/09 00:51:29  warneg
#
# Checkin for version 0.5.3
#
# Revision 1.8  2002/04/05 18:23:17  warneg
#
# - Updated contrast.lm to handle lme objects
# - Modified contrast.lm to compute confidence intervals even when
#   showall is true.
# - Added check and warning if conf.int is outside (0,1).  This will ensuere
#   that conf.int=TRUE does not cause nonsense results.
#
# Revision 1.7  2001/12/10 19:29:21  warneg
# incorrectly put contrast.coeff.R (now estimable.R) here.  Now correctly put contrast.factor.R back here
#
# Revision 1.1  2001/12/07 19:53:49  warneg
# Renamed 'contrast.lm.R' to 'contrast.lm.R' to highlight that this function only works on individual factors.
#
# Revision 1.5  2001/11/13 21:19:22  warneg
# - Fixed error that occured when a factor has 2 levels and only one
#   contrast is specified
#
# Revision 1.4  2001/09/18 14:14:34  warneg
#
# Fixed bug in make.contrasts.  There was leftover code expecting a
# parameter 'x' which is no lonber provided.
#
# Revision 1.3  2001/08/31 23:36:52  warneg
#
# Added make.contrasts() to allow S-Plus so that the remaining
# unspecified df for contrats are filled by contrasts orthogonal to the
# specified ones.  This results in getting the exact same value from
# computing contrasts in a one-way anova as would be obtained by
# directly computing the means and performing the approprial linalg.
#
# Revision 1.2  2001/08/31 20:46:55  warneg
# Previous version did not actually work.  This version now correctly
# computes contrasts.  It will also (in conjunction wiht RSCompat.S)
# work in S-Plus.
#

fit.contrast.lm <- function(model, varname, coeff, showall=FALSE,
                            conf.int=NULL, df=FALSE, ...)
{
  # check class of model
  if( !("lm" %in% class(model) ||
        "aov" %in% class(model) ||
        "lme" %in% class(model) ) )
    stop("contrast.lm can only be applied to objects inheriting from 'lm'",
         "and 'lme' (eg: lm,glm,aov,lme).")
  
  # make sure we have the NAME of the variable
  if(!is.character(varname))
     varname <- deparse(substitute(varname))

  # make coeff into a matrix 
  if(!is.matrix(coeff))
    {
       coeff <- matrix(coeff, nrow=1)
     }

  # make sure columns are labeled
  if (is.null(rownames(coeff)))
     {
       rn <- vector(length=nrow(coeff))
       for(i in 1:nrow(coeff))
          rn[i] <- paste(" c=(",paste(coeff[i,],collapse=" "), ")")
       rownames(coeff) <- rn
     }

  # now convert into the proper form for the contrast matrix
  cmat <- make.contrasts(coeff, ncol(coeff) )
  cn <- paste(" C",1:ncol(cmat),sep="")
  cn[1:nrow(coeff)] <- rownames(coeff)
  colnames(cmat) <- cn

  # recall fitting method with the specified contrast
  m <- model$call
  if(is.null(m$contrasts))
    m$contrasts <- list()
  m$contrasts[[varname]] <- cmat 
  if(is.R())
    r <- eval(m, parent.frame())
  else
    r <- eval(m)

  # now return the correct elements ....
  if( 'lme' %in% class(model) )
    {
      est <- r$coefficients$fixed
      se  <- sqrt(diag(r$varFix))
      tval <- est/se
      df.lme   <- r$fixDF$X
      retval <- cbind(
                      "Estimate"= est,
                      "Std. Error"= se,
                      "t-value"= tval,
                      "Pr(>|t|)"=  2 * (1 - pt(abs(tval), df.lme)),
                      "DF"=df.lme
                      )

    }
  else if('glm' %in% class(model))
    {
      smodel <- summary.glm(r)
      retval <- cbind(coef(smodel), "DF"=smodel$df[2])
    }
  else # lm, aov
    {
      smodel <- summary.lm(r)
      retval <- cbind(coef(smodel), "DF"=smodel$df[2])
    }

  if( !showall )
    {
      if( !is.R() && ncol(cmat)==1 )
        {
          retval <- retval[varname,,drop=FALSE]
          rownames(retval) <- rn
        }
      else
        {
          rn <- paste(varname,rownames(coeff),sep="")
          ind <- match(rn,rownames(retval))
          retval <- retval[ind,,drop=FALSE]
        }

    }

  if(!is.null(conf.int)) # add confidence intervals
    {
      alpha <- 1-conf.int
      retval <- cbind( retval,
                      "lower CI"=retval[,1] -
                      qt(1-alpha/2,retval[,5])*retval[,2],
                      "upper CI"=retval[,1] +
                      qt(1-alpha/2,retval[,5])*retval[,2] )
    }

  if(!df)
    return(retval[,-5,drop=FALSE])
  else
    return(retval)
}

# fit.contrast.lme is necessary because 'lme' objects do not inherit
# from 'lm'.
#
# **Make sure that the argument list *exactly* matches the one
# for fit.contrast.lm() above.**
#
fit.contrast.lme <- function(model, varname, coeff, showall=FALSE,
                            conf.int=NULL, df=FALSE, ...)
  {
    require(nlme)
    fit.contrast.lm(model, varname, coeff, showall, conf.int, df)
  }


fit.contrast <- function(model, varname, coeff, ...)
  UseMethod("fit.contrast")

foldchange <- function(num,denom)
  {
    ifelse(num >= denom, num/denom, -denom/num)
  }


# Compute foldchange from log-ratio values
logratio2foldchange <- function(logratio, base=2)
  {
    retval <- base^(logratio)
    retval <- ifelse(retval < 1, -1/retval, retval)
    retval
  }

# vice versa
foldchange2logratio <- function(foldchange, base=2)
  {
    retval <- ifelse( foldchange<0, 1/-foldchange, foldchange)
    retval <- log(retval,base)
    retval
  }

# $Id: glh.test.R,v 1.8 2002/09/24 19:12:42 warnes Exp $
#
# $Log: glh.test.R,v $
# Revision 1.8  2002/09/24 19:12:42  warnes
#
# - Fixed a typo.
#
# Revision 1.7  2002/04/09 00:51:29  warneg
#
# Checkin for version 0.5.3
#
# Revision 1.6  2002/03/26 21:22:26  warneg
#
# - Changed methods to include '...' to match the generic.
# - Updated for version 0.5.1
#
# Revision 1.5  2002/01/17 23:42:39  warneg
#
# - Fixed typo in code that resulted in an syntax error.
#
# Revision 1.4  2002/01/10 17:35:41  warneg
#
# - print.glh.test() was using cat() to printing the call.  This didn't work and
# generated an error.
#
# Revision 1.3  2001/12/19 20:05:27  warneg
#
# - Removed extra element of return object.
#
# Revision 1.2  2001/12/18 21:34:25  warneg
# - Modified to work correctly when obj is of class 'aov' by specifying
#   summary.lm instead of summary.  This ensures that the summary object
#   has the fields we need.
#
# - Moved detailed reporting of results from 'print' to 'summary'
#   function and added a simpler report to 'print'
#
# Revision 1.1  2001/12/18 00:43:23  warneg
#
# Initial checkin.
#
#

glh.test <- function( reg, cm, d=rep(0, nrow(cm)) )
{

  if( !is.matrix(cm) && !is.data.frame(cm) )
    cm <- matrix(cm, nrow=1)

  if ( !( "lm" %in% class(reg) ) )
    stop("Only defined for lm,glm objects")

  bhat <- summary.lm(reg)$coefficients[,1,drop=FALSE]
  XpX <- summary.lm(reg)$cov.unscaled
  df <- reg$df.residual
  msr <- summary.lm(reg)$sigma  # == SSE / (n-p)
  r <- nrow(cm)


  if ( ncol(cm) != length(bhat) ) stop(  
                   paste( "\n Dimension of ",
                         deparse( substitute( cm ) ), ": ",
                         paste( dim(cm), collapse="x" ),
                         ", not compatible with no of parameters in ",
                         deparse( substitute( reg ) ), ": ",
                         length(bhat), sep="" ) )


  #                        -1
  #     (Cb - d)' ( C (X'X)   C' ) (Cb - d) / r
  # F = ---------------------------------------
  #                 SSE / (n-p)
  #

  Fstat <- t(cm %*% bhat - d) %*% solve((cm %*% XpX %*% t(cm))) %*% (cm %*% bhat - d) / r / msr^2 

  p <- 1-pf(Fstat,r,df)

  retval <- list()
  retval$call <- match.call()
  retval$statistic <- c(F=Fstat)
  retval$parameter <- c(df1=r,df2=df)
  retval$p.value <- p
  retval$conf.int <- NULL
  retval$estimate <- cm%*%bhat
  retval$null.value <- d
  retval$method <- "Test of General Linear Hypothesis"
  retval$data.name <- deparse(substitute(reg))
  retval$matrix <- cm
  colnames(retval$matrix) <- names(reg$coef)
  
  class(retval) <- c("glh.test","htest")

  retval
}

print.glh.test <- function(x, digits = 4, ... )
{
    cat("\n")
    cat("\t",x$method, prefix = "\t")
    cat("\n")
    cat("Call:\n")
    print(x$call)
    
    if (!is.null(x$statistic)) 
        cat(names(x$statistic), " = ", format(round(x$statistic, 
            4)), ", ", sep = "")
    if (!is.null(x$parameter)) 
        cat(paste(names(x$parameter), " = ", format(round(x$parameter, 
            3)), ",", sep = ""), "")
    cat("p-value =",
        format.pval(x$p.value, digits = digits), 
        "\n")
    cat("\n")
  }


  
summary.glh.test <- function(object, digits = 4, ... )
{
    cat("\n")
    cat("\t",object$method, prefiobject = "\t")
    cat("\n")
    cat("Regression: ", object$data.name, "\n")
    cat("\n")
    cat("Null Hypothesis: C %*% Beta-hat = d \n")
    cat("\n")
    cat("C matrix: \n")
    print(object$matrix, digits=digits)
    cat("\n")
    cat("d vector: \n")
    print(object$null.value, digits=digits)
    cat("\n")
    cat("C %*% Beta-hat: \n")
    print(c(object$estimate))
    cat("\n")
    
    if (!is.null(object$statistic)) 
        cat(names(object$statistic), " = ", format(round(object$statistic, 
            4)), ", ", sep = "")
    if (!is.null(object$parameter)) 
        cat(paste(names(object$parameter), " = ", format(round(object$parameter, 
            3)), ",", sep = ""), "")
    cat("p-value =",
        format.pval(object$p.value, digits = digits), 
        "\n")
    cat("\n")
  }


  


## original Andy Liaw; modified RG, MM, GRW

heatmap.2 <- function (x,
                     
                     # dendrogram control
                     Rowv=NULL,
                     Colv=if(symm)"Rowv" else NULL,
                     distfun = dist,
                     hclustfun = hclust,
                     dendrogram = c("both","row","column","none"),
                     symm = FALSE,
                     
                     # data scaling
                     scale = c("none","row", "column"),
                     na.rm=TRUE,
                     
                     # image plot
                     revC = identical(Colv, "Rowv"),
                     add.expr,
                     breaks,
                     col="heat.colors",
                     
                     # block sepration
                     colsep,
                     rowsep,
                     sepcolor="white",
                     
                     # cell labeling
                     cellnote,
                     notecex=1.0,
                     notecol="cyan",
                     
                     # level trace
                     trace=c("column","row","both","none"),
                     tracecol="cyan",
                     hline=median(breaks),
                     vline=median(breaks),
                     linecol=tracecol,
                     
                     # Row/Column Labeling
                     margins = c(5, 5),
                     ColSideColors,
                     RowSideColors,
                     cexRow = 0.2 + 1/log10(nr),
                     cexCol = 0.2 + 1/log10(nc),
                     labRow = NULL,
                     labCol = NULL,
                     
                     # color key + density info
                     key = TRUE,
                     density.info=c("histogram","density","none"),
                     denscol=tracecol,
                     
                     # plot labels
                     main = NULL,
                     xlab = NULL,
                     ylab = NULL,
                     
                     # extras
                     ...
                     )
{
  scale01 <- function(x, low=min(x), high=max(x) )
    {
      x <- (x-low)/(high - low)
      x
    }
 
    scale <- if(symm && missing(scale)) "none" else match.arg(scale)
    dendrogram <- match.arg(dendrogram)
    trace <- match.arg(trace)
    density.info <- match.arg(density.info)

    if(!missing(breaks) && (scale!="none"))
      warning("Using scale=\"row\" or scale=\"column\" when breaks are",
              "specified can produce unpredictable results.",
              "Please consider using only one or the other.")

    # key & density don't make sense when data is not all on the same scale
#    if(scale!="none" &&  key==TRUE)
#      {
#        warning("Key cannot be plotted when scale!=\"none\".")
#        key=FALSE
#      }


    
    if(length(di <- dim(x)) != 2 || !is.numeric(x))
      stop("`x' must be a numeric matrix")

    nr <- di[1]
    nc <- di[2]

    if(nr <= 1 || nc <= 1)
      stop("`x' must have at least 2 rows and 2 columns")

    if(!is.numeric(margins) || length(margins) != 2)
      stop("`margins' must be a numeric vector of length 2")

    ## by default order by row/col mean
    if(is.null(Rowv)) Rowv <- rowMeans(x, na.rm = na.rm)
    if(is.null(Colv)) Colv <- colMeans(x, na.rm = na.rm)

    ## get the dendrograms and reordering indices

    if( dendrogram %in% c("both","row") )
      {
        if(inherits(Rowv, "dendrogram"))
          ddr <- Rowv
        else {
          hcr <- hclustfun(distfun(x))
          ddr <- as.dendrogram(hcr)
          if(!is.logical(Rowv) || Rowv)
	    ddr <- reorder(ddr, Rowv)
        }
        rowInd <- order.dendrogram(ddr)
        if(nr != length(rowInd))
          stop("row dendrogram ordering gave index of wrong length")
      }
    else
      {
        rowInd <- order(Rowv)
      }

    if( dendrogram %in% c("both","column") )
      {
        if(inherits(Colv, "dendrogram"))
          ddc <- Colv
        else if(identical(Colv, "Rowv")) {
          if(nr != nc)
            stop('Colv = "Rowv" but nrow(x) != ncol(x)')
          ddc <- ddr
        }
        else {
          hcc <- hclustfun(distfun(if(symm)x else t(x)))
          ddc <- as.dendrogram(hcc)
          if(!is.logical(Colv) || Colv)
	    ddc <- reorder(ddc, Colv)
        }
        colInd <- order.dendrogram(ddc)
        if(nc != length(colInd))
          stop("column dendrogram ordering gave index of wrong length")
      }
    else
      {
        colInd <- order(Colv)
      }

    ## reorder x
    x <- x[rowInd, colInd]
    x.unscaled <- x

    if(is.null(labRow)) 
      labRow <- if(is.null(rownames(x))) (1:nr)[rowInd] else rownames(x)
    else
      labRow <- labRow[rowInd]
  
    if(is.null(labCol))
      labCol <- if(is.null(colnames(x))) (1:nc)[colInd] else colnames(x)
    else
      labCol <- labCol[colInd]
  
    if(scale == "row") {
	x <- sweep(x, 1, rowMeans(x, na.rm = na.rm))
	sx <- apply(x, 1, sd, na.rm = na.rm)
	x <- sweep(x, 1, sx, "/")
    }
    else if(scale == "column") {
	x <- sweep(x, 2, colMeans(x, na.rm = na.rm))
	sx <- apply(x, 2, sd, na.rm = na.rm)
	x <- sweep(x, 2, sx, "/")
    }

    ## Set up breaks and force values outside the range into the endmost bins
    if(missing(breaks) || is.null(breaks) || length(breaks)<1 )
      if(missing(col))
        breaks <- 16
      else
        breaks <- length(col)+1
    if(length(breaks)==1)
      {
        breaks <- seq( min(x, na.rm=na.rm), max(x,na.rm=na.rm), length=breaks)
      }

    nbr <- length(breaks)
    ncol <- length(breaks)-1
    
    if(class(col)=="function")
      col <- col(ncol)
    else if(is.character(col) && length(col)==1)
      col <- do.call(col,list(ncol))
    
    min.breaks <- min(breaks)
    max.breaks <- max(breaks)

    x[] <- ifelse(x<min.breaks, min.breaks, x)
    x[] <- ifelse(x>max.breaks, max.breaks, x)
    
    ## Calculate the plot layout
    lmat <- rbind(4:3, 2:1)
    lwid <- lhei <- c(1, 4)
    if(!missing(ColSideColors)) { ## add middle row to layout
	if(!is.character(ColSideColors) || length(ColSideColors) != nc)
	    stop("'ColSideColors' must be a character vector of length ncol(x)")
	lmat <- rbind(lmat[1,]+1, c(NA,1), lmat[2,]+1)
	lhei <- c(lhei[1], 0.2, lhei[2])
    }
    if(!missing(RowSideColors)) { ## add middle column to layout
	if(!is.character(RowSideColors) || length(RowSideColors) != nr)
	    stop("'RowSideColors' must be a character vector of length nrow(x)")
	lmat <- cbind(lmat[,1]+1, c(rep(NA, nrow(lmat)-1), 1), lmat[,2]+1)
	lwid <- c(lwid[1], 0.2, lwid[2])
    }
    lmat[is.na(lmat)] <- 0

    ## Graphics `output' -----------------------

    op <- par(no.readonly = TRUE)
    on.exit(par(op))
    layout(lmat, widths = lwid, heights = lhei, respect = FALSE)

    ## draw the side bars
    if(!missing(RowSideColors)) {
	par(mar = c(margins[1],0, 0,0.5))
	image(rbind(1:nr), col = RowSideColors[rowInd], axes = FALSE)
    }
    if(!missing(ColSideColors)) {
	par(mar = c(0.5,0, 0,margins[2]))
	image(cbind(1:nc), col = ColSideColors[colInd], axes = FALSE)
    }
    ## draw the main carpet
    par(mar = c(margins[1], 0, 0, margins[2]))
    if(!symm || scale != "none") x <- t(x)
    if(revC) { # x columns reversed
        iy <- nr:1
        ddr <- rev(ddr)
        x <- x[,iy]
    } else iy <- 1:nr

    image(1:nc, 1:nr, x, xlim = 0.5+ c(0, nc), ylim = 0.5+ c(0, nr),
          axes = FALSE, xlab = "", ylab = "", col=col, breaks=breaks,
          ...)
    axis(1, 1:nc, labels= labCol, las= 2, line= -0.5, tick= 0, cex.axis= cexCol)
    if(!is.null(xlab)) mtext(xlab, side = 1, line = margins[1] - 1.25)
    axis(4, iy, labels= labRow, las= 2, line= -0.5, tick= 0, cex.axis= cexRow)
    if(!is.null(ylab)) mtext(ylab, side = 4, line = margins[2] - 1.25)

    if (!missing(add.expr))
	eval(substitute(add.expr))

    # show traces
    min.scale <- min(breaks)
    max.scale <- max(breaks)
    x.scaled  <- scale01(t(x), min.scale, max.scale)

    if(trace %in% c("both","column") )
      {
        for( i in colInd )
          {
            if(!is.null(vline))
              {
                vline.vals <- scale01(vline, min.scale, max.scale)
                abline(v=i-0.5 + vline.vals, col=linecol, lty=2)
              }
            xv <- rep(i, nrow(x.scaled)) + x.scaled[,i] - 0.5
            xv <- c(xv[1], xv)
            yv <- 1:length(xv)-0.5
            lines(x=xv, y=yv, lwd=1, col=tracecol, type="s")
          }
      }

    if(trace %in% c("both","row") )
      {
        for( i in rowInd )
          {
            if(!is.null(hline))
              {
                hline.vals <- scale01(hline, min.scale, max.scale)
                abline(h=i + hline, col=linecol, lty=2)
              }
            yv <- rep(i, ncol(x.scaled)) + x.scaled[i,] - 0.5
            yv <- rev(c(yv[1], yv))
            xv <- length(yv):1-0.5
            lines(x=xv, y=yv, lwd=1, col=tracecol, type="s")
          }
      }



    
    ## add 'background' colored spaces to visually separate sections
    if(!missing(colsep))
      for(csep in colsep)
        rect(xleft=csep+0.5,   ybottom=rep(0,length(csep)),
             xright=csep+0.55, ytop=rep(ncol(x)+1,csep),
             lty=1, lwd=1, col=sepcolor, border=sepcolor)

    if(!missing(rowsep))
      for(rsep in rowsep)
        rect(xleft=0,          ybottom=nrow(x)+1-rsep-0.5,
             xright=ncol(x)+1, ytop=nrow(x)+1-rsep-0.55,
             lty=1, lwd=1, col=sepcolor, border=sepcolor)

    if(!missing(cellnote))
      text(x=c(col(cellnote)),
           y=nrow(cellnote)+1-c(row(cellnote)), labels=c(cellnote),
           col=notecol, cex=notecex)
    
    ## the two dendrograms :
    par(mar = c(margins[1], 0, 0, 0))
    if( dendrogram %in% c("both","row") )
      {
        plot(ddr, horiz = TRUE, axes = FALSE, yaxs = "i", leaflab = "none")
      }
    else
      plot.new()

    par(mar = c(0, 0, if(!is.null(main)) 5 else 0, margins[2]))
    if( dendrogram %in% c("both","column") )
      {
        plot(ddc, axes = FALSE, xaxs = "i", leaflab = "none")
      }
    else
      plot.new()

    ## title 
    if(!is.null(main)) title(main, cex.main = 1.5*op[["cex.main"]])
    
    ## Add the color-key
    if(key)
      {
        par(mar = c(5, 4, 2, 1), cex=0.75)

        min.raw <- min(x, na.rm=TRUE) # Again, modified to use scaled or unscaled (SD 12/2/03)
        max.raw <- max(x, na.rm=TRUE)
        
        z <- seq(min.raw,max.raw,length=length(col))
        image(z=matrix(z, ncol=1),
              col=col, breaks=breaks,
              xaxt="n", yaxt="n" )

        par(usr=c(0,1,0,1))
        lv <- pretty(breaks)
        xv <- scale01(as.numeric(lv), min.raw, max.raw)
        axis(1, at=xv, labels=lv)
        if(scale=="row")
          mtext(side=1,"Row Z-Score", line=2)
        else if(scale=="column")
          mtext(side=1,"Column Z-Score", line=2)
        else
          mtext(side=1,"Value", line=2)

        if(density.info=="density")
          {
            # Experimental : also plot density of data
            dens <- density(x, adjust=0.25)
            omit <- dens$x < min(breaks) | dens$x > min(breaks)
            dens$x <- dens$x[-omit]
            dens$y <- dens$y[-omit]
            dens$x <- scale01(dens$x)
            lines(dens$x, dens$y / max(dens$y) * 0.95, col=denscol, lwd=1)
            axis(2, at=pretty(dens$y)/max(dens$y) * 0.95, pretty(dens$y) )
            title("Color Key\nand Density Plot")
            par(cex=0.5)
            mtext(side=2,"Density", line=2)
          }
        else if(density.info=="histogram")
          {
            h <- hist(x, plot=FALSE, breaks=breaks)
            hx <- scale01(breaks,min.raw,max.raw)
            hy <- c(h$counts, h$counts[length(h$counts)])
            lines(hx, hy/max(hy)*0.95, lwd=1, type="s", col=denscol)
            axis(2, at=pretty(hy)/max(hy) * 0.95, pretty(hy) )
            title("Color Key\nand Histogram")
            par(cex=0.5)
            mtext(side=2,"Count", line=2)
          }
        else
          title("Color Key")
        
      }
    else
      plot.new()

    invisible(list(rowInd = rowInd, colInd = colInd))
}
# $Id: hist2d.R,v 1.5 2003/03/07 15:48:35 warnes Exp $
#
# $Log: hist2d.R,v $
# Revision 1.5  2003/03/07 15:48:35  warnes
#
# - Minor changes to code to allow the package to be provided as an
#   S-Plus chapter.
#
# Revision 1.4  2002/09/23 13:59:30  warnes
# - Modified all files to include CVS Id and Log tags.
#
#

if(is.R()) 
hist2d <- function( x,y=NULL, nbins=200, same.scale=FALSE, na.rm=TRUE, show=TRUE, col=c("black", heat.colors(12)), ... )
  {
    if(is.null(y))
      {
        if(ncol(x) != 2) stop("If y is ommitted, x must be a 2 column matirx")
        y <- x[,2]
        x <- x[,1]
      }

    if(length(nbins)==1)
      nbins <- rep(nbins,2)

    nas <- is.na(x) | is.na(y)

    if(na.rm)
      {
        x <- x[!nas]
        y <- y[!nas]
      }
    else
      stop("missinig values not permitted if na.rm=FALSE")

    if(same.scale)
      {
        x.cuts <- seq( from=min(x,y), to=max(x,y), length=nbins[1]+1)
        y.cuts <- seq( from=min(x,y), to=max(x,y), length=nbins[2]+1)
      }
    else
      {
        x.cuts <- seq( from=min(x), to=max(x), length=nbins[1]+1)
        y.cuts <- seq( from=min(y), to=max(y), length=nbins[2]+1)
      }

      
    
    index.x <- cut( x, x.cuts, include.lowest=TRUE)
    index.y <- cut( y, y.cuts, include.lowest=TRUE)

    m <- matrix( 0, nrow=nbins[1], ncol=nbins[2], 
                dimnames=list( levels(index.x),
                               levels(index.y) ) )

    for( i in 1:length(index.x) )
      m[ index.x[i], index.y[i] ] <-  m[ index.x[i], index.y[i] ] + 1

    xvals <- x.cuts[1:nbins[1]]
    yvals <- y.cuts[1:nbins[2]]

    if(show)
      image( xvals,yvals, m, col=col,...)

    invisible(list(counts=m,x=xvals,y=yvals))
  }






# $Id: interleave.R,v 1.4 2003/05/20 16:03:02 warnes Exp $
#
# $Log: interleave.R,v $
# Revision 1.4  2003/05/20 16:03:02  warnes
#
# - Omit NULL variables.
#
# Revision 1.3  2002/09/23 13:59:30  warnes
# - Modified all files to include CVS Id and Log tags.
#
# Revision 1.2  2002/04/09 00:46:04  warneg
# - Properly handle case when some or all arguments are vectors.
#
# Revision 1.1  2002/02/20 21:42:23  warneg
# Initial checkin.
#
#

interleave <- function(..., append.source=TRUE, sep=": ")
  {
    sources <- list(...)

    sources[sapply(sources, is.null)] <- NULL

    sources <- lapply(sources, function(x)
                      if(is.matrix(x) || is.data.frame(x))
                      x else t(x) )

    nrows <- sapply( sources, nrow )
    mrows <- max(nrows)
    if(any(nrows!=mrows & nrows!=1 ))
      stop("Arguments have differening numbers of rows.")

    sources <- lapply(sources, function(x)
                      if(nrow(x)==1) x[rep(1,mrows),] else x )

    tmp <- do.call("rbind",sources)

    nsources <- length(sources)
    indexes <- outer( ( 0:(nsources-1) ) * mrows , 1:mrows, "+" )

    retval <- tmp[indexes,]

    if(append.source && !is.null(names(sources) ))
      if(!is.null(row.names(tmp)) )
        row.names(retval) <- paste(format(row.names(retval)),
                                   format(names(sources)),
                                   sep=sep)
      else
        row.names(retval) <- rep(names(sources), mrows)

    retval
  }
invalid <- function(x)
  {
    if( missing(x) || is.null(x) || length(x)==0 )
      return(TRUE)
    if(is.list(x)) 
      return(all(sapply(x,invalid)))
    else if(is.vector(x))
      return(all(is.na(x)))
    else
      return(FALSE)
  }

is.what <- function(object, verbose=FALSE)
{
  do.test <- function(test, object)
  {
    result <- all(get(test)(object))
    return(result)
  }

  # get all names starting with "is."
  is.names <- unlist(sapply( search(),
                            function(name) ls(name, pattern="^is\\.")))

  # narrow to functions 
  is.functions <- is.names[sapply( is.names, function(x) is.function(get(x)) )]
  
  not.using <- c("is.element", "is.empty.model", "is.loaded", "is.mts",
                 "is.na.data.frame", "is.na.POSIXlt", "is.na<-",
                 "is.na<-.default", "is.na<-.factor", "is.pairlist", "is.qr",
                 "is.R", "is.single", "is.unsorted",
                 "is.what")
  tests <- is.functions[!(is.functions %in% not.using)]
  names(tests) <- NULL
  results <- sapply(tests, do.test, object=object)
  names(results) <- tests

  
  if(verbose == FALSE)
  {
    output <- tests[results==TRUE & !is.na(results)]
  }
  else
  {
    results[results==TRUE] <- "T"
    results[results==FALSE] <- "."
    output <- data.frame(is=I(results))
  }

  return(output)
}

keep <- function(..., list=character(0), sure=FALSE)
{
  if(missing(...) && missing(list))
    stop("Keep something, or use rm(list=ls()) to clear workspace.")
  names <- as.character(substitute(list(...)))[-1]
  list <- c(list, names)
  keep.elements <- match(list, ls(1))

  if(sure == FALSE)
    return(ls(1)[-keep.elements])
  else
    rm(list=ls(1)[-keep.elements], pos=1)
}

ll <- function(pos=1, unit=c("KB","MB","bytes"), digits=0, dimensions=FALSE,
        function.dim="", ...)
{
  get.object.classname <- function(object.name, pos)
  {
    object <- get(object.name, pos=pos)
    classname <- class(object)[1]
    return(classname)
  }

  get.object.dimensions <- function(object.name, pos)
  {
    object <- get(object.name, pos=pos)
    if(class(object) == "function")
      dimensions <- function.dim
    else if(!is.null(dim(object)))
      dimensions <- paste(dim(object), collapse=" x ")
    else
      dimensions <- length(object)
    return(dimensions)
  }

  get.object.size <- function(object.name, pos)
  {
    object <- get(object.name, pos=pos)
    use.zero <- c("classRepresentation", "ClassUnionRepresentation", "grob")
    if(class(object)[1] %in% use.zero)
      size <- 0
    else
      size <- object.size(object)
    return(size)
  }

  unit <- match.arg(unit)
  denominator <- switch(unit, "KB"=1024, "MB"=1024^2, 1)

  if(is.character(pos))
    pos <- match(pos, search())

  if(length(ls(pos,...)) == 0)
  {
    object.frame <- data.frame()
  }
  else if(search()[pos] == "Autoloads")
  {
    object.frame <- data.frame(rep("function",length(ls(pos,...))),
                      rep(0,length(ls(pos,...))), row.names=ls(pos,...))
    if(dimensions == TRUE)
    {
      object.frame <- cbind(object.frame, rep("",nrow(object.frame)))
      names(object.frame) <- c("Class", unit, "Dimensions")
    }
    else
      names(object.frame) <- c("Class", unit)
  }
  else
  {
    class.vector <- sapply(ls(pos,...), get.object.classname, pos=pos)
    size.vector <- sapply(ls(pos,...), get.object.size, pos=pos)
    size.vector <- round(size.vector/denominator, digits)
    object.frame <- data.frame(class.vector=class.vector,
                      size.vector=size.vector, row.names=names(size.vector))
    names(object.frame) <- c("Class", unit)
    if(dimensions == TRUE)
      object.frame <- cbind(object.frame, Dim=sapply(ls(pos,...),
                        get.object.dimensions, pos=pos))
  }

  return(object.frame)
}

# $Id: logit.R,v 1.1 2003/05/23 18:14:12 warnes Exp $

logit <- function(x, min=0, max=1)
  {
    p <- (x-min)/(max-min)
    log(p/(1-p))
  }

inv.logit <- function(x, min=0, max=1)
  {
    p <- exp(x)/(1+exp(x))
    p <- ifelse( is.na(p) & !is.na(x), 1, p ) # fix problems with +Inf
    p * (max-min) + min
  }
                 
# $Id: lowess.R,v 1.9 2004/01/21 05:17:07 warnes Exp $
#
# $Log: lowess.R,v $
# Revision 1.9  2004/01/21 05:17:07  warnes
# Track R 1.9.0's move of 'lowess' from the base package to the (new)
# stats package.
#
# Revision 1.8  2003/03/07 15:43:44  warnes
# - Add 'NULL' as the last element of if statement that defines
#   lowess.default so that when the file is sourced, S-Plus doesn't
#   display the function definition.
#
# Revision 1.7  2003/03/07 15:41:44  warnes
#
# - Specify where the defualt lowess function should be found.
# - Use getFunction in S-Plus instead of 'get'
#
# Revision 1.6  2003/01/02 16:07:35  warnes
#
# - Added wrapper code so that R-specific fiddling won't be executed
#   under S-Plus.
#
# Revision 1.5  2002/09/23 13:59:30  warnes
# - Modified all files to include CVS Id and Log tags.
#
#

if(is.R())
  {
    # make original lowess into the default method
    if(R.version$major == 1 && R.version$minor < 9)
      lowess.default <- base::lowess
    else
      lowess.default <- stats::lowess

    lowess  <- function(x,...)
      UseMethod("lowess")

    # add "..." to the argument list to match the generic
    formals(lowess.default) <- c(formals(lowess.default),alist(...= ))

    NULL

  } else
  {

    # make original lowess into the default method
    lowess.default  <- getFunction("lowess",where="main")

    lowess  <- function(x,...)
      UseMethod("lowess")

    NULL
  }



"lowess.formula" <-  function (formula,
                               data = parent.frame(), subset, na.action, 
                               f=2/3,  iter=3,
                               delta=.01*diff(range(mf[-response])), ... )
{
  if (missing(formula) || (length(formula) != 3)) 
    stop("formula missing or incorrect")
  if (missing(na.action)) 
    na.action <- getOption("na.action")
  m <- match.call(expand.dots = FALSE)
  if (is.matrix(eval(m$data, parent.frame()))) 
    m$data <- as.data.frame(data)
  m$...  <- m$f <- m$iter <- m$delta <- NULL
  m[[1]] <- as.name("model.frame")
  mf <- eval(m, parent.frame())
  response <- attr(attr(mf, "terms"), "response")
  lowess.default(mf[[-response]], mf[[response]], f=f, iter=iter, delta=delta)
}
# $Id: make.contrasts.R,v 1.2 2003/01/30 14:58:31 warnes Exp $
#
# $Log: make.contrasts.R,v $
# Revision 1.2  2003/01/30 14:58:31  warnes
#
# - Added explicit check to ensure that the number of specified
#   contrasts is less than or equal to the ncol - 1.  Previously, this
#   failed with an obtuse error message when the contrast matrix had row
#   names, and silently dropped contrasts over ncol-1.
#
# Revision 1.1  2002/10/29 23:00:43  warnes
#
# - Moved make.contrasts to a separate file.
# - Enhanced make contrasts to better label contrast matrix, to give
#   how.many a default value, and to coerce vectors into row matrixes.
# - Added help page for make.contrasts.
# - Added link from contrasts.lm seealso to make.contrasts.
#
#


"make.contrasts" <-  function (contr, how.many=ncol(contr)) 
{
  if(!is.matrix(contr))
    contr <- matrix(contr,ncol=length(contr))

  if(nrow(contr)+1 > how.many)
    stop("Too many contrasts specified. Must be less than the number of factor levels (columns).")
  
  value <- as.matrix(ginv(contr))  # requires library(MASS)
  if (nrow(value) != how.many) 
    stop("wrong number of contrast matrix rows")
  n1 <- if (missing(how.many)) 
    how.many - 1
  else how.many
  nc <- ncol(value)
  if (nc < n1) {
    cm <- qr(cbind(1, value))
    if (cm$rank != nc + 1) 
      stop("singular contrast matrix")
    cm <- qr.qy(cm, diag(how.many))[, 2:how.many, drop=FALSE]
    cm[, 1:nc] <- value
  }
  else cm <- value[, 1:n1, drop = FALSE]

  colnames(cm) <- paste( "C", 1:ncol(cm), sep="")
  rownames(cm) <- paste( "V", 1:nrow(cm), sep="")
  
  if(!is.null(rownames(contr)))
    {
      namelist <- rownames(contr)
      colnames(cm)[1:length(namelist)] <- namelist
    }

  if(!is.null(colnames(contr)))
    rownames(cm) <- colnames(contr)

  cm
}
# select the columns which match/don't match a set of include/omit patterns.

matchcols <- function(object, with, without, method=c("and","or"), ...)
  {
    method <- match.arg(method)
    cols <- colnames(object)
    
    # include columns matching 'with' pattern(s)
    if(method=="and")
      for(i in 1:length(with))
        {
          if(length(cols)>0)
            cols <- grep(with[i], cols, value=TRUE, ...)
        }
    else
      if(!missing(with))
        if(length(cols)>0)
          cols <- sapply( with, grep, x=cols, value=TRUE, ...)

    # exclude columns matching 'without' pattern(s)
    if(!missing(without))
      for(i in 1:length(without))
        if(length(cols)>0)
          {
            omit <- grep(without[i], cols, ...)
            if(length(omit)>0)
              cols <- cols[-omit]
          }

    cols
  }
mixedorder <- function(x)
  {
    # - Split each each character string into an vector of strings and
    #   numbers 
    # - Separately rank numbers and strings
    # - Combine orders so that strings follow numbers


    numeric <- function(x)
      {
        optwarn = options("warn")
        on.exit( options(optwarn) )
        options(warn=-1)
        as.numeric(x)
      }

    nonnumeric <- function(x)
      {
        optwarn = options("warn")
        on.exit( options(optwarn) )
        options(warn=-1)

        ifelse(is.na(as.numeric(x)), x, NA)
      }

    
    x <- as.character(x)

    which.nas <- which(is.na(x))
    which.blanks <- which(x=="")

    x[ which.blanks ] <- -Inf
    x[ which.nas ] <- "\377"

    ####
    # - Convert each character string into an vector containing single
    #   character and  numeric values.
    ####
      
    # find and mark numbers in the form of +1.23e+45.67 
    delimited <- gsub("([+-]{0,1}[0-9\.]+([eE][\+\-]{0,1}[0-9\.]+){0,1})",
                      "$\\1$", x)

    # separate out numbers
    step1 <- strsplit(delimited, "\\$")

    # remove empty elements
    step1 <- sapply( step1, function(x) x[x>""] )

    # create numeric version of data
    step1.numeric <- sapply( step1, numeric )

    # create non-numeric version of data
    step1.character <- sapply( step1, nonnumeric )

    # now transpose so that 1st vector contains 1st element from each
    # original string
    maxelem <- max(sapply(step1, length))
    
    step1.numeric.t <- lapply(1:maxelem,
                              function(i)
                                 sapply(step1.numeric,
                                        function(x)x[i])
                              )
                            
    step1.character.t <- lapply(1:maxelem,
                              function(i)
                                 sapply(step1.character,
                                        function(x)x[i])
                              )
    
    # now order them
    rank.numeric   <- sapply(step1.numeric.t,rank)
    rank.character <- sapply(step1.character.t,
                             function(x) as.numeric(factor(x)))

    # and merge
    rank.numeric[!is.na(rank.character)] <- 0  # mask off string values
    rank.character <- t(t(rank.character) +
                        apply(rank.numeric,2,max,na.rm=TRUE))
    rank.overall <- ifelse(is.na(rank.character),rank.numeric,rank.character)

    order <- do.call("order",as.data.frame(rank.overall))

    return(order)
  }

mixedsort <- function(x) x[mixedorder(x)]


# $Id: nobs.R,v 1.4 2002/09/23 13:59:30 warnes Exp $
#
# $Log: nobs.R,v $
# Revision 1.4  2002/09/23 13:59:30  warnes
# - Modified all files to include CVS Id and Log tags.
#
# Revision 1.3  2002/03/26 19:29:31  warneg
#
# Added ... to methods.
#
# Revision 1.2  2002/03/26 14:28:02  warneg
# - Added CVS tags
#
#

nobs <- function(x,...)
  UseMethod("nobs",x)

nobs.default <- function(x, ...) sum( !is.na(x) )

nobs.data.frame <- function(x, ...)
  sapply(x, nobs.default)

nobs.lm <- function(x, ...)
  nobs.default(x$residuals)
# detect odd/even integers
odd <- function(x) x!=as.integer(x/2)*2
even <- function(x) x==as.integer(x/2)*2 
ooplot <- function(data, ...) UseMethod("ooplot")

ooplot.default <- function(data, width=1, space=NULL, names.arg=NULL, 
                           legend.text=NULL, horiz=FALSE, 
                           density=NULL, angle=45, kmg="fpnumkMGTP", 
                           kmglim=TRUE, 
                           type=c("xyplot", "linear", "barplot", "stackbar"), 
                           col=heat.colors(NC), prcol=NULL, 
                           border=par("fg"), main=NULL, sub=NULL, 
                           xlab=NULL, ylab=NULL, xlim=NULL, ylim=NULL, 
                           xpd=TRUE, log="", axes=TRUE, 
                           axisnames=TRUE, prval=TRUE, lm=FALSE,
                           cex.axis=par("cex.axis"), 
                           cex.names=par("cex.axis"),
                           cex.values=par("cex"),inside=TRUE, 
                           plot=TRUE, axis.lty=0, plot.ci=FALSE, 
                           ci.l=NULL, ci.u=NULL, ci.color="black", 
                           ci.lty="solid", ci.lwd=1, plot.grid=FALSE, 
                           grid.inc=NULL, grid.lty="dotted", 
                           grid.lwd=1, grid.col="black", add=FALSE, 
                           by.row=FALSE, ...)
{

  ##
  ##  oopplot function block
  ##
  ## this is the location of the helper functions for this method
  ##
  
  
  optlim <- function(lim, log=FALSE) {
    ## define xlim and ylim, adjusting for log-scale if needed
    factor <- 1.05
    if (log) {
      min <- 10^floor(log10(lim[1]/factor))
      max <- 10^ceiling(log10(factor*lim[2]))
    } else {
      range <- (factor*lim[2]-lim[1]/factor)
      if (range>0) {
        ## we know the range, now find the optimal start and endpoints
        scale <- 10^floor(log10(range))
        min <- scale*(floor((lim[1]/factor)/scale))
        max <- scale*(ceiling(factor*lim[2]/scale))
      } else {
        min=0
        max=1
      }
    }
    if (type=="barplot" || type=="stackbar") max=max+1
    
    return(c(min, max))
  }

  linearfitplot <- function(x,y,xlim,col) {
    ## calculate a linear fit through the datapoints and plot
    local <- data.frame(x=x,y=y)
    local.lm <- lm(y ~ x, data=local)
    summary <- summary(local.lm)
    xmin=min(x)
    xmax=max(x)
    if (xlim[1]<xmin) {
      x1=mean(xlim[1],xmin)
    } else {
      x1=max(min(x),xlim[1])
    }
    if (xlim[2]>xmax) {
      x2=mean(xlim[2],xmax)
    } else {
      x2=min(max(x),xlim[2])
    }
    y1=summary$coefficients[2,1]*x1+summary$coefficients[1,1]  
    y2=summary$coefficients[2,1]*x2+summary$coefficients[1,1]
    p1=c(x1,x2)
    p2=c(y1,y2)
    lines(p1,p2,col=col,lty=2,lwd=2)   
  }


  castNA <- function(matrix) {
    newmatrix <- matrix
    for (j in 1:ncol(matrix)) {
      for (i in 1:nrow(matrix)) {
        newmatrix[i,j] <- ifelse(is.na(matrix[i,j]),0,matrix[i,j])
      }
    }
    return(newmatrix)
  }
  ##
  ## End of function block
  ## 
  ## In R, most people think about data in columns rather than rows
  if(by.row)
    data <- as.matrix(data)
  else
    data <- t(as.matrix(data))

  ## make sure we only accept the supported plot options
  type <- match.arg(type)
  ## check data validity
  if( nrow(data) < 2 ) stop("At least 2 columns are required.") 
  ## check defaults 
  if (!missing(inside)) .NotYetUsed("inside", error=FALSE)# -> help(.)
  ## set the beside parameter
  if (type=="stackbar") {
    beside <- FALSE
    data <- castNA(data)
  }
  else 
    beside <- TRUE
  
  ## split the data into x and y values
  height <- data[-1,,drop=FALSE]

  
  heightscale <- ""
  heightsymbol <- ""
  ##
  if ((kmg!="") && (kmglim==TRUE)){   
    ##
    ## auto scale the parameters 
    ##
    ## now scale the data, valid factors are
    ## P : peta=1E15
    ## T : tera=1E12
    ## G : giga=1E09
    ## M : mega=1E06
    ## k : kilo=1E03
    ## m : milli=1E-03
    ## u : micro=1E-06
    ## n : nano=1E-09
    ## p : pico=1E-12
    ## f : femto=1E-15
    maxheight <- max(abs(height), na.rm=TRUE)
    heightfactor <- 1
    if (maxheight>10E15 && any(grep("P", kmg))) {
      heightfactor <- 1E15
      heightscale <- "Peta"
      heightsymbol <- "P"
    }
    else if (maxheight>1E12 && any(grep("T", kmg))) {
      heightfactor <- 1E12
      heightscale <- "Tera"
      heightsymbol <- "T"
    }
    else if (maxheight>1E09 && any(grep("G", kmg))) {
      heightfactor <- 1E09
      heightscale <- "Giga"
      heightsymbol <- "G"
    }
    else if (maxheight>1E06 && any(grep("M", kmg))) {
      heightfactor <- 1E06
      heightscale <- "Mega"
      heightsymbol <- "M"
    }
    else if (maxheight>1E03 && any(grep("k", kmg))) {
      heightfactor <- 1E03
      heightscale <- "Kilo"
      heightsymbol <- "k"
    }
    else if (maxheight<1E-15 && any(grep("f", kmg))) {
      heightfactor <- 1E-15
      heightscale <- "Femto"
      heightsymbol <- "f"
    }
    else if (maxheight<1E-12 && any(grep("p", kmg))) {
      heightfactor <- 1E-12
      heightscale <- "Pico"
      heightsymbol <- "p"
    }
    else if (maxheight<1E-09 && any(grep("n", kmg))) {
      heightfactor <- 1E-09
      heightscale <- "Nano"
      heightsymbol <- "n"
    }
    else if (maxheight<1E-06 && any(grep("u", kmg))) {
      heightfactor <- 1E-06
      heightscale <- "Micro"
      heightsymbol <- "u"
    }
    else if (maxheight<1E-03 && any(grep("m", kmg))) {
      heightfactor <- 1E-03
      heightscale <- "Milli"
      heightsymbol <- "m"
    }

    height <- height/heightfactor
    
  }

  ## fill the xaxis data set and matching rownames
  xaxis <- data[1, ]
  rownames <- rownames(data)

  ##
  if (missing(space))
    space <- if (is.matrix(height) && beside) c(0, 1) else 0.2  
  space <- space * mean(width)
  ##
  if (plot && axisnames && missing(names.arg)) {
    if (type=="xyplot")
      names.arg <- colnames(height)
    else if (type=="linear")
      names.arg <- xaxis
    else if (type=="barplot")
      names.arg <- xaxis
    else if (type=="stackbar")
      names.arg <- xaxis
    else names.arg <- xaxis
  }


  ## set the legend text if it is null
  if(is.logical(legend.text))
    legend.text <- rownames[-c(1)]

  ##  if(legend.text && is.matrix(height)) rownames(height) 
  ##  else colnames(height)
  
  if (is.vector(height) || is.array(height))
    {
      height <- rbind(height)
    }
  else if (!is.matrix(height))
    stop("`height' must be a vector or a matrix")

  ## Check for log scales
  logx <- FALSE
  logy <- FALSE
  
  if (log !="")
    {
      if (any(grep("x", log)))
        logx <- TRUE
      if (any(grep("y", log)))
        logy <- TRUE
    }
  
  ## Cannot "hatch" with rect() when log scales used
  if ((logx || logy) && !is.null(density))
    stop("Cannot use shading lines in bars when log scale is used")

  ## Set the size of Rows and Columns
  NR <- nrow(height)
  NC <- ncol(height)
  ## w.r, w.l, w.m are the x-axis coordinates

  if (type=="barplot") {
    if (NR<1) NR <- 1
    if (NC<1) NC <- 1
    if (length(space)==2) {
      space <- rep.int(c(space[2], rep.int(space[1], NR - 1)), NC)
    }
    width <- rep(width, length=NR * NC)
  } else if (type=="stackbar") {
    width <- rep(width, length=NC)
  }  else {
    width <- rep(width, length=NC)  
    delta <- width / 2
  }

  ## set the proper x-axis scale
  ## linear is a switch between equidistant and scaled
  if (type=="xyplot") {
    delta <- 0
    w.m <- xaxis
    w.r <- w.m + delta
    w.l <- w.m - delta
  } 
  else if (type=="linear") {
    w.r <- cumsum(width)
    w.m <- w.r - delta
    w.l <- w.m - delta
    xaxis <- w.m
  }
  else if(type=="barplot" || type=="stackbar") {
    delta <- width / 2
    w.r <- cumsum(space + width)
    w.m <- w.r - delta
    w.l <- w.m - delta
    ##if graphic will be stacked bars, do not plot ci
    if (beside && (NR > 1) && plot.ci)
      plot.ci=FALSE
  }
  else stop("Unkown plot type")


  ## error check ci arguments
  if (plot && plot.ci)
    {
      if ((missing(ci.l)) || (missing(ci.u)))
        stop("confidence interval values are missing")
      
      if (is.vector(ci.l))
        ci.l <- cbind(ci.l)
      else if (is.array(ci.l) && (length(dim(ci.l))==1))
        ci.l <- rbind(ci.l)
      else if (!is.matrix(ci.l))
        stop("`ci.l' must be a vector or a matrix")
      
      if (is.vector(ci.u))
        ci.u <- cbind(ci.u)
      else if (is.array(ci.u) && (length(dim(ci.u))==1))
        ci.u <- rbind(ci.u)
      else if (!is.matrix(ci.u))
        stop("`ci.u' must be a vector or a matrix")
      
      if ( any(dim(height) !=dim(ci.u) ) )
        stop("'height' and 'ci.u' must have the same dimensions.")
      else if ( any( dim(height) !=dim(ci.l) ) )
        stop("'height' and 'ci.l' must have the same dimensions.")
    }
  
  if (beside)
    w.m <- matrix(w.m, nc=NC)

  ## check height/ci.l if using log scale to prevent log(<=0) error
  ## adjust appropriate ranges and bar base values
  if ((logx && horiz) || (logy && !horiz))
    {
      if (min(height, na.rm=TRUE) <=0)
        stop("log scale error: at least one 'height' value <=0")
      
      if (plot.ci && (min(ci.l) <=0))
        stop("log scale error: at least one lower c.i. value <=0")
      
      if (logx && !is.null(xlim) && (xlim[1] <=0))
        stop("log scale error: 'xlim[1]' <=0")
      
      if (logy && !is.null(ylim) && (ylim[1] <=0))
        stop("'log scale error: 'ylim[1]' <=0")
      
      ## arbitrary adjustment to display some of bar for min(height)
      ## or min(ci.l)
      if (plot.ci)
        rectbase <- rangeadj <- (0.9 * min(c(height, ci.l)))
      else
        rectbase <- rangeadj <- (0.9 * min(height))
      
      ## if axis limit is set to < above, adjust bar base value
      ## to draw a full bar
      if (logy && !is.null(ylim))
        rectbase <- ylim[1]
      else if (logx && !is.null(xlim))
        rectbase <- xlim[1]
      
      ## if stacked bar, set up base/cumsum levels, adjusting for log scale
      if (type=="stackbar") {
        heightdata <- height ## remember the original values
        height <- rbind(rectbase, apply(height, 2, cumsum))
      }
      
      ## if plot.ci, be sure that appropriate axis limits are set to
      ## include range(ci)
      lim <-
        if (plot.ci)
          c(height, ci.l, ci.u)
        else
          height
    }
  else
    {
      ## Use original bar base value
      rectbase <- 0
      
      ## if stacked bar, set up base/cumsum levels
      if (type=="stackbar") {
        heightdata <- height ## remember the original values

        height <- rbind(rectbase, apply(height, 2, cumsum))
      }

      ## if plot.ci, be sure that appropriate axis limits are set to
      ## include range(ci)
      lim <-
        if (plot.ci)
          c(height, ci.l, ci.u)
        else
          height
      
      ## use original range adjustment factor
      rangeadj <- (-0.01 * lim)
    }
  

  ## calculate the ranges ourselves
  if (missing(xlim)) xlim <- optlim(range(w.l, w.r, na.rm=TRUE), logx)
  if (missing(ylim)) ylim <- optlim(range(height, na.rm=TRUE), logy)
  ##
  if(plot) ##-------- Plotting :
    {
      if (type=="barplot" || type=="stackbar") {
        if (horiz)

          rectbase <- xlim[1] ## make sure the xlimit and the rectbase
                              ## are in sync
        else
          rectbase <- ylim[1] ## make sure the ylimit and rectbase are in sync
        opar <-
          if (horiz)
            par(xaxs="i", xpd=xpd)
          else
            par(yaxs="i", xpd=xpd)
      } else {
        opar <-
          if (horiz)
            par(xaxs="r", yaxs="r", xpd=xpd)
          else
            par(xaxs="r", yaxs="r", xpd=xpd)
      }
      on.exit(par(opar))
      
      ## If add=FALSE open new plot window
      ## else allow for adding new plot to existing window
      if (!add)
        {
          plot.new()
          plot.window(xlim, ylim, log=log, ...)
        }
      
      ## Set plot region coordinates
      usr <- par("usr")
      
      ## adjust par("usr") values if log scale(s) used
      if (logx)
        {
          usr[1] <- 10 ^ usr[1]
          usr[2] <- 10 ^ usr[2]
        }
      
      if (logy)
        {
          usr[3] <- 10 ^ usr[3]
          usr[4] <- 10 ^ usr[4]
        }
      
      ## if prcol specified, set plot region color
      if (!missing(prcol))
        rect(usr[1], usr[3], usr[2], usr[4], col=prcol)
      
      ## if plot.grid, draw major y-axis lines if vertical or x axis
      ## if horizontal R V1.6.0 provided axTicks() as an R equivalent
      ## of the C code for CreateAtVector.  Use this to determine
      ## default axis tick marks when log scale used to be consistent
      ## when no grid is plotted.  Otherwise if grid.inc is specified, 
      ## use pretty()
      
      if (plot.grid)
      {
        par(xpd=FALSE)

        if (is.null(grid.inc))
        {
          if (horiz)
          {
            grid <- axTicks(1)
            abline(v=grid, lty=grid.lty, lwd=grid.lwd, col=grid.col)
          }
          else
          {
            grid <- axTicks(2)
            abline(h=grid, lty=grid.lty, lwd=grid.lwd, col=grid.col)
          }
        }
        else
        {
          if (horiz)
          {
            grid <- pretty(xlim, n=grid.inc)
            abline(v=grid, lty=grid.lty, lwd=grid.lwd, col=grid.col)
          }
          else
          {
            grid <- pretty(ylim, n=grid.inc)
            abline(h=grid, lty=grid.lty, lwd=grid.lwd, col=grid.col)
          }
        }

         par(xpd=xpd)
      }

      ##
      ## end of the general setup, now get ready to plot the elements
      ##  
      ## cycle through all the sets and plot the lines
      if (type=="xyplot" || type=="linear") {
        pch <- c()
        for (i in 1:NR) {
          list <- is.finite(height[i, ])
          xrange <- xaxis[list]
          yrange <- height[i, list]
          lines(xrange, yrange, col=col[i])
          if ((type=="xyplot") && (lm==TRUE)) linearfitplot(xrange,yrange,xlim,col[i])
          symbol=21+(i %% 5)
          points(xrange, yrange, pch=symbol, bg=col[i], col=col[i])
          pch <- c(pch, symbol)
          ##
          if (prval) {
            values=paste(format(yrange, digits=4), heightsymbol, sep="")
            chh <- par()$cxy[2]
            chw <- par()$cxy[1]
            if (logy) {
              factor=(yrange+chh)/yrange
              text(xrange, 1.1*yrange, labels=values, adj=c(0, 1), 
	           srt=0, cex=cex.values, bg="white", col="navy") 
              } else {
              text(xrange, yrange+chh, labels=values, adj=c(0, 1), 
	           srt=0, cex=cex.values, bg="white", col="navy")
            }
          }
        }
      }
      
      if (type=="barplot" || type=="stackbar") {
        xyrect <- function(x1, y1, x2, y2, horizontal=TRUE, ...)
          {
            if(horizontal)
              rect(x1, y1, x2, y2, ...)
            else
              rect(y1, x1, y2, x2, ...)
          }

        chh <- par()$cxy[2]
        chw <- par()$cxy[1]
        if (type=="barplot") {
          xyrect(rectbase, w.l, c(height), w.r, horizontal=horiz, 
                 angle=angle, density=density, col=col, border=border)
          if (prval) {
            values=paste(format(c(height), digits=4), heightsymbol, sep="")
                if (horiz) {
                   text(c(height)+chw, w.l, labels=values, adj=c(0, 1),
                        srt=0, cex=cex.values, bg="white", col="navy")
                } else {
                  text(w.l, c(height)+chh, labels=values, adj=c(0, 1),
                       srt=0, cex=cex.values, bg="white", col="navy")
         }
          }
        } else if (type=="stackbar") {
          for (i in 1:NC) {
            xyrect(height[1:NR, i], w.l[i], height[-1, i], w.r[i], 
                   horizontal=horiz, angle=angle, density=density, 
                   col=col, border=border)
            if (prval) {
              for (j in 1:NR) {
                values=paste(format(c(heightdata[j, i]), digits=4), 
                  heightsymbol, sep="")
                if (horiz) {
                   text(c(height[j, i])+chw, w.m[i], labels=values, 
                        adj=c(0, 1), srt=0, cex=cex.values, bg="white",
                        col="navy")
                } else {
                text(w.m[i], c(height[j, i])+chh, labels=values, adj=c(0, 1), 
                     srt=0, cex=cex.values, bg="white", col="navy")
              }
              }
            }
          }
        }
      }
      
      if (plot.ci)
        {
          ## CI plot width=barwidth / 2
          ci.width=width / 4
          
          if (horiz)
            {
              segments(ci.l, w.m, ci.u, w.m, col=ci.color, lty=ci.lty, 
                       lwd=ci.lwd)
              segments(ci.l, w.m - ci.width, ci.l, w.m + ci.width, 
                       col=ci.color, lty=ci.lty, lwd=ci.lwd)
              segments(ci.u, w.m - ci.width, ci.u, w.m + ci.width, 
                       col=ci.color, lty=ci.lty, lwd=ci.lwd)
            }
          else
            {
              segments(w.m, ci.l, w.m, ci.u, col=ci.color, lty=ci.lty, 
                       lwd=ci.lwd)
              segments(w.m - ci.width, ci.l, w.m + ci.width, ci.l, 
                       col=ci.color, lty=ci.lty, lwd=ci.lwd)
              segments(w.m - ci.width, ci.u, w.m + ci.width, ci.u, 
                       col=ci.color, lty=ci.lty, lwd=ci.lwd)
            }
        }
      
      
      if (axisnames && !is.null(names.arg)) # specified or from {col}names
        {
          at.l <-
            if (length(names.arg) !=length(w.m))
              {
                if (length(names.arg)==NC) # i.e. beside (!)
                  colMeans(w.m)
                else
                  if ((type=="barplot") && (NR==1)) {
                    median(w.m)
                  } else if ((type=="linear") && (NR==1)) {
                    median(w.m) 
                  } else {
                    stop("incorrect number of names now")
                  }
              }
            else w.m

          ##axis(if(horiz) 2 else 1, at=at.l, labels=names.arg, 
          ##     lty=axis.lty, cex.axis=cex.names, ...)
        }
      
      
      if(!is.null(legend.text))
        {
          legend.col <- rep(col, length=length(legend.text))
          
          if((horiz & beside) || (!horiz & !beside))
            {
              legend.text <- rev(legend.text)
              legend.col <- rev(legend.col)
              density <- rev(density)
              angle <- rev(angle)
            }
          
          ## adjust legend x and y values if log scaling in use
          if (logx)
            legx <- usr[2] - ((usr[2] - usr[1]) / 10)
          else
            legx <- usr[2] - xinch(0.1)
          
          if (logy)
            legy <- usr[4] - ((usr[4] - usr[3]) / 10)
          else
            legy <- usr[4] - yinch(0.1)

          if (type=="barplot" || type=="stackbar") {
            legend(legx, legy, 
                   legend=legend.text, angle=angle, density=density, 
                   fill=legend.col, xjust=1, yjust=1)
          } else {
            legend(legx, legy, 
                   legend=legend.text, 
                   ##angle=angle, 
                   ##density=density, 
                   lwd=1, 
                   col=legend.col, 
                   pch=pch, 
                   pt.bg=legend.col, 
                   xjust=1, 
                   yjust=1)
          }
        }
      
      title(main=main, sub=sub, xlab=xlab, ylab=paste(heightscale, ylab), ...)
      
      ## if axis is to be plotted, adjust for grid "at" values
      if (axes)
        {
          par(lab=c(10, 10, 7))
          if (type=="barplot" || type=="stackbar") {
            axis(if(horiz) 2 else 1, at=at.l, labels=names.arg, lty=axis.lty, 
                 cex.axis=cex.names, ...)
            if(plot.grid)
              axis(if(horiz) 1 else 2, at=grid, cex.axis=cex.axis, ...)
            else
              axis(if(horiz) 1 else 2, cex.axis=cex.axis, ...)
          }
          else if (type=="linear") {
            if(plot.grid) {
              axis(if(horiz) 1 else 2, at=grid, cex.axis=cex.axis, ...) 
              axis(if(horiz) 2 else 1, at=at.l, labels=names.arg, 
                   lty=axis.lty, cex.axis=cex.names, ...) 
            } else {
              axis(if(horiz) 1 else 2, pos=0, cex.axis=cex.axis, ...)
              axis(if(horiz) 2 else 1, pos=ylim[1], at=at.l, labels=names.arg, 
                   cex.axis=cex.names, ...)
            }
          }
          else if (type=="xyplot") {
            if(plot.grid) {
              axis(if(horiz) 1 else 2, at=grid, cex.axis=cex.axis, ...)
              axis(if(horiz) 2 else 1, cex.axis=cex.axis, )
            } else {
              if (horiz) {
                axis(1, pos=ylim[1], cex.axis=cex.axis, ...)
                axis(2, pos=xlim[1], cex.axis=cex.axis, ...)
              } else {
                axis(2, pos=xlim[1], cex.axis=cex.axis, ...)
                axis(1, pos=ylim[1], cex.axis=cex.axis, ...)
              }
            }
          }
        }
      
      invisible(w.m)
      
    }
  
    else w.m
}

panel.overplot <- function(formula, data, subset, col, lty, ...)
  {
    m <- match.call()

    m[[1]] <- graphics:::plot.formula
    eval(m, parent.frame() )

    m[[1]] <- as.name('lowess.formula')
    tmp <- eval(m, parent.frame() )
          
    lines( tmp, col=col, lwd=2, lty=lty )
  }

overplot <- function (formula, data = parent.frame(),
                      same.scale=FALSE,
                      xlab, ylab,
                      xlim, ylim,
                      min.y, max.y,
                      log='',
                      panel='panel.overplot',
                      subset,
                      plot=TRUE, groups,
                      main,
                      f=2/3,
                      ... )
{
  ###
  # check that the formula had the right form
  ###
  if(
     length(formula)!=3 ||
     length(formula[[3]]) != 3  ||
     formula[[3]][[1]] != as.name("|")
     )
    stop("Formula must be of the form y ~ x1 | x2")

  if(!missing(subset))
    {
      flag <- eval(substitute(subset), envir=data)
      data <- data[flag,]
    }
  
  ###
  # Get the actual formula values
  ###
  cond <- eval(formula[[3]][[3]], envir=data, parent.frame())
  x <- eval(formula[[3]][[2]], envir=data, parent.frame())
  y <- eval(formula[[2]], envir=data, parent.frame())

  #print(data.frame(cond,x,y))
  
  y.all.min <- min(y, na.rm=TRUE)
  y.all.max <- max(y, na.rm=TRUE)

  x.all.min <- min(x, na.rm=TRUE)
  x.all.max <- max(x, na.rm=TRUE)

  if(y.all.min==y.all.max) browser()
  
  if (length(cond) == 0) {
    cond <- list(as.factor(rep(1, length(x))))
  }
  if (!is.factor(cond))
    {
      cond <- factor(cond)
    }


  ###
  # create a new call to the requested function
  ###
  mycall <- match.call(expand.dots=FALSE)
  mycall$panel <- mycall$plot <- mycall$groups <- mycall$same.scale <- NULL
  mycall$min.y <- mycall$max.y <- NULL
  mycall$data <- data
    
  # remove condition from formula
  mycall$formula[3] <- formula[[3]][2]

  # function name
  if(is.character(panel))
    panel <- as.name(panel)
  else
    panel <- deparse(substitute(panel))
  
  mycall[[1]] <- as.name(panel)

  # ylim
  if(same.scale)
    {
      if(missing(ylim))
        mycall$ylim <- range(y[y>0],na.rm=TRUE)
    }

  # xlim is always the same for all graphs
  if(missing(xlim))
    if(log %in% c("x","xy"))
      mycall$xlim <- range(x[x>0],na.rm=TRUE)
    else
      mycall$xlim <- range(x,na.rm=TRUE)


  
  ###
  # Only plot groups with non-na entries
  ##
  tmp <- na.omit(data.frame(x,y,cond))
  leveln <- sapply( split(tmp, tmp$cond), nrow )

  if(missing(groups)) groups <- names(leveln)
  ngroups <- length(groups)

  if(!missing(min.y) && length(min.y==1))
    min.y <- rep(min.y, length=ngroups)
  
  if(!missing(max.y) && length(max.y==1))
    max.y <- rep(max.y, length=ngroups)
  
  ###
  # Set up the plot regions
  ###
  oldpar <- par()['mar']
  on.exit(par(oldpar))
  
  par(mar=par("mar") +c(0,ngroups*2.5,0,0))
  
  ###
  # Ok. Now iterate over groups
  ## 

    i <- 1
    for(level in groups)
      {

        if(i>1)
          {
            par(new=TRUE)
          }

        mycall$subset <- (cond==level)

        mycall$ylab <- ''
        mycall$xlab <- ''
        mycall$xaxt = 'n'
        mycall$yaxt = 'n'
        mycall$pch = i
        mycall$col = i
        mycall$lty = i

        tmp.y <- y[mycall$subset & cond==level]
        min.tmp.y <- min(tmp.y, na.rm=TRUE)
        max.tmp.y <- max(tmp.y, na.rm=TRUE)
                
        if( !missing(min.y) || !missing(max.y) )
          {
            if(!missing(min.y) && !missing(max.y) )
              {
                mycall$ylim <- c(min.y[i], max.y[i])                
              }
            else if(missing(min.y) && !missing(max.y))
              {
                if(same.scale)
                  mycall$ylim <- c(y.all.min, max.y[i])
                else
                  mycall$ylim <- c(min.tmp.y, max.y[i])
              }
            else # !missing(min.y) && missing(max.y)
              {
                if(same.scale)
                  mycall$ylim <- c(min.y[i], y.all.max)
                else
                  mycall$ylim <- c(min.y[i], max.tmp.y)
              }
          }

        if(plot)
          {
            status <- try(
                eval(mycall, parent.frame())
                )
            if('try-error' %in% class(status))
              break;
            
          }

        usr <- par("usr")

        axis(side=2,line=2.5*(i-1),lty=i,col=i,lwd=2) 
        title(ylab=level, line=2.5*(i-1))
        
        if(i == 1)
          axis(side=1)

        i <- i+1
        
      }


  if(missing(xlab))
      xlab <- as.character(formula[[3]][[2]])
  
  if(missing(ylab))
      ylab <- as.character(formula[[2]])

  if(missing(main))
    main <- paste("plot of ", xlab, "vs.", ylab, "by",
                  as.character(formula[[3]][[3]]))
  
  if(i>1)
    {
      title(main=main,xlab=xlab)
      title(ylab=ylab, line=2.5*(i-1))
    }

  par(oldpar)
  
  invisible( split(data.frame(x,y),cond) )
}


permute <- function(x) sample( x, size=length(x), replace=FALSE )

#!!! BROKEN !!!

#plot.lm <-
#function(x, which = 1:5,
#         caption = c("Residuals vs Fitted", "Normal Q-Q plot",
#           "Scale-Location plot", "Cook's distance plot"),
#         panel = panel.smooth,
#         sub.caption = deparse(x$call), main = "",
#         ask = interactive() && nb.fig < length(which)
#         	&& .Device != "postscript",
#         ...,
#         id.n = 3, labels.id = names(residuals(x)), cex.id = 0.75,
#         band=TRUE,rug=TRUE)
#{
#    if (!inherits(x, "lm"))
#	stop("Use only with 'lm' objects")
#    show <- rep(FALSE, 5)
#    if(!is.numeric(which) || any(which < 1) || any(which > 5))
#        stop("`which' must be in 1:5")
#    show[which] <- TRUE
#    r <- residuals(x)
#    n <- length(r)
#    yh <- predict(x) # != fitted() for glm
#    if (any(show[2:4]))
#        s <- if(inherits(x, "rlm")) x$s else sqrt(deviance(x)/df.residual(x))
#    if (any(show[2:3])) {
#        ylab23 <- if(inherits(x, "glm"))
#          "Std. deviance resid." else "Standardized residuals"
#        hii <- lm.influence(x)$hat
#        w <- weights(x)
#        # r.w := weighted.residuals(x):
#        r.w <- if(is.null(w)) .Alias(r) else (sqrt(w)*r)[w!=0]
#        rs <- r.w/(s * sqrt(1 - hii))
#    }
#    if (any(show[c(1,3)]))
#        l.fit <- if(inherits(x,"glm"))
#            "Predicted values" else "Fitted values"
#    if (is.null(id.n))
#	id.n <- 0
#    else {
#	id.n <- as.integer(id.n)
#	if(id.n < 0 || id.n > n)
#	    stop(paste("`id.n' must be in { 1,..,",n,"}"))
#    }
#    if(id.n > 0) {
#        if(is.null(labels.id))
#            labels.id <- paste(1:n)
#        iid <- 1:id.n
#	show.r <- order(-abs(r))[iid]
#        if(any(show[2:3]))
#            show.rs <- order(-abs(rs))[iid]
#        text.id <- function(x,y, ind, adj.x = FALSE)
#            text(x - if(adj.x) strwidth(" ")*cex.id else 0, y, labels.id[ind],
#                 cex = cex.id, xpd = TRUE, adj = if(adj.x) 1)
#    }
#    nb.fig <- prod(par("mfcol"))
#    one.fig <- prod(par("mfcol")) == 1
#    if (ask) {
#	op <- par(ask = TRUE)
#	on.exit(par(op))
#    }
#    ##---------- Do the individual plots : ----------
#    if (show[1]) {
#	ylim <- range(r)
#	if(id.n > 0)
#	    ylim <- ylim + c(-1,1)* 0.08 * diff(ylim)
#	plot(yh, r, xlab = l.fit, ylab = "Residuals", main = main,
#	     ylim = ylim, type = "n", ...)
#	panel(yh, r, ...)
#        if(rug)  rug(yh)                    ## GRW 2001-06-08
#        if(band) bandplot(yh,r,add=TRUE,lty=3) ## GRW 2001-06-08
#	if (one.fig)
#	    title(sub = sub.caption, ...)
#	mtext(caption[1], 3, 0.25)
#	if(id.n > 0) {
#	    y.id <- r[show.r]
#	    y.id[y.id < 0] <- y.id[y.id < 0] - strheight(" ")/3
#	    text.id(yh[show.r], y.id, show.r, adj.x = TRUE)
#	}
#	abline(h = 0, lty = 3, col = "gray")
#    }
#    if (show[2]) {
#	ylim <- range(rs)
#	ylim[2] <- ylim[2] + diff(ylim) * 0.075
#	qq <- qqnorm(rs, main = main, ylab = ylab23, ylim = ylim, ...)
#        qqline(rs)
#	if (one.fig)
#	    title(sub = sub.caption, ...)
#	mtext(caption[2], 3, 0.25)
#	if(id.n > 0)
#	    text.id(qq$x[show.rs], qq$y[show.rs], show.rs, adj.x = TRUE)
#    }
#    if (show[3]) {
#	sqrtabsr <- sqrt(abs(rs))
#	ylim <- c(0, max(sqrtabsr))
#	yl <- as.expression(substitute(sqrt(abs(YL)), list(YL=as.name(ylab23))))
#        yhn0 <- if(is.null(w)) .Alias(yh) else yh[w!=0]
#	plot(yhn0, sqrtabsr, xlab = l.fit, ylab = yl, main = main,
#	     ylim = ylim, type = "n", ...)
#	panel(yhn0, sqrtabsr, ...)

#        abline(h=mean(sqrtabsr),lty = 3, col = "gray")
#        if(rug)  rug(yh)                             ## GRW 2001-06-08
#        if(band) bandplot(yhn0,sqrtabsr,add=TRUE,lty=3) ## GRW 2001-06-08

#	if (one.fig)
#	    title(sub = sub.caption, ...)
#	mtext(caption[3], 3, 0.25)
#	if(id.n > 0)
#	    text.id(yhn0[show.rs], sqrtabsr[show.rs], show.rs, adj.x = TRUE)
#    }
#    if (show[4]) {
#	cook <- cooks.distance(x, sd=s)
#	if(id.n > 0) {
#	    show.r <- order(-cook)[iid]# index of largest `id.n' ones
#	    ymx <- cook[show.r[1]] * 1.075
#	} else ymx <- max(cook)
#	plot(cook, type = "h", ylim = c(0, ymx), main = main,
#	     xlab = "Obs. number", ylab = "Cook's distance", ...)
#	if (one.fig)
#	    title(sub = sub.caption, ...)
#	mtext(caption[4], 3, 0.25)
#	if(id.n > 0)
#	    text.id(show.r, cook[show.r] + 0.4*cex.id * strheight(" "), show.r)
#    }
#    if (show[5])
#      {
#        ## plot residuals against each predictor ##
#        data  <- model.frame(x)
#        for( i in 2:ncol(data) )
#          {
#            plot( data[,i], r, xlab=names(data)[i], ylab="Residuals", type="n")
#            panel( data[,i], r, ... )
#            if(rug)  rug(data[,i])
#            if(band) bandplot(data[,i],r,add=TRUE,lty=3)
#            abline(h=0,lty = 3, col = "gray")
#          }
#      }
#    if (!one.fig && par("oma")[3] >= 1)
#	mtext(sub.caption, outer = TRUE, cex = 1.25)
#    invisible()
#}
# $Id: plotCI.R,v 1.9 2004/05/24 23:43:46 warnes Exp $
#
# $Log: plotCI.R,v $
# Revision 1.9  2004/05/24 23:43:46  warnes
# Modified to use invalid() to check arguments instead of missing().
# This fixes some build errors under R-1.9.0-Patched.
#
# Revision 1.8  2003/12/02 16:54:40  warnes
#
# - Add '...' parameter to call to text to allow user to control size/color/etc.
#
# Revision 1.7  2002/04/09 00:51:30  warneg
#
# Checkin for version 0.5.3
#
# Revision 1.6  2002/03/20 04:16:32  warneg
# - Changes to add compatibility with S-Plus 2000.
#
# Revision 1.5  2002/02/04 19:12:27  warneg
#
# - When err="x", 'col' was being used to plot bars, rather than 'barcol'.
#
# Revision 1.4  2002/02/04 19:09:53  warneg
#
# - fixed typo, when err="x", lty was 'slty' causing an error.
#
# Revision 1.3  2001/10/16 23:00:19  warneg
#
# - Added minbar and maxbar parameters
# - Added cvs id and log tags to header
#
#


plotCI <- function (x, y = NULL,
                    uiw, liw = uiw,   
                    ui, li, 
                    err='y', 
                    col=par("col"),
                    ylim=NULL,
                    xlim=NULL,
                    barcol=col,
                    sfrac = 0.01,
                    gap=1,
                    lwd=par("lwd"),
                    lty=par("lty"),
                    labels=FALSE,
                    add=FALSE,
                    xlab,
                    ylab,
                    minbar,
                    maxbar,
                    ...
                    )
{

  if (is.list(x)) { 
    y <- x$y 
    x <- x$x 
  }

  if(invalid(xlab))
    xlab <- deparse(substitute(x))
  
  if(invalid(ylab))
    {
      if(is.null(y))
        {
          xlab  <- ""
          ylab <- deparse(substitute(x))
        }
      else
        ylab <- deparse(substitute(y))
    }

  if (is.null(y)) { 
    if (is.null(x)) 
      stop("both x and y NULL") 
    y <- as.numeric(x) 
    x <- seq(along = x) 
  }

  
  if(err=="y")
    z  <- y
  else
    z  <- x
  
  if(invalid(ui))
    ui <- z + uiw
  if(invalid(li)) 
    li <- z - liw

  if(!invalid(minbar))
    li <- ifelse( li < minbar, minbar, li)

  if(!invalid(maxbar))
    ui <- ifelse( ui > maxbar, maxbar, ui)
  
   if(err=="y")
     {
       if(is.null(ylim))
         ylim <- range(c(y, ui, li), na.rm=TRUE)
       if(is.null(xlim) && !is.R() )
         xlim <- range( x, na.rm=TRUE)
     }
   else if(err=="x")
     {
       if(is.null(xlim))
         xlim <- range(c(x, ui, li), na.rm=TRUE)
       if(is.null(ylim) && !is.R() )
         ylim <- range( x, na.rm=TRUE)
     }
    
  if(!add)
    {
      if(invalid(labels) || labels==FALSE )
        plot(x, y, ylim = ylim, xlim=xlim, col=col,
             xlab=xlab, ylab=ylab, ...)
      else
        {
          plot(x, y, ylim = ylim, xlim=xlim, col=col, type="n",
               xlab=xlab, ylab=ylab,  ...)
          text(x, y, label=labels, col=col, ... )
        }
    }
  if(is.R())
    myarrows <- function(...) arrows(...)
  else
    myarrows <- function(x1,y1,x2,y2,angle,code,length,...)
      {
        segments(x1,y1,x2,y2,open=TRUE,...)
        if(code==1)
          segments(x1-length/2,y1,x1+length/2,y1,...)
        else
          segments(x2-length/2,y2,x2+length/2,y2,...)
      }
 
  if(err=="y")
    {
      if(gap!=FALSE)
        gap <- strheight("O") * gap
      smidge <- par("fin")[1] * sfrac

      
      # draw upper bar
      if(!is.null(li))
          myarrows(x , li, x, pmax(y-gap,li), col=barcol, lwd=lwd,
                 lty=lty, angle=90, length=smidge, code=1)
      # draw lower bar
      if(!is.null(ui))
          myarrows(x , ui, x, pmin(y+gap,ui), col=barcol,
                 lwd=lwd, lty=lty, angle=90, length=smidge, code=1)
    }
  else
    {
      if(gap!=FALSE)
        gap <- strwidth("O") * gap
      smidge <- par("fin")[2] * sfrac

      # draw left bar
      if(!is.null(li))
        myarrows(li, y, pmax(x-gap,li), y, col=barcol, lwd=lwd,
                 lty=lty, angle=90, length=smidge, code=1)
      if(!is.null(ui))
        myarrows(ui, y, pmin(x+gap,ui), y, col=barcol, lwd=lwd,
                 lty=lty, angle=90, length=smidge, code=1)
      
    }
      
    

invisible(list(x = x, y = y)) 
} 
# $Id: plotmeans.R,v 1.14 2004/05/24 23:43:46 warnes Exp $
#
# $Log: plotmeans.R,v $
# Revision 1.14  2004/05/24 23:43:46  warnes
# Modified to use invalid() to check arguments instead of missing().
# This fixes some build errors under R-1.9.0-Patched.
#
# Revision 1.13  2003/11/10 22:11:13  warnes
#
# - Add files contributed by Arni Magnusson
#   <arnima@u.washington.edu>. As well as some of my own.
#
# Revision 1.12  2003/04/22 17:28:56  warnes
#
# - Fixeed warning messing caused when 'connect' is a vector.
#
# Revision 1.11  2002/09/24 14:53:56  warnes
#
# - Changed digits=options("digits") which produces a list of length 1
#   to digits=getOption("digits") which returns a vector of length one.
#   The former was causing an error when passed to round().
#
# Revision 1.10  2002/04/09 00:51:30  warneg
#
# Checkin for version 0.5.3
#
# Revision 1.9  2002/03/20 04:17:05  warneg
# - Changes to add compatibility with S-Plus 2000
#
# Revision 1.8  2002/03/05 16:44:24  warneg
# - Replace "TRUE" with "TRUE".  Problems arrive when there is a variable named "TRUE".
#
# Revision 1.7  2001/12/05 19:49:29  warneg
# - Added ability to use the t-distribution to compute confidence
#   intervals.  This is controlled using the 'use.t' parameter.
#
# Revision 1.6  2001/10/26 02:48:19  warneg
#
# Added correct handling of 'xaxt="n"'.
#
# Revision 1.5  2001/10/16 22:59:27  warneg
#
# Added cvs id and log tag to file header
#
#


# Plot means with confidence intervals for groups defined by right
# side of formulae
#
# example:
#
# data  <-  data.frame(y=rnorm(100), x=factor(rep(c("A","C","FALSE","Z"),25)))
# plotmeans( y ~ x, data=data, connect=FALSE )

plotmeans  <- function (formula, data = NULL, subset, na.action,
                         bars=TRUE, p=0.95,
                         minsd=0, minbar=NULL, maxbar=NULL,
                         xlab=names(mf)[2], ylab=names(mf)[1],
                         mean.labels=FALSE, ci.label=FALSE, n.label=TRUE,
                         digits=getOption("digits"), col="black",
                         barwidth=1,
                         barcol="blue",
                         connect=TRUE,
                         ccol=col,
                         legends=names(means),
                         xaxt,
                         use.t = TRUE,
                         ...)
{
  is.R <- get("is.R")
  if(is.null(is.R)) is.R <- function(x) FALSE

  if(!is.R())
    {
      if(col=="black")
        col <- 1
      if(barcol=="blue")
        barcol <- 2
    }
  
    if (invalid(formula) || (length(formula) != 3)) 
        stop("formula missing or incorrect")
    if (invalid(na.action)) 
        na.action <- options("na.action")
    m <- match.call(expand.dots = FALSE)
    if(is.R())
      {
        if (is.matrix(eval(m$data, parent.frame()))) 
          m$data <- as.data.frame(data)
      }
    else
      {
        if (is.matrix(eval(m$data, FALSE)))
          m$data <- as.data.frame(data)
      }
    m$... <- m$bars <- m$barcol <- m$p <- NULL
    m$minsd <- m$minbar <- m$maxbar <- NULL
    m$xlab <- m$ylab <-  NULL
    m$col  <- m$barwidth  <- NULL
    m$digits  <- m$mean.labels  <- m$ci.label  <- m$n.label <- NULL
    m$connect  <- m$ccol  <-  m$legends <- m$labels<- NULL
    m$xaxt <- m$use.t <- NULL
    m[[1]] <- as.name("model.frame")
    mf <- eval(m, parent.frame())
    response <- attr(attr(mf, "terms"), "response")
    means  <-  sapply(split(mf[[response]], mf[[-response]]), mean, na.rm=TRUE)
    xlim  <-  c(0.5, length(means)+0.5)
    
    if(!bars)
      {
        plot( means, ..., col=col, xlim=xlim)
      }
    else
      {

    myvar  <-  function(x) var(x[!is.na(x)])
        
    vars <- sapply(split(mf[[response]], mf[[-response]]), myvar)
    ns   <- sapply( sapply(split(mf[[response]], mf[[-response]]), na.omit,
                           simplify=FALSE), length )

    # apply minimum variance specified by minsd^2
    vars <- ifelse( vars < (minsd^2), (minsd^2), vars)

    if(use.t)
      ci.width  <- qt( (1+p)/2, ns-1 ) * sqrt(vars/ns)
    else
      ci.width  <- qnorm( (1+p)/2 ) * sqrt(vars/ns)

    if(length(mean.labels)==1)
      {
        if (mean.labels==TRUE)
          mean.labels  <-  format( round(means, digits=digits ))
        else if (mean.labels==FALSE)
          mean.labels  <- NULL
      }

    plotCI(x=1:length(means), y=means, uiw=ci.width, xaxt="n",
           xlab=xlab, ylab=ylab, labels=mean.labels, col=col, xlim=xlim,
           lwd=barwidth, barcol=barcol, minbar=minbar, maxbar=maxbar, ... )
    if(invalid(xaxt) || xaxt!="n")
      axis(1, at = 1:length(means), labels = legends)
    
    if(ci.label)
      {
        ci.lower <- means-ci.width
        ci.upper <- means+ci.width

        if(!invalid(minbar))
          ci.lower <- ifelse(ci.lower < minbar, minbar, ci.lower)
        if(!invalid(maxbar))
          ci.upper <- ifelse(ci.upper > maxbar, maxbar, ci.upper)
        
        labels.lower <- paste( " \n", format(round(ci.lower, digits=digits)),
                              sep="")
        labels.upper <- paste( format(round(ci.upper, digits=digits)), "\n ",
                              sep="")

        text(x=1:length(means),y=ci.lower, labels=labels.lower, col=col)
        text(x=1:length(means),y=ci.upper, labels=labels.upper, col=col)
      }
    
  }
    
    
    if(n.label)
      if(is.R())
        text(x=1:length(means),y=par("usr")[3],
             labels=paste("n=",ns,"\n",sep=""))
      else
        {
          axisadj <- (par("usr")[4] - (par("usr")[3]) )/75
          text(x=1:length(means),y=par("usr")[3] + axisadj,
               labels=paste("n=",ns,"\n",sep=""))
        }
    
    if(!invalid(connect) & !identical(connect,FALSE))
      {
        if(is.list(connect))
          {
            if(length(ccol)==1)
              ccol  <-  rep(ccol, length(connect) )
            
            for(which in 1:length(connect))
              lines(x=connect[[which]],y=means[connect[[which]]],col=ccol[which])
          }
        else  
          lines(means, ..., col=ccol)
      }
    

    
}

# $Id: qqnorm.aov.R,v 1.3 2003/03/07 15:48:35 warnes Exp $
#
# $Log: qqnorm.aov.R,v $
# Revision 1.3  2003/03/07 15:48:35  warnes
#
# - Minor changes to code to allow the package to be provided as an
#   S-Plus chapter.
#
# Revision 1.2  2003/01/02 15:40:09  warnes
# - Renamed first parameter to match qqnorm generic.
#
# Revision 1.1  2002/12/31 19:50:58  warnes
# Initial checkin of qqnorm.aov function and documentation submitted by
# Kjetil Halvorsen <kjetilh@jupiter.umsanet.edu.bo>.
#
#

if(is.R())
qqnorm.aov <- function (y, full = FALSE, label = FALSE, omit = NULL,
                        xlab = paste(if(full) "" else "Half", " Normal plot"),
                        ylab = "Effects", ...)
{
    r <- y$rank
    eff <- if (full)
        effects(y, set.sign = TRUE)[1:r]
    else abs(effects(y))[1:r]
    na <- names(eff)
    int <- match("(Intercept)", na)
    if (!is.null(omit)) {
        if (is.character(omit)) {
            int <- c(int, match(omit, na))
        }
        else int <- c(int, omit)
    }
    int <- int[!is.na(int)]
    if (length(int))
        eff <- eff[-int]
    n <- length(eff)
    if (n <= 0)
        stop("Not enough effects")
    ord <- order(eff)
    na <- names(eff)
    P <- if (full)
        ppoints(n)
    else ((1:n) + n)/(2 * n + 1)
    Q <- qnorm(P)
    plot(x = Q, y = eff[ord], xlab = xlab, ylab = ylab, ...)
    if (label && dev.interactive())
        identify(Q, eff[ord], names(eff)[ord])
}


# $Id: quantcut.R,v 1.1 2002/03/26 14:29:08 warneg Exp $
#
# $Log: quantcut.R,v $
# Revision 1.1  2002/03/26 14:29:08  warneg
#
# Initial checkin.
#
#

quantcut <- function(x, q=c(0,1/4,2/4,3/4,4/4), na.rm=TRUE,
                     include.lowest=TRUE, ... )
  {
    quant <- quantile(x, q, na.rm=na.rm)
    cut( x, quant, include.lowest=include.lowest,  ... )
  }
read.xls <- function(xls, sheet = 1, verbose=FALSE, ...)
{
  ###
  # directories
  package.dir <- .path.package('gregmisc')
  perl.dir <- file.path(package.dir,'perl')
  #
  ###

  ###
  # files
  xls <- dQuote(xls) # dQuote in case of spaces in path
  xls2csv <- file.path(perl.dir,'xls2csv.pl')
  csv <- paste(tempfile(), "csv", sep = ".")
  #
  ###

  ###
  # execution command
  cmd <- paste("perl", xls2csv, xls, dQuote(csv), sheet, sep=" ")
  #
  ###

  ###
  # do the translation
  if(verbose)  cat("Executing ", cmd, "... \n")
  #
  results <- system(cmd, intern=!verbose)
  #
  if (verbose) cat("done.\n")
  #
  ###
  
  # now read the csv file
  out <- read.csv(csv, ...)

  # clean up
  file.remove(csv)
  
  return(out)
}
# $Id: rename.vars.R,v 1.5 2004/04/01 20:23:14 warnes Exp $
#
# $Log: rename.vars.R,v $
# Revision 1.5  2004/04/01 20:23:14  warnes
# Add function remove.vars().
#
# Revision 1.4  2002/04/09 00:51:30  warneg
#
# Checkin for version 0.5.3
#
# Revision 1.3  2001/12/08 01:54:19  warneg
# Changed 'TRUE' to 'TRUE' in parameter list.
#
# Revision 1.2  2001/12/07 22:55:32  warneg
#
# Added attribution.
#
# Revision 1.1  2001/12/07 21:40:54  warneg
#
# Initial checkin
#
#

# submitted by Don MacQueen <macq@llnl.gov>

rename.vars <- function(data,from='',to='',info=TRUE) {

   dsn <- deparse(substitute(data))
   dfn <- names(data)

   if ( length(from) != length(to)) {
     cat('--------- from and to not same length ---------\n')
     stop()
   }

   if (length(dfn) < length(to)) {
     cat('--------- too many new names ---------\n')
     stop()
   }

   chng <- match(from,dfn)

   frm.in <- from %in% dfn
   if (!all(frm.in) ) {
     cat('---------- some of the from names not found in',dsn,'\n')
     stop()
   }

   if (length(to) != length(unique(to))) {
     cat('---------- New names not unique\n')
     stop()
   }

   dfn.new <- dfn
   dfn.new[chng] <- to
   if (info) cat('\nChanging in',dsn)
   tmp <- rbind(from,to)
   dimnames(tmp)[[1]] <- c('From:','To:')
   dimnames(tmp)[[2]] <- rep('',length(from))
   if (info)
     {
       print(tmp,quote=FALSE)
       cat("\n")
     }
   names(data) <- dfn.new
   data
}


# GRW 2004-04-01
remove.vars <- function( data, names, info=TRUE)
  {
    for( i in names )
      {
        if(info)
          cat("Removing variable '", i, "'\n", sep="")
        data[[i]] <- NULL
      }
    data
  }
# $Id: reorder.R,v 1.4 2004/01/21 12:06:26 warnes Exp $
#
# $Log: reorder.R,v $
# Revision 1.4  2004/01/21 12:06:26  warnes
# - Add ... argument to match generic provided in mva.
#
# Revision 1.3  2003/04/22 15:42:33  warnes
#
# - The mva package (which is part of recommended) now provides a
#   generic 'reorder' function.  Consequently, the 'reorder' function
#   here has been renamed to 'reorder.factor'.
#
# - Removed check of whether the argument is a factor object.
#
# Revision 1.2  2003/03/03 17:24:21  warnes
#
# - Added handling of factor level names in addition to numeric indexes.
#
# Revision 1.1  2002/08/01 18:06:41  warnes
#
# Added reorder() function to reorder the levels of a factor.
#
#

# Reorder the levels of a factor.

reorder.factor <- function( x, order, ... )
  {
    if(is.numeric(order))
      factor( x, levels=levels(x)[order] )
    else
      factor( x, levels=order )
  }
# $Id: residplot.R,v 1.3 2002/09/23 13:59:30 warnes Exp $
#
# $Log: residplot.R,v $
# Revision 1.3  2002/09/23 13:59:30  warnes
# - Modified all files to include CVS Id and Log tags.
#
#

residplot  <-  function(model, formula, ...)
  {
    data  <- expand.model.frame( model, formula, na.expand=TRUE)

    newform  <- eval(parse( text=paste("as.call(", "resid(model) ~",
                        formula[-1],")" )))
    
    plot( newform, data=data, ylab="Residuals")
    lines(lowess( newform, data=data ), col="red")
    bandplot(newform,data=data)
  }
              
rich.colors <- function(n, palette=c("temperature", "blues"),
                        rgb.matrix=FALSE, plot.colors=FALSE)
{
  if(n <= 0)
    return(character(0))

  palette <- match.arg(palette)
  x <- seq(0, 1, length=n)

  if(palette == "temperature") {
    r <- 1 / (1+exp(20-35*x))
    g <- pmin(pmax(0,-0.8+6*x-5*x^2), 1)
    b <- dnorm(x,0.25,0.15) / max(dnorm(x,0.25,0.15))
  }
  else {
    r <-        0.6*x + 0.4*x^2
    g <-        1.5*x - 0.5*x^2
    b <- 0.36 + 2.4*x - 2.0*x^2
    b[x>0.4] <- 1
  }

  rgb.m <- matrix(c(r,g,b), ncol=3,
                  dimnames=list(as.character(seq(length=n)),
                    c("red","green","blue")))
  rich.vector <- apply(rgb.m, 1, function(v) rgb(v[1],v[2],v[3]))

  if(rgb.matrix)
    attr(rich.vector, "rgb.matrix") <- rgb.m

  if(plot.colors) {
    opar <- par("fig", "plt")
    par(fig=c(0,1,0,0.7), plt=c(0.15,0.9,0.2,0.95))
    plot(NA, xlim=c(-0.01,1.01), ylim=c(-0.01,1.01), xlab="Spectrum",
         ylab="", xaxs="i", yaxs="i", axes=FALSE)
    title(ylab="Value", mgp=c(3.5,0,0))
    matlines(x, rgb.m, col=colnames(rgb.m), lty=1, lwd=3)
    matpoints(x, rgb.m, col=colnames(rgb.m), pch=16)
    axis(1, at=0:1)
    axis(2, at=0:1, las=1)
    par(fig=c(0,1,0.75,0.9), plt=c(0.08,0.97,0,1), new=TRUE)
    midpoints <- barplot(rep(1,n), col=rich.vector, border=FALSE,
                         space=FALSE, axes=FALSE)
    axis(1, at=midpoints, labels=1:n, lty=0, cex.axis=0.6)
    par(opar)
  }
  return(rich.vector)
}
# $Id: running.R,v 1.9 2004/04/27 14:33:15 warnes Exp $
#

"running" <- function(X, Y=NULL,
                      fun=mean,
                      width=min(length(X), 20),
                      allow.fewer=FALSE, pad=FALSE,
                      align=c("right", "center", "left"),
                      simplify=TRUE,
                      ...)
{
  align=match.arg(align)
  
  n <- length(X)

  if (align=="left")
    {
      from  <-  1:n 
      to    <-  pmin( (1:n)+width-1, n)
    }
  else if (align=="right")
    {
      from  <-  pmax( (1:n)-width+1, 1)
      to    <-  1:n 
    }
  else #align=="center"
    {
      from <-  pmax((2-width):n,1)
      to   <-  pmin(1:(n+width-1),n)
      if(!odd(width)) stop("width must be odd for center alignment")
      
    }

  elements  <- apply(cbind(from, to), 1, function(x) seq(x[1], x[2]) )

  if(is.matrix(elements))
    elements <- as.data.frame(elements) # ensure its a list!

  names(elements) <- paste(from,to,sep=':')

  if(!allow.fewer)
    {
      len <- sapply(elements, length)
      skip <- (len < width)
    }
  else
    {
      skip <- 0
    }
  

  run.elements  <- elements[!skip]
  
  if(is.null(Y))  # univariate 
    {
      funct <- function(which,what,fun,...) fun(what[which],...)

      if(simplify)
        Xvar <- sapply(run.elements, funct, what=X, fun=fun, ...)
      else
        Xvar <- lapply(run.elements, funct, what=X, fun=fun, ...)        
    }
  else # bivariate
    {
      funct <- function(which,XX,YY,fun,...) fun(XX[which],YY[which], ...)
      
      if(simplify)
        Xvar <- sapply(run.elements, funct, XX=X, YY=Y, fun=fun, ...)
      else
        Xvar <- lapply(run.elements, funct, XX=X, YY=Y, fun=fun, ...)
    }

  
  if(allow.fewer || !pad)
      return(Xvar)

  if(simplify)
    if(is.matrix(Xvar))
      {
        wholemat <- matrix( new(class(Xvar[1,1]), NA),
                           ncol=length(to), nrow=nrow(Xvar) )
        colnames(wholemat) <- paste(from,to,sep=':')
        wholemat[,-skip] <- Xvar
        Xvar <- wholemat
      }
    else
      {
        wholelist <- rep(new(class(Xvar[1]),NA),length(from))
        names(wholelist) <-  names(elements)
        wholelist[ names(Xvar) ] <- Xvar
        Xvar <- wholelist
      }
  
  return(Xvar)}

sinkplot <- function(operation=c("start","plot","cancel"),...)
  {
    operation <- match.arg(operation)

    if( operation=="start" )
      {
        if (exists(".sinkplot.conn", env=globalenv()) &&
            get(".sinkplot.conn", env=globalenv()) )
          stop("sinkplot already in force")

        
        .sinkplot.conn <- textConnection(".sinkplot.data", "w", local=FALSE)
        assign(x=".sinkplot.conn", value=.sinkplot.conn, envir=globalenv())

        on.exit(sink())
        sink(.sinkplot.conn)
        on.exit()
      }
    else
      {
        if (!exists(".sinkplot.conn", env=globalenv()) || !.sinkplot.conn )
          stop("No sinkplot currently in force")

        sink()

        data <- get(".sinkplot.data", env=globalenv())
        
        if( operation=="plot" )
            textplot( paste( data, collapse="\n"), ... )
        
        close(get(".sinkplot.conn", env=globalenv()))

        if(exists(".sinkplot.data", env=globalenv()))
          rm(".sinkplot.data", pos=globalenv())

        if(exists(".sinkplot.conn", env=globalenv()))
          rm(".sinkplot.conn", pos=globalenv())

        invisible(data)
      }
  }
smartlegend <- function(x=c("left","center","right"),
                        y=c("top","center","bottom"),
                        ..., inset=0.05 )
  {
    x <- match.arg(x)
    y <- match.arg(y)

    usr <- par("usr")
    inset.x <- inset * (usr[2] - usr[1])
    inset.y <- inset * (usr[4] - usr[3])
    
    if(x=="left")
      {
        x.pos <- usr[1] + inset.x
        xjust = 0
      }
    else if(x=="center")
      {
        x.pos <- (usr[1] + usr[2])/2
        xjust = 0.5
      }
    else # y=="right"
      {
        x.pos <- usr[2] - inset.x
        xjust = 1
      }

    if(y=="bottom")
      {
        y.pos <- usr[3] + inset.y
        yjust = 0
      }
    else if(y=="center")
      {
        y.pos <- (usr[3] + usr[4])/2
        yjust = 0.5
      }
    else
      {
        y.pos <- usr[4] - inset.y
        yjust = 1
      }

    
    if(par("xlog")) x.pos <- 10^x.pos
    if(par("ylog")) y.pos <- 10^y.pos 
    
    legend( x=x.pos, y=y.pos, ..., xjust=xjust, yjust=yjust)
  }


# $Id: space.R,v 1.6 2004/04/13 13:42:19 warnes Exp $
#
# $Log: space.R,v $
# Revision 1.6  2004/04/13 13:42:19  warnes
# Add ability to space points along 'y' direction.
#
# Revision 1.5  2003/11/10 22:11:13  warnes
#
# - Add files contributed by Arni Magnusson
#   <arnima@u.washington.edu>. As well as some of my own.
#
# Revision 1.4  2002/04/09 00:51:31  warneg
#
# Checkin for version 0.5.3
#
# Revision 1.3  2002/02/04 23:20:29  warneg
#
# - Add na.rm parameter and make the default TRUE.
#
# Revision 1.2  2001/12/07 22:24:30  warneg
#
# - Added cvs tags.
#
#

# When there are two or more points with the same (x,y) value (or
# within x+-s[1] and x+-s[2]), space these out in the y direction so
# that the points are separated by at least distance s.


space <-  function(x,y,s=1/50, na.rm=TRUE, direction="x")
  {
    if(direction!='x')
      {
        tmp <- y
        y <- x
        x <- tmp
      }
    
    if(na.rm)
      {
        ind <- is.na(x) | is.na(y)
        x <- x[!ind]
        y <- y[!ind]
      }
    
    if (length(s)==1) s <- c(s,s)
    
    spacing.x <- (max(x) - min(x))*s[1]
    spacing.y <- (max(y) - min(y))*s[2]

    within <- function(x,y,delta) { abs(x-y) < delta }

    # sort x,y so we can do the work
    ord <- order(x,y)
    undo <- order(ord)
    
    x <- x[ord]
    y <- y[ord]

    startsame <- 1
    same.x <- x[1]
    same.y <- y[1]
    
    for( i in 1:length(x) )
      {
        if(i>1 &&
           within(x[i],same.x,spacing.x) &&
           within(y[i],same.y,spacing.y) )
          {
            if(x[startsame] == same.x )
              x[startsame] <- x[startsame] 

            cumrun <- i - startsame
            
            x[i] <- x[i] + (-1)^(cumrun+1) * floor((cumrun+1) /2) * spacing.x
          } else {
            startsame <- i
            same.x <- x[i]
            same.y <- y[i]
          }
      }

    if(direction!='x')
      {
        tmp <- y
        y <- x
        x <- tmp
      }

    return( list(x=x[undo], y=y[undo]) )
}    
      
# $Id: textplot.R,v 1.5 2004/03/30 19:04:53 warnes Exp $
#
# $Log: textplot.R,v $
# Revision 1.5  2004/03/30 19:04:53  warnes
# - Fix bug in textplot() reported by Wright, Kevin <kevin.d.wright@pioneer.com>.
#
# Revision 1.4  2004/03/26 22:16:44  warnes
# Misc changes.
#
# Revision 1.3  2004/01/21 04:31:25  warnes
# - Mark sprint() as depreciated.
# - Replace references to sprint with capture.output()
# - Use match.arg for halign and valign arguments to textplot.default.
# - Fix textplot.character so that a vector of characters is properly
#   displayed. Previouslt, character vectors were plotted on top of each
#   other.
#
# Revision 1.2  2003/04/04 13:41:42  warnes
#
#
# - Added textplot.character to handle character strings.
# - Moved test for vector and matrix arguments to textplot.default.
# - Renamed arguments "col.margin" and "row.margin" to "cmar", and
#   "rmar" so that specifying "col='red'" is possible.
#
# Revision 1.1  2003/04/02 22:29:53  warnes
#
# - Added textplot function and friends, as well as documentation.
#
#

textplot <- function(object, halign="center", valign="center", cex, ... )
  UseMethod('textplot')


textplot.default <- function(object,
                             halign=c("center","left","right"),
                             valign=c("center","top","bottom"),
                             cex, ... )
{

  if (is.matrix(object) || (is.vector(object) && length(object)>1) )
    return(textplot.matrix(object, halign, valign, cex, ... ))

  halign <- match.arg(halign)
  valign <- match.arg(valign)

  #if(is.character(object))
  #  object <- paste(object, collapse="\n")
  #else
  #  object <- capture.output(print(object))

  textplot.character(object, halign,  valign, cex, ...)
}


textplot.data.frame <- function(object,
                             halign=c("center","left","right"),
                             valign=c("center","top","bottom"),
                             cex, ... )
    textplot.matrix(object, halign, valign, cex, ... )  


textplot.matrix <- function(object,
                            halign=c("center","left","right"),
                            valign=c("center","top","bottom"),
                            cex, cmar=2, rmar=0.5,
                            show.rownames=TRUE, show.colnames=TRUE,
                            hadj=1,
                            vadj=1,
                            mar= c(0,0,3,0)+0.1,
                            ... )
{

  if(is.vector(object))
    object <- t(as.matrix(object))
  else
    object <- as.matrix(object)

  halign=match.arg(halign)
  valign=match.arg(valign)
  
  opar <- par()[c("mar","xpd","cex")]
  on.exit( par(opar) )
  mar=c(0,0,3,0)+0.1
  par(mar=mar, xpd=FALSE )
    
  # setup plot area
  plot.new()
  plot.window(xlim=c(0,1),ylim=c(0,1), log = "", asp=NA)


  
  # add 'r-style' row and column labels if not present
  if( is.null(colnames(object) ) )
    colnames(object) <- paste( "[,", 1:ncol(object), "]", sep="" )
  if( is.null(rownames(object)) )
    rownames(object) <- paste( "[", 1:nrow(object), ",]", sep="")

  
  # extend the matrix to include them
  if( show.rownames )
    {
      object <- cbind( rownames(object), object )
    }
  if( show.colnames )
    {
        
      object <- rbind( colnames(object), object )
    }

  # set the character size
  if( missing(cex) )
    {
      cex <- 1.0
      lastloop <- FALSE
    }
  else
    {
      lastloop <- TRUE
    }

  for (i in 1:20)
    {
      oldcex <- cex
      
      width  <- sum(
                    apply( object, 2,
                          function(x) max(strwidth(x,cex=1) ) )
                    ) +
                      strwidth('M', cex=cex) * cmar * ncol(object)
      
      height <- strheight('M', cex=cex) * nrow(object) * (1 + rmar)
      
      if(lastloop) break
      
      if (abs(oldcex - cex) < 0.001)
        {
          lastloop <- TRUE
          cex <- oldcex * (1/max(width, height))
        }
    }

  
  # compute the individual row and column heights
  rowheight<-strheight("W",cex=cex) * (1 + rmar)
  colwidth<- apply( object, 2, function(XX) max(strwidth(XX, cex=cex)) ) + 
               strwidth("W")*cmar

  
  width  <- sum(colwidth)
  height <- rowheight*nrow(object)

  # setup x alignment
  if(halign=="left")
    xpos <- 0
  else if(halign=="center")
    xpos <- 0 + (1-width)/2
  else #if(halign=="right")
    xpos <- 0 + (1-width)

  # setup y alignment
  if(valign=="top")
    ypos <- 1
  else if (valign=="center") 
    ypos <- 1 - (1-height)/2
  else #if (valign=="bottom") 
    ypos <- 0 + height

  x <- xpos
  y <- ypos
  
  # iterate across elements, plotting them
  xpos<-x
  for(i in 1:ncol(object)) {
    xpos <- xpos + colwidth[i]
    for(j in 1:nrow(object)) {
      ypos<-y-(j-1)*rowheight
      if( (show.rownames && i==1) || (show.colnames && j==1) )
        text(xpos, ypos, object[j,i], adj=c(hadj,vadj), cex=cex, font=2, ... )
      else
        text(xpos, ypos, object[j,i], adj=c(hadj,vadj), cex=cex, font=1, ... )
    }
  }

  par(opar)
}

textplot.character <- function (object, 
                                halign = c("center", "left", "right"), 
                                valign = c("center", "top", "bottom"), 
                                cex, fixed.width=TRUE,
                                cspace=1,
                                lspace=1,
                                mar=c(0,0,3,0)+0.1,
                                ...)
  {
    # when fixed.withd=TRUE, we simulate a fixed width font by
    # explicitly placing characters at appropriate x,y positions


    object <- paste(object,collapse="\n",sep="")
  
    halign = match.arg(halign)
    valign = match.arg(valign)
    plot.new()

    opar <- par()[c("mar","xpd","cex")]
    on.exit( par(opar) )

    par(mar=mar,xpd=FALSE )
    plot.window(xlim = c(0, 1), ylim = c(0, 1), log = "", asp = NA)

    slist   <- lapply(object, function(x) unlist(strsplit(x,'\n')))[[1]]

    if(length(slist)>1)
      slist   <- lapply(slist, function(x) unlist(strsplit(x,'')))
    else
      slist   <- lapply(slist, function(x) strsplit(x,''))

    slen    <- sapply(slist, length)
    slines  <- length(slist)

    if (missing(cex))
      {
        lastloop <- FALSE
        cex <- 1
      }
    else
      lastloop <- TRUE

    
    for (i in 1:20)
      {
        oldcex <- cex
        
        cwidth  <- max(sapply(unlist(slist), strwidth,  cex=cex)) * cspace
        cheight <- max(sapply(unlist(slist), strheight, cex=cex)) *
                   (lspace * 1.5)

        if(fixed.width)
          {
            width  <- max(slen) * cwidth
            height <- slines * cheight
          }
        else
          {
            width <- strwidth(object, cex=cex)
            height <- strheight(object, cex=cex)
          }

        if(lastloop) break

        if (abs(oldcex - cex) < 0.001)
          {
            lastloop <- TRUE
            cex <- oldcex * (1/max(width, height))
          }

      }

    if (halign == "left")
        xpos <- 0
    else if (halign == "center")
        xpos <- 0 + (1 - width)/2
    else xpos <- 0 + (1 - width)

    if (valign == "top")
        ypos <- 1
    else if (valign == "center")
        ypos <- 1 - (1 - height)/2
    else ypos <- 1 - (1 - height)

    if(fixed.width)
      {
        
        for(line in 1:slines)
          {
            cpos <- ((1:length(slist[[line]]))-1) * cwidth
            if(length(slist[[line]]))
              text(x = xpos + cpos, y = ypos - (line-1)*cheight,
                   labels=slist[[line]], adj = c(0, 1), cex = cex, ...)
          }
      }
    else
      {
        text(x=xpos, y=ypos, labels=object, adj=c(0,1),
             cex=cex, ...)
      }

    par(opar)
    invisible(cex)
}
# $Id: trim.R,v 1.1 2003/05/20 13:16:47 warnes Exp $
#
# $Log: trim.R,v $
# Revision 1.1  2003/05/20 13:16:47  warnes
#
# - Added function trim() and assocated docs.
#
#

trim <- function(s)
  {
    s <- sub("^ +","",s)
    s <- sub(" +$","",s)
    s
  }
# $Id: wapply.R,v 1.9 2003/11/10 22:11:13 warnes Exp $
#
# $Log: wapply.R,v $
# Revision 1.9  2003/11/10 22:11:13  warnes
#
# - Add files contributed by Arni Magnusson
#   <arnima@u.washington.edu>. As well as some of my own.
#
# Revision 1.8  2003/01/20 17:13:04  warnes
#
# - Updated wapply.R to allow specification of evaluation points when
#   method is 'width' or 'range' using the 'pts' argument.
# - Updated wapply.Rd to add 'pts' argument
# - Fixed typos, spelling errors, gramatical errors and lack of clarity
#   in wapply.Rd help text.
#
# Revision 1.7  2002/08/01 19:37:14  warnes
#
# - Corrected documentation mismatch for ci, ci.default.
#
# - Replaced all occurences of '_' for assignment with '<-'.
#
# - Replaced all occurences of 'T' or 'F' for 'TRUE' and 'FALSE' with
#   the spelled out version.
#
# - Updaded version number and date.
#
# Revision 1.6  2002/04/09 00:51:31  warneg
#
# Checkin for version 0.5.3
#
# Revision 1.5  2002/02/16 17:58:59  warneg
#
# - Fixed Bug: When method=="range", the absolute range of x was being
#   used to compute the relative width instead of the (correct) relative
#   range.
#
# Revision 1.4  2002/02/16 17:23:33  warneg
#
# - Corrected problem removing missing values: The missing values of $x
#   and $y were being removed indepdendently, leaving an uneven number
#   of points in the result.
#
# Revision 1.3  2001/12/05 19:29:50  warneg
#
# - Added a better default for "width" when method="nobs".  For this case,
#   width=max(5, length(x)/10).
#
# - Allow omission of x values which result in missing y values via
#   'drop.na' parameter.
#
# Revision 1.2  2001/08/31 23:45:45  warneg
# Used wrong character in header (% instead of #).  Fixed.
#
# Revision 1.1  2001/08/25 05:48:54  warneg
# Initial checkin.
#
#
"wapply" <- function( x, y, fun=mean, method="range",
                    width, n=50, drop.na=TRUE, pts, ...)
{
  method <- match.arg(method, c("width","range","nobs","fraction"))
  if(missing(width))
    if( method=="nobs" ) width <- max(5, length(x)/10 )
  else
    width <- 1/10

  if(method == "width" || method == "range" )
    {
      if(method=="range")
        width <- width * diff(range(x))

      if(missing(pts))
        pts <- seq(min(x),max(x),length.out=n)
      
      result <- sapply( pts, function(pts,y,width,fun,XX,...)
                      {
                        low <- min((pts-width/2),max(XX))
                        high <- max((pts+width/2), min(XX))
                        return (fun(y[(XX>= low) & (XX<=high)],...))
                      },
                      y=y,
                      width=width,
                      fun=fun,
                      XX = x,
                      ...)
      if(drop.na)
        {
          missing <- is.na(pts) & is.na(result)
          pts <- pts[!missing]
          result <- result[!missing]
        }
      
      return(list(x=pts,y=result))
    }
  else # method=="nobs" || method=="fraction"
    {
      if( method=="fraction")
        width <- floor(length(x) * width)

      ord <- order(x)
      x  <- x[ord]
      y  <- y[ord]

      n  <- length(x)
      center  <- 1:n
      below <- sapply(center - width/2, function(XX) max(1,XX) )
      above <- sapply(center + width/2, function(XX) min(n,XX) )

      retval  <- list()
      retval$x  <- x
      retval$y  <- apply(cbind(below,above), 1,
                         function(x) fun(y[x[1]:x[2]],...) )
                         
      if(drop.na)
        {
          missing <- is.na(retval$x) | is.na(retval$y)
          retval$x <- retval$x[!missing]
          retval$y <- retval$y[!missing]
        }


      return(retval)
    }
      
}


.First.lib <- function(libname, pkgname)
  {
    if(is.R())
      {
        library(stats)
        library(MASS)   
      }
  }
