pal <- function(n, alpha=1)
{
  if (n==2) {
    val <- hcl(seq(15,375,len=4), l=60, c=150, alpha=alpha)[c(1,3)]
  } else val <- hcl(seq(15,375,len=n+1), l=60, c=150, alpha=alpha)[1:n]
  val
}
toplegend <- function(...) {
  if (par("oma")[3]==0) {
    x <- mean(par("usr")[1:2])
    yy <- transform.coord(par("usr")[3:4], par("plt")[3:4])
    y  <- mean(c(yy[2],par("usr")[4]))
    legend(x, y, xpd=NA, bty="n", xjust=0.5, yjust=0.5, ...)    
  } else {
    g <- par("mfrow")
    xx <- transform.coord(par("usr")[1:2], par("plt")[1:2])
    yy <- transform.coord(par("usr")[3:4], par("plt")[3:4])
    xxx <- transform.coord(xx, c(g[2]-1,g[2])/g[2])
    yyy <- transform.coord(yy, c(g[1]-1,g[1])/g[1])
    yyyy <- transform.coord(yyy, par("omd")[3:4])
    legend(mean(xxx), mean(c(yyy[2],yyyy[2])), xpd=NA, bty="n", xjust=0.5, yjust=0.5, ...)
  }
}
rightlegend <- function(...) {
  if (par("oma")[3]==0) {
    y <- mean(par("usr")[3:4])
    xx <- transform.coord(par("usr")[1:2], par("plt")[1:2])
    x <- mean(c(xx[2],par("usr")[2]))
    legend(x, y, xpd=NA, bty="n", xjust=0.5, yjust=0.5, ...)    
  } else {
    g <- par("mfrow")
    xx <- transform.coord(par("usr")[1:2], par("plt")[1:2])
    yy <- transform.coord(par("usr")[3:4], par("plt")[3:4])
    xxx <- transform.coord(xx, c(g[2]-1,g[2])/g[2])
    yyy <- transform.coord(yy, c(g[1]-1,g[1])/g[1])
    yyyy <- transform.coord(yyy, par("omd")[3:4])
    legend(mean(xxx), mean(c(yyy[2],yyyy[2])), xpd=NA, bty="n", xjust=0.5, yjust=0.5, ...)
  }  
}
transform.coord <- function(x,p) {
  ba <- (x[2]-x[1])/(p[2]-p[1])
  a <- x[1]-p[1]*ba
  b <- a + ba
  c(a,b)
}
dnplot <- function(M, labs, lwd=3, lty=1, col=pal(p), ylab="Density", xlab="", las=1, dens.par=NULL, ...) {
  if (!is.matrix(M)) M <- matrix(M, ncol=1)
  p <- ncol(M)
  X <- Y <- NULL
  dots <- list(...)
  if ("xlim" %in% names(dots)) {
    dens.par$from <- dots$xlim[1]
    dens.par$to <- dots$xlim[2]
  }
  mode <- numeric(p)
  for (i in 1:p) {
    dens.args <- list(x=M[,i])
    if (length(dens.par)) dens.args[names(dens.par)] <- dens.par
    d <- do.call("density", dens.args)
    X <- cbind(X, d$x)
    Y <- cbind(Y, d$y)
    mode[i] <- d$x[which.max(d$y)]
  }
  matplot(X, Y, lty=lty, lwd=lwd, type="l", col=col, ylab=ylab, xlab=xlab, las=1, ...)
  if (!missing(labs)) toplegend(legend=labs, col=col, lwd=lwd, ncol=min(p, 4))
  invisible(mode)
}
displayProgressBar <- function(i,nI,increments=30)
  {
    if (nI < increments) increments <- nI
    if (i==1)
      {
        cat("Progress:\n")
        cat("|",rep(" ",increments),"|\n ",sep="")
      }
    display.iterates <- ceiling((1:increments) * (nI/increments))
    if (is.element(i,display.iterates)) cat("*")
    if (Sys.info()["sysname"]=="Windows") flush.console()
    if (i==nI) cat("\n")
  }
CIplot <- function(obj, labels=rownames(B), pxlim, xlim, ylim, sub, diff=(ncol(B)==4), n.ticks=6, mar=c(5,nn/3+.5,2,4), axis=TRUE, trans, p.label=FALSE, ...)
{
  B <- obj
  nn <- max(nchar(rownames(B)))
  par(mar=mar)
  n <- nrow(B)
  if (!missing(trans)) B[,1:3] <- trans(B[,1:3])
  
  if (missing(pxlim)) {
    pxlim <- if (missing(xlim)) pretty(range(B[,2:3], na.rm=TRUE),n=n.ticks-1) else pretty(xlim, n=n.ticks-1)
  }
  if (missing(ylim)) ylim <- c(0.5,n+0.5)
  
  plot(B[n:1,1],1:n,xlim = range(pxlim), ylim=ylim, ylab="", axes=FALSE,pch=19,...)
  for (i in 1:n) {
    lines(c(B[i,2:3]),c(n-i+1,n-i+1),lwd=2)
    if (diff) {
      p <- formatP(B[,4], label=p.label)
      p[is.na(B[,4])] <- ""
      mtext(at=n-i+1,p[i],line=1,side=4,las=1,cex=0.7,adj=0)
    }
  }
  if (axis) axis(1, pxlim)
  if (diff) {
    null <- 0
    if (!missing(trans)) null <- trans(0)
    abline(v=null,col="gray")
  }
  if (!missing(sub)) mtext(sub,3,0,cex=0.8)
  
  ind <- !is.na(B[,1])
  text(x=par("usr")[1], adj=1, y=(n:1)[ind], labels=labels[ind], xpd=TRUE, cex=.8)
  if (sum(!ind) > 0) {
    a <- diff(par("usr")[1:2])/diff(par("plt")[1:2])
    b <- par("usr")[1] - a*par("plt")[1]
    text(x=b+a*.01, adj=0, y=(n:1)[!ind], labels=labels[!ind], xpd=TRUE, cex=.8)    
  }
  invisible(B)
}
drawBand <- function(x, l, u, border=NA, col="gray90", ...) {
  polygon(c(x, rev(x)), c(l, rev(u)), border=border, col=col, ...)
}
psm <- function(x, labels, level=.95) {
  if (is.null(dim(x))) x <- as.matrix(x)
  alpha <- 1-level
  probs <- c(.5, alpha/2, 1-alpha/2)
  P <- t(apply(x, 2, quantile, probs=probs))
  if (!missing(labels)) rownames(P) <- labels
  P
}
ilogit <- binomial()$linkinv
rdirichlet <- function(n, alpha) {
  l <- length(alpha)
  x <- matrix(rgamma(l * n, alpha), ncol = l, byrow = TRUE)
  sm <- x %*% rep(1, l)
  x/as.vector(sm)
}