"arrayP" <-
function(data, nperm=100, method=c("fc","pv"), verbose=TRUE, global=TRUE, design=NULL, contrasts=1)
{
  # Calculate original DE statistic
    method <- match.arg(method)
    if(method=="fc") t <- abs(contrasts.fit(lmFit(data, design=design), contrasts=contrasts)$coefficients)
    else if(method=="pv") t <- 1 - eBayes(contrasts.fit(lmFit(data, design=design), contrasts=contrasts))$p.value

  # Generate random DE statistics
    t.rand <- matrix(nrow=nrow(data), ncol=nperm)
    for(i in 1:nperm) {
      data.rand <- data[,sample(ncol(data))]
      if(method=="fc") t.rand[,i] <- abs(contrasts.fit(lmFit(data.rand, design=design), contrasts=contrasts)$coefficients)
      else if(method=="pv") t.rand[,i] <- 1 - eBayes(contrasts.fit(lmFit(data.rand, design=design), contrasts=contrasts))$p.value
      if(verbose && !(i %% 10)) cat(paste(i, " ", sep=""))
    }
    if(verbose) cat("\n")

  # Compute p-values
    t.pvalues <- rep(NA,nrow(data))
    for(j in 1:nrow(data)) {
      if(global) t.pvalues[j] <- min(1, sum(t.rand >= t[j])/nperm)
      else t.pvalues[j] <- sum(t.rand[j,] >= t[j])/nperm
      if(verbose && !(j %% 100)) cat(paste(j, " ", sep=""))
    }
    if(verbose) cat("\n")

    return(t.pvalues)
}