`cov2modules` <-
function(cov.mat, is.DE=rep(FALSE,nrow(cov.mat)), bin.cutoff, clust.method=NULL, mod.cutoff=NULL, plot=FALSE, plot.half=FALSE, new=TRUE)
{
  # Sanity check
  if(nrow(cov.mat) != ncol(cov.mat)) stop("Invalid covariance matrix!")
  N <- nrow(cov.mat)
  
  # Binarize
  require(biclust, quietly=TRUE)
  cov.mat <- binarize(abs(cov.mat), bin.cutoff)
  diag(cov.mat) <- 1

  # Hclust
  if(!is.null(clust.method)) {
    clust.order <- hclust(d=as.dist(1-cov.mat), method=clust.method)$order
    cov.mat <- cov.mat[clust.order, clust.order]
    is.DE <- is.DE[clust.order]
  }

  # Pruning
  if(!is.null(mod.cutoff)) {
    which.module <- list()
    for(i in 1:(N-mod.cutoff+1)) {
      if(sum(cov.mat[i:(i+mod.cutoff-1), i:(i+mod.cutoff-1)]) == mod.cutoff^2) {
        
        # Module extension
        mod.ext.cutoff <- mod.cutoff
        while((i+mod.ext.cutoff) <= N) {
          if(sum(cov.mat[i:(i+mod.ext.cutoff), i:(i+mod.ext.cutoff)]) == (mod.ext.cutoff+1)^2) mod.ext.cutoff <- mod.ext.cutoff + 1
          else break
        }

        # Accumulate modules
        which.module <- c(which.module, list(c(i,i+mod.ext.cutoff-1)))
      }
    }

    # Modify matrix
    cov.mat[] <- 0
    for(mod in which.module) {
      cov.mat[mod[1]:mod[2], mod[1]:mod[2]] <- 1
    }
    cat(paste(sum(base:::rowSums(cov.mat) > 0), "\n", sep=""))
  }

  # DE gene labeling
  if(any(is.DE)) {
    cov.mat[is.DE, is.DE][cov.mat[is.DE, is.DE] == 1] <- 2
    col <- c("white","black","red")
  }
  else col <- c("white","black")
    

  # Plot
  if(plot) {
    cov.mat.plot <- cov.mat
    if(plot.half) cov.mat.plot[lower.tri(cov.mat.plot)] <- NA
    if(new) x11()
    image(x=1:N, y=1:N, z=cov.mat.plot, ylim=c(N, 1), col=col)
  }

  return(cov.mat)
}