DoubletDeconSeurat: Seurat wrappers with DoubletDecon

#' DDHeatmap
#'
#' This function is a modified version of heatmap.2 from the gplots package (https://cran.r-project.org/package=gplots), corrected by ChristophH (https://github.com/ChristophH/gplots/blob/master/R/heatmap.2.R) to correct the issue of dendrogram generation even with row and column clustering off and dendrogram set to "none".
#'
#' @keywords heatmap heatmap.2
#'
#' @import gplots
#'
#' @export
#'
DDheatmap <- function (
  x,

  ## dendrogram control
  Rowv = TRUE,
  Colv = if (symm) 'Rowv' else TRUE,
  distfun = dist,
  hclustfun = hclust,
  dendrogram = c('both', 'row', 'column', 'none'),
  reorderfun = function(d, w) reorder(d, w),
  symm = FALSE,

  ## data scaling
  scale = c('none', 'row', 'column'),
  na.rm = TRUE,

  ## image plot
  revC = identical(Colv, 'Rowv'),
  add.expr,

  ## mapping data to colors
  breaks,
  symbreaks = any(x < 0, na.rm = TRUE) || scale != 'none',

  ## colors
  col = 'heat.colors',

  ## block sepration
  colsep,
  rowsep,
  sepcolor = 'white',
  sepwidth = c(0.05,0.05),

  ## cell labeling
  cellnote,
  notecex = 1.0,
  notecol = 'cyan',
  na.color = par('bg'),

  ## 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,
  srtRow = NULL,
  srtCol = NULL,
  adjRow = c(0,NA),
  adjCol = c(NA,0),
  offsetRow = 0.5,
  offsetCol = 0.5,
  colRow = NULL,
  colCol = NULL,

  ## color key + density info
  key = TRUE,
  keysize = 1.5,
  density.info = c('histogram', 'density', 'none'),
  denscol = tracecol,
  symkey = any(x < 0, na.rm = TRUE) || symbreaks,
  densadj = 0.25,
  key.title = NULL,
  key.xlab = NULL,
  key.ylab = NULL,
  key.xtickfun = NULL,
  key.ytickfun = NULL,
  key.par = list(),

  ## plot labels
  main = NULL,
  xlab = NULL,
  ylab = NULL,

  ## plot layout
  lmat = NULL,
  lhei = NULL,
  lwid = NULL,

  ## extras
  extrafun = NULL,
  ...
) {
  scale01 <- function(x, low = min(x), high = max(x)) {
    x <- (x - low)/(high - low)
    x
  }

  retval <- list()

  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 (length(x = col) == 1 && is.character(col))
    col <- get(col, mode = 'function')

  if (!missing(breaks) && any(duplicated(breaks)))
    stop('breaks may not contian duplicate values')

  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.")

  if (is.null(Rowv) || is.na(Rowv))
    Rowv <- FALSE
  if (is.null(Colv) || is.na(Colv))
    Colv <- FALSE
  else if (all(Colv == 'Rowv'))
    Colv <- Rowv

  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")

  if (missing(cellnote))
    cellnote <- matrix("", ncol = ncol(x), nrow = nrow(x))

  if (!inherits(Rowv, 'dendrogram')) {
    ## Check if Rowv and dendrogram arguments are consistent
    if (((is.logical(Rowv) && !isTRUE(Rowv)) || (is.null(Rowv))) && (dendrogram %in% c('both', 'row'))) {
      warning("Discrepancy: Rowv is FALSE, while dendrogram is `", dendrogram, "'. Omitting row dendogram.")
      if (dendrogram == 'both')
        dendrogram <- 'column'
      else
        dendrogram <- 'none'
    }
  }

  if (!inherits(Colv, 'dendrogram')) {
    ## Check if Colv and dendrogram arguments are consistent
    if (((is.logical(Colv) && !isTRUE(Colv)) || (is.null(Colv))) && (dendrogram %in% c('both', 'column'))) {
      warning("Discrepancy: Colv is FALSE, while dendrogram is `", dendrogram, "'. Omitting column dendogram.")
      if (dendrogram == 'both')
        dendrogram <- 'row'
      else
        dendrogram <- 'none'
    }
  }

  ## 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") )
  ##  { ## dendrogram option is used *only* for display purposes
  if (inherits(Rowv, 'dendrogram')) {
    ddr <- Rowv ## use Rowv 'as-is', when it is dendrogram
    rowInd <- order.dendrogram(ddr)
    if (length(x = rowInd) > nr || any(rowInd < 1 | rowInd > nr))
      stop("Rowv dendrogram doesn't match size of x")
    if (length(x = rowInd) < nr)
      nr <- length(x = rowInd)
  } else if (is.integer(Rowv)) {
    ## Compute dendrogram and do reordering based on given vector
    distr <- distfun(x)
    hcr <- hclustfun(distr)
    ddr <- as.dendrogram(hcr)
    ddr <- reorderfun(ddr, Rowv)

    rowInd <- order.dendrogram(ddr)
    if (nr != length(x = rowInd))
      stop("row dendrogram ordering gave index of wrong length")
  } else if (isTRUE(Rowv)) { ## If TRUE, compute dendrogram and do reordering based on rowMeans
    Rowv <- rowMeans(x, na.rm = na.rm)
    distr <- distfun(x)
    hcr <- hclustfun(distr)
    ddr <- as.dendrogram(hcr)
    ddr <- reorderfun(ddr, Rowv)

    rowInd <- order.dendrogram(ddr)
    if (nr != length(x = rowInd))
      stop("row dendrogram ordering gave index of wrong length")
  } else {
    rowInd <- nr:1
  }

  if (inherits(Colv, 'dendrogram')) {
    ddc <- Colv ## use Colv 'as-is', when it is dendrogram
    colInd <- order.dendrogram(ddc)
    if (length(x = colInd) > nc || any(colInd < 1 | colInd > nc))
      stop("Colv dendrogram doesn't match size of x")
    if (length(colInd) < nc)
      nc <- length(x = colInd)
  } else if (identical(Colv, 'Rowv')) {
    if (nr != nc)
      stop('Colv = "Rowv" but nrow(x) != ncol(x)')
    if (exists('ddr')) {
      ddc <- ddr
      colInd <- order.dendrogram(ddc)
    } else
      colInd <- rowInd
  } else if (is.integer(Colv)) {## Compute dendrogram and do reordering based on given vector
    distc <- distfun(if (symm) x else t(x))
    hcc <- hclustfun(distc)
    ddc <- as.dendrogram(hcc)
    ddc <- reorderfun(ddc, Colv)

    colInd <- order.dendrogram(ddc)
    if (nc != length(x = colInd))
      stop("column dendrogram ordering gave index of wrong length")
  } else if (isTRUE(Colv)) {## If TRUE, compute dendrogram and do reordering based on rowMeans
    Colv <- colMeans(x, na.rm = na.rm)
    distc <- distfun(if (symm) x else t(x))
    hcc <- hclustfun(distc)
    ddc <- as.dendrogram(hcc)
    ddc <- reorderfun(ddc, Colv)

    colInd <- order.dendrogram(ddc)
    if(nc != length(x = colInd))
      stop("column dendrogram ordering gave index of wrong length")
  } else {
    colInd <- 1:nc
  }

  retval$rowInd <- rowInd
  retval$colInd <- colInd
  retval$call <- match.call()


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

  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 (!is.null(colRow))
    colRow <- colRow[rowInd]

  if (!is.null(colCol))
    colCol <- colCol[colInd]

  if (scale == 'row') {
    retval$rowMeans <- rm <- rowMeans(x, na.rm = na.rm)
    x <- sweep(x, 1, rm)
    retval$rowSDs <-  sx <- apply(x, 1, sd, na.rm = na.rm)
    x <- sweep(x, 1, sx, '/')
  } else if (scale == 'column') {
    retval$colMeans <- rm <- colMeans(x, na.rm = na.rm)
    x <- sweep(x, 2, rm)
    retval$colSDs <-  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) ||  is.function(col))
      breaks <- 16
    else
      breaks <- length(x = col) + 1
  }

  if (length(x = breaks) == 1) {
    if(!symbreaks)
      breaks <- seq(min(x, na.rm = na.rm), max(x, na.rm = na.rm), length = breaks)
    else {
      extreme <- max(abs(x), na.rm = TRUE)
      breaks <- seq( -extreme, extreme, length = breaks )
    }
  }

  nbr <- length(x = breaks)
  ncol <- length(x = breaks)-1

  if (class(col) == 'function')
    col <- col(ncol)

  min.breaks <- min(breaks)
  max.breaks <- max(breaks)

  x[x < min.breaks] <- min.breaks
  x[x > max.breaks] <- max.breaks

  ## Calculate the plot layout
  if (missing(lhei) || is.null(lhei))
    lhei <- c(keysize, 4)

  if (missing(lwid) || is.null(lwid))
    lwid <- c(keysize, 4)

  if (missing(lmat) || is.null(lmat)) {
    lmat <- rbind(4:3, 2:1)

    if (!missing(ColSideColors)) { ## add middle row to layout
      if (!is.character(ColSideColors) || length(x = 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(x = lmat) - 1), 1), lmat[, 2] + 1)
      lwid <- c(lwid[1], 0.2, lwid[2])
    }

    lmat[is.na(lmat)] <- 0
  }

  if (length(lhei) != nrow(lmat))
    stop("lhei must have length = nrow(lmat) = ", nrow(x = lmat))

  if (length(lwid) != ncol(lmat))
    stop("lwid must have length = ncol(lmat) =", ncol(x = lmat))

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

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

  plot.index <- 1

  ## 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)
    plot.index <- plot.index + 1
  }
  if (!missing(ColSideColors)) {
    par(mar = c(0.5, 0, 0,margins[2]))
    image(cbind(1:nc), col = ColSideColors[colInd], axes = FALSE)
    plot.index <- plot.index + 1
  }
  ## draw the main carpet
  par(mar = c(margins[1], 0, 0, margins[2]))
  #if(scale != "none" || !symm)
  #  {
  x <- t(x)
  cellnote <- t(cellnote)
  #  }
  if (revC) { ## x columns reversed
    iy <- nr:1
    if(exists('ddr'))
      ddr <- rev(ddr)
    x <- x[, iy]
    cellnote <- cellnote[, iy]
  }
  else iy <- 1:nr

  ## display the main carpet
  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,
        ...)
  retval$carpet <- x
  if (exists('ddr'))
    retval$rowDendrogram <- ddr
  if (exists('ddc'))
    retval$colDendrogram <- ddc
  retval$breaks <- breaks
  retval$col <- col

  ## fill 'na' positions with na.color
  if (!gtools::invalid(na.color) & any(is.na(x))) {
    mmat <- ifelse(is.na(x), 1, NA)
    image(1:nc, 1:nr, mmat, axes = FALSE, xlab = '', ylab = '', col = na.color, add = TRUE)
  }

  ## add column labels
  if (is.null(srtCol) && is.null(colCol))
    axis(1,
         1:nc,
         labels= labCol,
         las= 2,
         line = -0.5 + offsetCol,
         tick = 0,
         cex.axis = cexCol,
         hadj = adjCol[1],
         padj = adjCol[2]
    )
  else
  {
    if (is.null(srtCol) || is.numeric(srtCol)) {
      if (missing(adjCol) || is.null(adjCol))
        adjCol = c(1, NA)

      if(is.null(srtCol))
        srtCol <- 90

      xpd.orig <- par('xpd')
      par(xpd = NA)
      xpos <- axis(1, 1:nc, labels = rep('', nc), las = 2, tick = 0)
      text(x = xpos,
           y = par('usr')[3] - (1.0 + offsetCol) * strheight('M'),
           labels = labCol,
           ##pos=1,
           adj = adjCol,
           cex = cexCol,
           srt = srtCol,
           col = colCol
      )
      par(xpd = xpd.orig)
    }
    else
      warning("Invalid value for srtCol ignored.")
  }

  ## add row labels
  if(is.null(srtRow) && is.null(colRow)) {
    axis(4,
         iy,
         labels = labRow,
         las = 2,
         line = -0.5 + offsetRow,
         tick = 0,
         cex.axis = cexRow,
         hadj = adjRow[1],
         padj = adjRow[2]
    )
  } else {
    if (is.null(srtRow) || is.numeric(srtRow)) {
      xpd.orig <- par('xpd')
      par(xpd = NA)
      ypos <- axis(4, iy, labels = rep('', nr), las = 2, line = -0.5, tick = 0)
      text(x = par('usr')[2] + (1.0 + offsetRow) * strwidth('M'),
           y = ypos,
           labels = labRow,
           adj = adjRow,
           cex = cexRow,
           srt = srtRow,
           col = colRow
      )
      par(xpd = xpd.orig)
    }
    else
      warning("Invalid value for srtRow ignored.")
  }

  ## add row and column headings (xlab, ylab)
  if (!is.null(xlab)) mtext(xlab, side = 1, line = margins[1] - 1.25)
  if (!is.null(ylab)) mtext(ylab, side = 4, line = margins[2] - 1.25)

  ## perform user-specified function
  if (!missing(add.expr))
    eval(substitute(add.expr))

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

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

  ## 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')) {
    retval$vline <- vline
    vline.vals <- scale01(vline, min.scale, max.scale)
    for(i in 1:length(colInd)) {
      if (!is.null(vline)) {
        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(x = xv) - 0.5
      lines(x = xv, y = yv, lwd = 1, col = tracecol, type = 's')
    }
  }

  if (trace %in% c('both', 'row')) {
    retval$hline <- hline
    hline.vals <- scale01(hline, min.scale, max.scale)
    for (i in 1:length(rowInd)) {
      if (!is.null(hline)) {
        abline(h = i - 0.5 + hline.vals, 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')
    }
  }

  if (!missing(cellnote))
    text(x = c(row(cellnote)),
         y = c(col(cellnote)),
         labels = c(cellnote),
         col = notecol,
         cex = notecex)

  plot.index <- plot.index + 1

  ## increment plot.index and then do
  ##   latout_set( lmat, plot.index )
  ## to set to the correct plot region, instead of
  ## relying on plot.new().

  ## the two dendrograms :
  par(mar = c(margins[1], 0, 0, 0))
  if (dendrogram %in% c('both', 'row')) {
    flag <- try(
      plot.dendrogram(ddr, horiz = TRUE, axes = FALSE, yaxs = 'i', leaflab = 'none')
    )
    if ('try-error' %in% class(flag)) {
      cond <- attr(flag, 'condition')
      if (!is.null(cond) && conditionMessage(cond) == 'evaluation nested too deeply: infinite recursion / options(expressions=)?')
        stop('Row dendrogram too deeply nested, recursion limit exceeded.  Try increasing option("expressions"=...).')
    }
  } else
    plot.new()

  par(mar = c(0, 0, if (!is.null(main)) 5 else 0, margins[2]))

  if (dendrogram %in% c('both', 'column')) {
    flag <- try(
      plot.dendrogram(ddc, axes = FALSE, xaxs = 'i', leaflab = 'none')
    )
    if ('try-error' %in% class(flag)) {
      cond <- attr(flag, 'condition')
      if(!is.null(cond) && conditionMessage(cond) == 'evaluation nested too deeply: infinite recursion / options(expressions=)?')
        stop('Column dendrogram too deeply nested, recursion limit exceeded.  Try increasing option("expressions"=...).')
    }
  } else
    plot.new()

  ## title
  if (!is.null(main)) title(main, cex.main = 1.5*op[['cex.main']])

  ## Add the color-key
  if (key) {
    mar <- c(5, 4, 2, 1)
    if (!is.null(key.xlab) && is.na(key.xlab))
      mar[1] <- 2
    if (!is.null(key.ylab) && is.na(key.ylab))
      mar[2] <- 2
    if (!is.null(key.title) && is.na(key.title))
      mar[3] <- 1
    par(mar = mar, cex = 0.75, mgp = c(2, 1, 0))
    if (length(key.par) > 0)
      do.call(par, key.par)
    tmpbreaks <- breaks

    if (symkey) {
      max.raw <- max(abs(c(x, breaks)), na.rm=TRUE)
      min.raw <- -max.raw
      tmpbreaks[1] <- -max(abs(x), na.rm = TRUE)
      tmpbreaks[length(x = tmpbreaks)] <- max(abs(x), na.rm = TRUE)
    } else {
      min.raw <- min.breaks
      max.raw <- max.breaks
    }

    z <- seq(min.raw, max.raw, by = min(diff(breaks)/100))
    image(z = matrix(z, ncol = 1),
          col = col, breaks = tmpbreaks,
          xaxt = 'n', yaxt = 'n')

    par(usr = c(0,1,0,1))
    if (is.null(key.xtickfun)) {
      lv <- pretty(breaks)
      xv <- scale01(as.numeric(lv), min.raw, max.raw)
      xargs <- list(at=xv, labels=lv)
    } else {
      xargs <- key.xtickfun()
    }
    xargs$side <- 1
    do.call(axis, xargs)
    if (is.null(key.xlab)) {
      if (scale == 'row')
        key.xlab <- 'Row Z-Score'
      else if (scale == 'column')
        key.xlab <- 'Column Z-Score'
      else
        key.xlab <- 'Value'
    }
    if (!is.na(key.xlab)) {
      mtext(side = 1, key.xlab, line = par('mgp')[1], padj = 0.5, cex = par('cex') * par('cex.lab'))
    }

    if (density.info == 'density') {
      dens <- density(x, adjust = densadj, na.rm = TRUE, from = min.scale, to = max.scale)
      omit <- dens$x < min(breaks) | dens$x > max(breaks)
      dens$x <- dens$x[!omit]
      dens$y <- dens$y[!omit]
      dens$x <- scale01(dens$x, min.raw, max.raw)
      lines(dens$x, dens$y/max(dens$y)*0.95, col = denscol, lwd = 1)
      if (is.null(key.ytickfun)) {
        yargs <- list(at = pretty(dens$y)/max(dens$y)*0.95, labels = pretty(dens$y))
      } else {
        yargs <- key.ytickfun()
      }
      yargs$side <- 2
      do.call(axis, yargs)
      if (is.null(key.title))
        key.title <- 'Color Key\nand Density Plot'
      if (!is.na(key.title))
        title(key.title)
      par(cex = 0.5)
      if (is.null(key.ylab))
        key.ylab <- 'Density'
      if (!is.na(key.ylab))
        mtext(side = 2,key.ylab, line = par('mgp')[1], padj = 0.5, cex = par('cex')*par('cex.lab'))
    } 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)
      if (is.null(key.ytickfun)) {
        yargs <- list(at = pretty(hy)/max(hy)*0.95, labels = pretty(hy))
      } else {
        yargs <- key.ytickfun()
      }
      yargs$side <- 2
      do.call(axis, yargs)
      if (is.null(key.title))
        key.title <- 'Color Key\nand Histogram'
      if (!is.na(key.title))
        title(key.title)
      par(cex = 0.5)
      if (is.null(key.ylab))
        key.ylab <- 'Count'
      if (!is.na(key.ylab))
        mtext(side = 2, key.ylab, line = par('mgp')[1], padj = 0.5, cex = par('cex')*par('cex.lab'))
    }
    else
      if (is.null(key.title))
        title("Color Key")

    if(trace %in% c('both', 'column')) {
      vline.vals <- scale01(vline, min.raw, max.raw)
      if (!is.null(vline)) {
        abline(v = vline.vals, col = linecol, lty = 2)
      }
    }

    if (trace %in% c('both', 'row')) {
      hline.vals <- scale01(hline, min.raw, max.raw)
      if(!is.null(hline)) {
        abline(v = hline.vals, col = linecol, lty = 2)
      }
    }
  } else {
    par(mar = c(0, 0, 0, 0))
    plot.new()
  }
  ## Create a table showing how colors match to (transformed) data ranges
  retval$colorTable <- data.frame(
    low = retval$breaks[-length(retval$breaks)],
    high = retval$breaks[-1],
    color = retval$col
  )

  # Store layout information, suggested by Jenny Drnevich
  retval$layout <- list(lmat = lmat,
                        lhei = lhei,
                        lwid = lwid
  )

  ## If user has provided an extra function, call it.
  if (!is.null(extrafun))
    extrafun()

  invisible(retval)
}