R/utils.R
d1da987e
 
d3ce90c5
 # environment that contains global variables
d1da987e
 INDEX_ENV = new.env()
 
0623db7b
 INDEX_ENV$I_FIGURE = 0
d1da987e
 INDEX_ENV$I_HEATMAP = 0
 INDEX_ENV$I_ANNOTATION = 0
 INDEX_ENV$I_ROW_ANNOTATION = 0
3bb5581b
 INDEX_ENV$I_COLOR_MAPPING = 0
d1da987e
 
0623db7b
 get_figure_index = function() {
     INDEX_ENV$I_FIGURE
 }
 
 increase_figure_index = function() {
     INDEX_ENV$I_FIGURE = INDEX_ENV$I_FIGURE + 1
 }
 
d1da987e
 get_heatmap_index = function() {
 	INDEX_ENV$I_HEATMAP
 }
 
 increase_heatmap_index = function() {
 	INDEX_ENV$I_HEATMAP = INDEX_ENV$I_HEATMAP + 1
 }
 
 get_annotation_index = function() {
 	INDEX_ENV$I_ANNOTATION
 }
 
 increase_annotation_index = function() {
 	INDEX_ENV$I_ANNOTATION = INDEX_ENV$I_ANNOTATION + 1
 }
 
 get_row_annotation_index = function() {
 	INDEX_ENV$I_ROW_ANNOTATION
 }
 
 increase_row_annotation_index = function() {
 	INDEX_ENV$I_ROW_ANNOTATION = INDEX_ENV$I_ROW_ANNOTATION + 1
 }
 
3bb5581b
 get_color_mapping_index = function() {
     INDEX_ENV$I_COLOR_MAPPING
 }
 
 increase_color_mapping_index = function() {
     INDEX_ENV$I_COLOR_MAPPING = INDEX_ENV$I_COLOR_MAPPING + 1
 }
d1da987e
 
 # default colors for matrix or annotations
 # this function should be improved later
 default_col = function(x, main_matrix = FALSE) {
 
     if(is.factor(x)) {
         x = as.vector(x)
     }
 
     if(length(unique(x)) == 1) {
         x = as.character(x)
     }
 
     attributes(x) = NULL
 
efbb6324
     x = x[!is.na(x)]
 
d1da987e
     if(is.character(x)) {  # discrete
         levels = unique(x)
6ae76409
         #colors = hsv(runif(length(levels)), 1-runif(1)/2, 1-runif(1)/2)
         colors = rand_color(length(levels), luminosity = sample(c("bright", "light", "dark", "random"), 1))
d1da987e
         names(colors) = levels
         return(colors)
     } else if(is.numeric(x)) {
         if(main_matrix) {
2a88897d
             p = sum(x > 0)/sum(x != 0)
9c4d56c4
             if(p > 0.25 & p < 0.75) {
d90404d6
                 if(ht_opt$verbose) {
                     cat("This matrix has both negative and positive values, use a color mapping symmetric to zero\n")
                 }
                 if(length(unique(x)) >= 100) {
                     q1 = quantile(abs(x), 0.99)
                     col_fun = colorRamp2(c(-q1, 0, q1), c("blue", "#EEEEEE", "red"))
a93897a0
                 } else {
d90404d6
                     q1 = max(abs(x))
                     col_fun = colorRamp2(c(-q1, 0, q1), c("blue", "#EEEEEE", "red"))
a93897a0
                 }
d1da987e
             } else {
d90404d6
                 if(length(unique(x)) >= 100) {
                     q1 = quantile(x, 0.01)
                     q2 = quantile(x, 0.99)
                     if(length(unique(x[x > q1 & x < q2])) == 1) {
                          col_fun = colorRamp2(seq(min(x), max(x), length = 3), c("blue", "#EEEEEE", "red"))
                     } else {
                         col_fun = colorRamp2(seq(q1, q2, length = 3), c("blue", "#EEEEEE", "red"))
                     }
                 } else {
                     col_fun = colorRamp2(seq(min(x), max(x), length = 3), c("blue", "#EEEEEE", "red"))
                 }
d1da987e
             }
         } else {
6ae76409
             #col_fun = colorRamp2(range(min(x), max(x)), c("white", hsv(runif(1), 1, 1)))
             col_fun = colorRamp2(range(min(x), max(x)), c("white", rand_color(1, luminosity = sample(c("bright", "dark"), 1))))
d1da987e
         }
         return(col_fun)
     }
 }
 
 # == title
e27480b9
 # Calculate Pairwise Distance from a Matrix
d1da987e
 #
 # == param
64d651fe
 # -x A matrix or a list. If it is a matrix, the distance is calculated by rows.
e27480b9
 # -pairwise_fun A function which calculates distance between two vectors.
 # -... Pass to `stats::as.dist`.
d1da987e
 #
 # == detail
 # You can construct any type of distance measurements by defining a pair-wise distance function.
b08d7ba4
 # The function is implemented by two nested ``for`` loops, so the efficiency may not be so good.
d1da987e
 #
 # == value
 # A `stats::dist` object.
 #
 # == author
 # Zuguang Gu <z.gu@dkfz.de>
 #
64d651fe
 # == example
 # lt = lapply(1:10, function(i) {
 #     sample(letters, sample(6:10, 1))
 # })
 # dist2(lt, function(x, y) {
 #     length(intersect(x, y))/length(union(x, y))
 # })
 dist2 = function(x, pairwise_fun = function(x, y) sqrt(sum((x - y)^2)), ...) {
d1da987e
 
64d651fe
     if(is.matrix(x)) {
         if(nrow(x) < 2) {
             stop_wrap("`x` should have at least two rows.")
         }
d1da987e
 
64d651fe
         nr = nrow(x)
         mat2 = matrix(NA, nrow = nr, ncol = nr)
         rownames(mat2) = colnames(mat2) = rownames(x)
d1da987e
 
64d651fe
         for(i in 2:nr) {
             for(j in 1:(nr-1)) {
                 mat2[i, j] = pairwise_fun(x[i, ], x[j, ])
             }
         }
d1da987e
 
64d651fe
         as.dist(mat2, ...)
     } else if(is.list(x)) {
         if(length(x) < 2) {
             stop_wrap("`x` should have at least length of 2.")
d1da987e
         }
 
64d651fe
         nr = length(x)
         mat2 = matrix(NA, nrow = nr, ncol = nr)
         rownames(mat2) = colnames(mat2) = names(x)
 
         for(i in 2:nr) {
             for(j in 1:(nr-1)) {
                 mat2[i, j] = pairwise_fun(x[[i]], x[[j]])
             }
         }
         as.dist(mat2, ...)
     } else {
         stop_wrap("`x` can be a matrix or a list.")
     }
d1da987e
 }
 
 
 get_dist = function(matrix, method) {
     if(is.function(method)) {
         nargs = length(as.list(args(method)))
         if(nargs == 2) { # a distance function
             dst = method(matrix)
         } else if(nargs == 3) {
             dst = dist2(matrix, method)
         } else {
bfaaa9dd
             stop_wrap("Since your distance method is a function, it can only accept one or two arguments.")
d1da987e
         }
     } else if(method %in% c("euclidean", "maximum", "manhattan", "canberra", "binary", "minkowski")) {
900f036d
         # if(any(is.na(matrix))) {
         #     dst = get_dist(matrix, function(x, y) {
         #         l = is.na(x) | is.na(y)
         #         x = x[!l]
         #         y = y[!l]
         #         as.vector(dist(rbind(x, y), method = method))
         #     })
         #     warning("NA exists in the matrix, calculating distance by removing NA values.")
         # } else {
acb6c95b
             dst = dist(matrix, method = method)
900f036d
         # }
d1da987e
     } else if(method %in% c("pearson", "spearman", "kendall")) {
acb6c95b
         if(any(is.na(matrix))) {
             dst = get_dist(matrix, function(x, y) {
                     l = is.na(x) | is.na(y)
                     x = x[!l]
                     y = y[!l]
                     1 - cor(x, y, method = method)
                 })
6f036283
             warning_wrap("NA exists in the matrix, calculating distance by removing NA values.")
acb6c95b
         } else {
             dst = switch(method,
                          pearson = as.dist(1 - cor(t(matrix), method = "pearson")),
                          spearman = as.dist(1 - cor(t(matrix), method = "spearman")),
                          kendall = as.dist(1 - cor(t(matrix), method = "kendall")))
         }
d1da987e
     }
     return(dst)
 }
 
12e85497
 get_dend_order = function(x) {
d1da987e
     switch(class(x),
         hclust = x$order,
         dendrogram = order.dendrogram(x))
 }
 
 recycle_gp = function(gp, n = 1) {
     for(i in seq_along(gp)) {
         x = gp[[i]]
         gp[[i]] = c(rep(x, floor(n/length(x))), x[seq_len(n %% length(x))])
     }
     return(gp)
 }
 
 check_gp = function(gp) {
     if(!inherits(gp, "gpar")) {
6f036283
         stop_wrap("Graphic parameters should be specified by `gpar()`.")
d1da987e
     }
     return(gp)
 }
 
 
e27480b9
 # == title
64d651fe
 # Subset a gpar Object
e27480b9
 #
 # == param
 # -gp A `gpar` object.
30d2c5b2
 # -i A vector of indices.
e27480b9
 #
30d2c5b2
 # == value
 # A `grid::gpar` object.
 #
 # == example
 # gp = gpar(col = 1:10, fill = 1)
 # subset_gp(gp, 1:5)
 subset_gp = function(gp, i) {
acb6c95b
     gp = lapply(gp, function(x) {
         if(length(x) == 1) x
30d2c5b2
         else x[i]
acb6c95b
     })
d1da987e
     class(gp) = "gpar"
     return(gp)
 }
4818ea96
 
 
 get_text_just = function(rot, side) {
     rot = rot %% 360
     if(! rot %in% c(0, 90, 270)) {
6f036283
         stop_wrap("Only support horizontal or vertical rotations for text.\n")
4818ea96
     }
     if(side == "left") {
         if(rot == 0) {
             return(c(1, 0.5))
         } else if(rot == 90) {
             return(c(0.5, 0))
         } else if(rot == 270) {
             return(c(0.5, 1))
         }
     } else if(side == "right") {
         if(rot == 0) {
             return(c(0, 0.5))
         } else if(rot == 90) {
             return(c(0.5, 1))
         } else if(rot == 270) {
             return(c(0.5, 0))
         }
     } else if(side == "top") {
         if(rot == 0) {
             return(c(0.5, 0))
         } else if(rot == 90) {
             return(c(0, 0.5))
         } else if(rot == 270) {
             return(c(1, 0.5))
         }
     } else if(side == "bottom") {
         if(rot == 0) {
             return(c(0.5, 1))
         } else if(rot == 90) {
             return(c(1, 0.5))
         } else if(rot == 270) {
             return(c(0, 0.5))
         }
     }
5ab2cb81
 }
 
 c.list = function(lt, ..., list = NULL) {
     if(length(lt) == 0) lt = list()
 
     if(is.null(list)) {
         lt_add = list(...)
 
         n = length(lt)
         for(i in seq_along(lt_add)) {
             lt[[n+i]] = lt_add[[i]]
         }
     } else {
         lt = c(lt, list)
     }
     return(lt)
 }
 
 rep.list = function(x, n) {
     lt = vector("list", n)
     for(i in seq_len(n)) {
         lt[i] = list(x)
     }
     return(lt)
 }
 
1a56796e
 # == title
 # List All Heatmap Components
 #
 # == value
 # A vector of viewport names.
 #
d9834c56
 list_components = function() {
b8dacfcc
     vp = grid.ls(viewports = TRUE, grobs = FALSE, flatten = FALSE, print = FALSE)
933e808c
     vp = unlist(vp)
     attributes(vp) = NULL
     vp = vp[!grepl("^\\d+$", vp)]
     vp = vp[!grepl("GRID.VP", vp)]
a4788e64
     # unique(vp)
     vp
076c1ff2
 }
 
 # == title
e27480b9
 # Maximum Width of Text
076c1ff2
 #
 # == param
30d2c5b2
 # -text A vector of text.
ad35494a
 # -gp Graphic parameters for text.
076c1ff2
 #
793b75de
 # == details
64d651fe
 # It simply calculates maximum width of a list of `grid::textGrob` objects.
793b75de
 #
30d2c5b2
 # Note it ignores the text rotation.
 #
12e85497
 # == value
30d2c5b2
 # A `grid::unit` object which is in "mm".
12e85497
 #
 # == author
 # Zuguang Gu <z.gu@dkfz.de>
 #
30d2c5b2
 # == seealso
 # `max_text_height` calculates the maximum height of a text vector.
 #
12e85497
 # == example
793b75de
 # x = c("a", "bb", "ccc")
 # max_text_width(x, gp = gpar(fontsize = 10))
12e85497
 #
261bfec2
 max_text_width = function(text, gp = gpar()) {
ce800d69
     if(is.null(text)) {
         return(unit(0, "mm"))
     }
261bfec2
     n = length(text)
     gp = recycle_gp(gp, n)
 
d7a3c7af
     u = max(do.call("unit.c", lapply(seq_len(n), function(i) grobWidth(textGrob(text[i], gp = subset_gp(gp, i))))))
     convertWidth(u, "mm")
076c1ff2
 }
 
 # == title
e27480b9
 # Maximum Height of Text
076c1ff2
 #
 # == param
30d2c5b2
 # -text A vector of text.
ad35494a
 # -gp Graphic parameters for text.
076c1ff2
 #
793b75de
 # == details
64d651fe
 # It simply calculates maximum height of a list of `grid::textGrob` objects.
793b75de
 #
30d2c5b2
 # Note it ignores the text rotation.
 #
12e85497
 # == value
 # A `grid::unit` object.
 #
 # == author
 # Zuguang Gu <z.gu@dkfz.de>
 #
30d2c5b2
 # == seealso
 # `max_text_width` calculates the maximum width of a text vector.
 #
12e85497
 # == example
793b75de
 # x = c("a", "b\nb", "c\nc\nc")
 # max_text_height(x, gp = gpar(fontsize = 10))
12e85497
 #
261bfec2
 max_text_height = function(text, gp = gpar()) {
ce800d69
     if(is.null(text)) {
         return(unit(0, "mm"))
     }
261bfec2
     n = length(text)
     gp = recycle_gp(gp, n)
 
d7a3c7af
     u = max(do.call("unit.c", lapply(seq_len(n), function(i) grobHeight(textGrob(text[i], gp = subset_gp(gp, i))))))
     convertHeight(u, "mm")
261bfec2
 }
 
 dev.null = function(...) {
     pdf(file = NULL, ...)
 }
 
 stop_wrap = function (...) {
cc35650a
     x = paste0(...)
261bfec2
     x = paste(strwrap(x), collapse = "\n")
933e808c
     stop(x, call. = FALSE)
076c1ff2
 }
261bfec2
 
 warning_wrap = function (...) {
cc35650a
     x = paste0(...)
261bfec2
     x = paste(strwrap(x), collapse = "\n")
933e808c
     warning(x, call. = FALSE)
261bfec2
 }
 
 message_wrap = function (...) {
cc35650a
     x = paste0(...)
261bfec2
     x = paste(strwrap(x), collapse = "\n")
     message(x)
 }
 
d7a3c7af
 generate_param_list_fun = function(default) {
     if(!is.list(default)) {
6f036283
         stop_wrap("`default` needs to be a list.")
d7a3c7af
     }
     lt = default
     function(..., list = NULL) {
         if(missing(list)) {
             lt2 = list(...)
         } else {
             lt2 = list
         }
         for(nm in intersect(names(lt), names(lt2))) {
             lt[[nm]] = lt2[[nm]]
         }
         return(lt)
     }
 }
 
 add_vp_name = function(vpname) {
     grid.text(vpname, 0, 1, just = c("left", "top"), gp = gpar(fontsize = 6, col = "red"))
 }
 
 upViewport = function(...) {
a93897a0
     if(ht_global_opt$show_vp) {
402ff791
         grid.rect(gp = gpar(fill = "transparent", col = "black", lty = 3))
         vpname = current.viewport()$name
1ee53830
         if(!grepl("^GRID.VP", vpname)) {
             add_vp_name(vpname)
         }
402ff791
     }
d7a3c7af
     grid::upViewport(...)
 }
 
 popViewport = function(...) {
a93897a0
     if(ht_global_opt$show_vp) {
402ff791
         grid.rect(gp = gpar(fill = "transparent", col = "black", lty = 3))
         vpname = current.viewport()$name
1ee53830
         if(!grepl("^GRID.VP", vpname)) {
             add_vp_name(vpname)
         }
402ff791
     }
d7a3c7af
     grid::popViewport(...)
 }
 
 
 dev.off2 = function () {
     i1 = dev.prev()
     i2 = dev.cur()
f0a25823
 
     if (i1 == 2) {
d7a3c7af
         dev.set(i1)
f0a25823
     } else if(i1 > 2) {
         i11 = dev.prev(i1)
         if(names(i11) == "RStudioGD") {
             dev.set(i11)
         } else {
             dev.set(i1)
         }
     }
d7a3c7af
     dev.off(i2)
 }
402ff791
 
 unit.c = function(...) {
     lt = list(...)
     lt = lt[!sapply(lt, is.null)]
     do.call(grid::unit.c, lt)
 }
 
30d2c5b2
 ">.unit" = function(x, y) {
     if(!identical(attr(x, "unit"), "mm")) {
6f036283
         stop_wrap("x should be in mm unit")
30d2c5b2
     }
     if(!identical(attr(y, "unit"), "mm")) {
6f036283
         stop_wrap("y should be in mm unit")
30d2c5b2
     }
     x[[1]] > y[[1]]
 }
 
 "<.unit" = function(x, y) {
     if(!identical(attr(x, "unit"), "mm")) {
6f036283
         stop_wrap("x should be in mm unit")
30d2c5b2
     }
     if(!identical(attr(y, "unit"), "mm")) {
6f036283
         stop_wrap("y should be in mm unit")
30d2c5b2
     }
     x[[1]] < y[[1]]
 }
402ff791
 
 normalize_graphic_param_to_mat = function(x, nc, nr, name) {
     if(is.matrix(x)) {
         if(nrow(x) == nr && ncol(x) == nc) {
             return(x)
         } else {
c4a66bf9
             stop_wrap(paste0(name, "needs to be a matrix with ", nc, " columns and ", nr, " rows."))
402ff791
         }
     } else {
         if(length(x) == nc) {
ad35494a
             return(matrix(rep(x, each = nr), ncol = nc))
402ff791
         } else if(length(x) == nr) {
ad35494a
             return(matrix(rep(x, times = nc), ncol = nc))
402ff791
         } else if(length(x) == 1) {
ad35494a
             return(matrix(x, ncol = nc, nrow = nr))
402ff791
         } else {
c4a66bf9
             stop_wrap(paste0("Since ", name, " is a vector, it should have length of ", nc, " or ", nr, "."))
402ff791
         }
     }
 }
1bba04a0
 
d9a1645d
 recycle_param = function(x, all_names, default, as.list = FALSE) {
1bba04a0
     n = length(all_names)
d9a1645d
     if(length(x) == 0) {
         if(as.list) {
             rep(list(default), n)
         } else {
             rep(default, n)
         }
     } else if(length(x) == n) {
         if(as.list) {
bfaaa9dd
             x = lapply(1:n, function(i) x[i])
d9a1645d
         }
1bba04a0
         return(x)
     } else {
         nm = names(x)
         if(length(intersect(nm, all_names)) == 0) {
             nm = NULL
         }
         if(is.null(nm)) {
             if(length(x) == 1) {
d9a1645d
                 if(as.list) {
bfaaa9dd
                     x = rep(list(x), n)
d9a1645d
                 } else {
                     x = rep(x, n)
                 }
1bba04a0
             } else {
                 if(length(x) > n) {
                     x = x[1:n]
d9a1645d
                     if(as.list) {
bfaaa9dd
                         x = lapply(1:n, function(i) x[i])
d9a1645d
                     }
1bba04a0
                 } else {
d9a1645d
                     if(as.list) {
                         x = c(lapply(seq_along(x), function(i) x[i], 
                               rep(list(default), n - length(x))))
                     } else {
                         x = c(x, rep(default, n - length(x)))
                     }
1bba04a0
                 }
             }
         } else {
d9a1645d
             if(as.list) {
                 x2 = rep(list(default), n)
                 names(x2) = all_names
                 for(cn in intersect(nm, all_names)) {
                     x2[[cn]] = x[cn]
                 }
                 x = x2
             } else {
                 x2 = structure(rep(default, n), names = all_names)
                 x2[intersect(nm, all_names)] = x[intersect(nm, all_names)]
                 x = x2
             }
1bba04a0
         }
         return(x)
     }
 }
6f11a271
 
 # == title
64d651fe
 # Convert XY in a Parent Viewport
6f11a271
 #
 # == param
64d651fe
 # -u A list of two units which correspond to x and y.
 # -vp_name The name of the parent viewport.
6f11a271
 #
64d651fe
 # == details
 # It converts a coordinate measured in current viewport to the coordinate in a parent viewport.
 #
 # In the conversion, all units are recalculated as absolute units, so if you change the size
 # of the interactive graphic window, you need to rerun the function.
 #
 # == value
 # A list of two units.
 # 
6f11a271
 # == example
 # grid.newpage()
 # pushViewport(viewport(x = 0.5, y = 0.5, width = 0.5, height = 0.5, just = c("left", "bottom")))
64d651fe
 # grid.rect()
 # grid.points(x = unit(2, "cm"), y = unit(2, "cm"), pch = 1)
6f11a271
 # u = list(x = unit(2, "cm"), y = unit(2, "cm"))
64d651fe
 # u2 = getXY_in_parent_vp(u)
 # popViewport()
 # grid.rect(gp = gpar(col = "red"))
 # grid.points(x = u2$x, u2$y, pch = 2)
 getXY_in_parent_vp = function(u, vp_name = "ROOT") {
     if(inherits(u, "unit")) {
         if(length(u) == 2) {
             u = list(x = u[1], y = u[2])
         } else {
             stop_wrap("If `u` is a unit vector, it must have length of 2.")
         }
     }
     if(length(u) != 2) {
         stop_wrap("`u` should be a list of length of 2 (two elements: `x` and `y`).")
     }
     if(is.null(names(u))) {
         names(u) = c("x", "y")
     }
 
6f11a271
     vp = current.viewport()
     current_vp_name = vp$name
     original_vp_name = current_vp_name
64d651fe
     on.exit(seekViewport(original_vp_name))
 
     if(current_vp_name == "ROOT") {
         return(u)
     }
6f11a271
     while(current_vp_name != vp_name) {
 
         if(current_vp_name == "ROOT") {
64d651fe
             stop_wrap(qq("Cannot find a parent viewport with name \"@{vp_name}\"."))
6f11a271
         }
 
         u$x = convertX(u$x, "mm")
         u$y = convertX(u$y, "mm")
64d651fe
         # vp is measured in parent vp
         current_vp_x = vp$x - vp$width*vp$valid.just[1]
         current_vp_y = vp$y - vp$height*vp$valid.just[2]
6f11a271
 
         upViewport(1)
         offset_x = convertX(current_vp_x, "mm")
         offset_y = convertY(current_vp_y, "mm")
         u$x = u$x + offset_x
         u$y = u$y + offset_y
 
         vp = current.viewport()
         current_vp_name = vp$name
     }
 
     return(u)
 }
 
 # == title
64d651fe
 # Get Values in a Matrix by Pair-wise Indices
6f11a271
 #
 # == param
64d651fe
 # -m A matrix or a 3-dimension array.
 # -i Row indices or the indices in the first dimension.
 # -j Column indicies or the indices in the second dimension.
 #
 # == value
 # If ``m`` is a matrix, the value returned is a vector ``c(m[i1, j1], m[i2, j2], ...)```. 
 #
 # If ``m`` is an array, the value returned is a matrix ``rbind(m[i1, j1, ], m[i2, j2, ], ...)```.
6f11a271
 #
 # == example
 # m = matrix(rnorm(100), 10)
 # m2 = m[m > 0]
 # ind = do.call("rbind", lapply(1:10, function(ci) {
 #     i = which(m[, ci] > 0)
 #     cbind(i = i, j = rep(ci, length(i)))
 # }))
 # pindex(m, ind[, 1], ind[, 2])
 # identical(pindex(m, ind[, 1], ind[, 2]), m[m > 0])
 #
 # # 3d array
 # arr = array(1:27, dim = c(3, 3, 3))
 # pindex(arr, 1:2, 2:3)
 # identical(pindex(arr, 1:2, 2:3),
 #    rbind(arr[1, 2, ], arr[2, 3, ]))
 pindex = function(m, i, j) {
64d651fe
 
     if(length(i) == 1) i = rep(i, length(j))
     if(length(j) == 1) j = rep(j, length(i))
     if(length(i) != length(j)) {
         stop_wrap("Length of index i and j should be the same.")
     }
 
6f11a271
     nr = nrow(m)
     nc = ncol(m)
     ind = (j-1)*nr + i
     dm = dim(m)
     if(length(dm) == 2) {
         v = as.vector(m)
         v[ind]
     } else if(length(dm) == 3) {
         v = m
         dim(v) = c(dm[1]*dm[2], dm[3])
         v[ind, , drop = FALSE]
     } else {
c4a66bf9
         stop_wrap("dimension of `m` can only be 2 and 3.")
6f11a271
     }
 }
 
9c4d56c4
 # == title
 # Restore the index vector to index matrix in layer_fun
 #
 # == param
 # -j Column indices directly from ``layer_fun``.
 # -i Row indices directly from ``layer_fun``.
 # -x Position on x-direction directly from ``layer_fun``.
 # -y Position on y-direction directly from ``layer_fun``.
 #
 # == details
 # The values that are sent to ``layer_fun`` are all vectors (for the vectorization
 # of the grid graphic functions), however, the heatmap slice where
 # ``layer_fun`` is applied to, is still represented by a matrix, thus, it would be
 # very convinient if all the arguments in ``layer_fun`` can be converted to the
 # sub-matrix for the current slice. Here, as shown in above example,
 # `restore_matrix` does the job. `restore_matrix` directly accepts the first
 # four argument in ``layer_fun`` and returns an index matrix, where rows and
 # columns correspond to the rows and columns in the current slice, from top to
 # bottom and from left to right. The values in the matrix are the natural order
 # of e.g. vector ``j`` in current slice.
 #
 # For following code:
 #
 #     Heatmap(small_mat, name = "mat", col = col_fun,
 #         row_km = 2, column_km = 2,
 #         layer_fun = function(j, i, x, y, w, h, fill) {
 #             ind_mat = restore_matrix(j, i, x, y)
 #             print(ind_mat)
 #         }
 #     )
 #
 # The first output which is for the top-left slice:
 # 
 #          [,1] [,2] [,3] [,4] [,5]
 #     [1,]    1    4    7   10   13
 #     [2,]    2    5    8   11   14
 #     [3,]    3    6    9   12   15
 #
 # As you see, this is a three-row and five-column index matrix where the first
 # row corresponds to the top row in the slice. The values in the matrix
 # correspond to the natural index (i.e. 1, 2, ...) in ``j``, ``i``, ``x``, ``y``,
 # ... in ``layer_fun``. Now, if we want to add values on the second column in the
 # top-left slice, the code which is put inside ``layer_fun`` would look like:
 #
 #     for(ind in ind_mat[, 2]) {
 #         grid.text(small_mat[i[ind], j[ind]], x[ind], y[ind], ...)
 #     }
 #
 # == example
 # set.seed(123)
 # mat = matrix(rnorm(81), nr = 9)
 # Heatmap(mat, row_km = 2, column_km = 2,
 #     layer_fun = function(j, i, x, y, width, height, fill) {
 #        ind_mat = restore_matrix(j, i, x, y)
 #        print(ind_mat)
 # })
 #
 # set.seed(123)
 # mat = matrix(round(rnorm(81), 2), nr = 9)
 # Heatmap(mat, row_km = 2, column_km = 2,
 #     layer_fun = function(j, i, x, y, width, height, fill) {
 #        ind_mat = restore_matrix(j, i, x, y)
 #        ind = unique(c(ind_mat[2, ], ind_mat[, 3]))
 #        grid.text(pindex(mat, i[ind], j[ind]), x[ind], y[ind])
 # })
 restore_matrix = function(j, i, x, y) {
     x = as.numeric(x)
     y = as.numeric(y)
     od = order(x, rev(y))
     ind = seq_along(i)
     j = j[od]
     i = i[od]
     x = x[od]
     y = y[od]
     ind = ind[od]
     
     nr = length(unique(i))
     nc = length(unique(j))
     # I = matrix(i, nrow = nr, ncol = nc)
     # J = matrix(j, nrow = nr, ncol = nc)
     IND = matrix(ind, nrow = nr, ncol = nc)
     return(IND)
 }
 
bb20f132
 
 unit_with_vp = function(..., vp = current.viewport()$name) {
     u = unit(...)
     attr(u, "viewport") = vp
     return(u)
 }
 
 
028e61a2
 # == title
 # Draw a Single Boxplot
 #
 # == param
1a56796e
 # -value A vector of numeric values.
 # -pos Position of the boxplot.
64d651fe
 # -outline Whether draw outlines?
1a56796e
 # -box_width width of the box.
 # -pch Point type.
 # -size Point size.
 # -gp Graphic parameters.
 # -direction Whether the box is vertical or horizontal.
 #
 # == details
 # All the values are measured with ``native`` coordinate.
 #
64d651fe
 # == example
 # lt = list(rnorm(100), rnorm(100))
 # grid.newpage()
 # pushViewport(viewport(xscale = c(0.5, 2.5), yscale = range(lt)))
 # grid.boxplot(lt[[1]], pos = 1, gp = gpar(fill = "red"))
 # grid.boxplot(lt[[2]], pos = 2, gp = gpar(fill = "green"))
 # popViewport()
cc35650a
 grid.boxplot = function(value, pos, outline = TRUE, box_width = 0.6,
     pch = 1, size = unit(2, "mm"), gp = gpar(fill = "#CCCCCC"), 
028e61a2
     direction = c("vertical", "horizontal")) {
cc35650a
 
028e61a2
     direction = match.arg(direction)[1]
cc35650a
     boxplot_stats = boxplot(value, plot = FALSE)$stats
 
     if(direction == "vertical") {
028e61a2
         grid.rect(x = pos, y = boxplot_stats[2, 1], 
             height = boxplot_stats[4, 1] - boxplot_stats[2, 1], width = 1*box_width, just = "bottom", 
cc35650a
             default.units = "native", gp = gp)
         
028e61a2
         grid.segments(pos - 0.5*box_width, boxplot_stats[5, 1],
                       pos + 0.5*box_width, boxplot_stats[5, 1], 
cc35650a
                       default.units = "native", gp = gp)
028e61a2
         grid.segments(pos, boxplot_stats[5, 1],
                       pos, boxplot_stats[4, 1], 
cc35650a
                       default.units = "native", gp = gp)
028e61a2
         grid.segments(pos, boxplot_stats[1, 1],
                       pos, boxplot_stats[2, 1], 
cc35650a
                       default.units = "native", gp = gp)
028e61a2
         grid.segments(pos - 0.5*box_width, boxplot_stats[1, 1],
                       pos + 0.5*box_width, boxplot_stats[1, 1], 
cc35650a
                       default.units = "native", gp = gp)
028e61a2
         grid.segments(pos - 0.5*box_width, boxplot_stats[3, 1],
                       pos + 0.5*box_width, boxplot_stats[3, 1], 
cc35650a
                       default.units = "native", gp = gp)
         if(outline) {   
028e61a2
             l1 = value > boxplot_stats[5, 1]
8ef41188
             if(sum(l1)) grid.points(x = rep(pos, sum(l1)), y = value[l1], 
028e61a2
                 default.units = "native", gp = gp, pch = pch, size = size)
             l2 = value < boxplot_stats[1, 1]
8ef41188
             if(sum(l2)) grid.points(x = rep(pos, sum(l2)), y = value[l2], 
028e61a2
                 default.units = "native", gp = gp, pch = pch, size = size) 
         }
     } else {
         grid.rect(y = pos, x = boxplot_stats[2, 1], 
             width = boxplot_stats[4, 1] - boxplot_stats[2, 1], height = 1*box_width, just = "left", 
             default.units = "native", gp = gp)
         
         grid.segments(boxplot_stats[5, 1], pos - 0.5*box_width,
                       boxplot_stats[5, 1], pos + 0.5*box_width,
                       default.units = "native", gp = gp)
         grid.segments(boxplot_stats[5, 1], pos,
                       boxplot_stats[4, 1], pos,
                       default.units = "native", gp = gp)
         grid.segments(boxplot_stats[1, 1], pos,
                       boxplot_stats[2, 1], pos,
                       default.units = "native", gp = gp)
         grid.segments(boxplot_stats[1, 1], pos - 0.5*box_width,
                       boxplot_stats[1, 1], pos + 0.5*box_width,
                       default.units = "native", gp = gp)
         grid.segments(boxplot_stats[3, 1], pos - 0.5*box_width,
                       boxplot_stats[3, 1], pos + 0.5*box_width,
                       default.units = "native", gp = gp)
         if(outline) {   
             l1 = value > boxplot_stats[5, 1]
8ef41188
             if(sum(l1)) grid.points(y = rep(pos, sum(l1)), x = value[l1], 
028e61a2
                 default.units = "native", gp = gp, pch = pch, size = size)
             l2 = value < boxplot_stats[1, 1]
8ef41188
             if(sum(l2)) grid.points(y = rep(pos, sum(l2)), x = value[l2], 
028e61a2
                 default.units = "native", gp = gp, pch = pch, size = size) 
cc35650a
         }
028e61a2
     }
 }
 
36c3a378
 random_str = function(k = 1, len = 10) {
     sapply(seq_len(k), function(i) paste(sample(c(letters, LETTERS, 0:9), len), collapse = ""))
991774bb
 }
028e61a2
 
d9a1645d
 
 
 to_unit_str = function(unit) {
     as.character(unit)
 }
 
 to_unit = function(str) {
     d = gsub("[^\\d]+$", "", str, perl = TRUE)
     u = gsub("[\\d.]", "", str, perl = TRUE)
     unit(as.numeric(d), u)
 }
 
4379b65f
 
 resize_matrix = function(mat, nr, nc) {
     w_ratio = nc/ncol(mat)
     h_ratio = nr/nrow(mat)
7185434b
     mat[ ceiling(1:nr / h_ratio), ceiling(1:nc / w_ratio), drop = FALSE]
4379b65f
 }
1f0f3def
 
 
 # == title
 # Adjust positions of rectanglar shapes
 #
 # == param
 # -start position which corresponds to the start (bottom or left) of the rectangle-shapes.
 # -end position which corresponds to the end (top or right) of the rectanglar shapes.
 # -range data ranges (the minimal and maximal values)
 # -range_fixed Whether the range is fixed for ``range`` when adjust the positions?
 #
 # == details
 # This is an improved version of the `circlize::smartAlign`.
 #
 # It adjusts the positions of the rectangular shapes to make them do not overlap
 #
 # == example
 # require(circlize)
 # make_plot = function(pos1, pos2, range) {
 #     oxpd = par("xpd")
 #     par(xpd = NA)
 #     plot(NULL, xlim = c(0, 4), ylim = range, ann = FALSE)
 #     col = rand_color(nrow(pos1), transparency = 0.5)
 #     rect(0.5, pos1[, 1], 1.5, pos1[, 2], col = col)
 #     rect(2.5, pos2[, 1], 3.5, pos2[, 2], col = col)
 #     segments(1.5, rowMeans(pos1), 2.5, rowMeans(pos2))
 #     par(xpd = oxpd)
 # }
 #
 # range = c(0, 10)
 # pos1 = rbind(c(1, 2), c(5, 7))
 # make_plot(pos1, smartAlign2(pos1, range = range), range)
 #
 # range = c(0, 10)
 # pos1 = rbind(c(-0.5, 2), c(5, 7))
 # make_plot(pos1, smartAlign2(pos1, range = range), range)
 #
 # pos1 = rbind(c(-1, 2), c(3, 4), c(5, 6), c(7, 11))
 # pos1 = pos1 + runif(length(pos1), max = 0.3, min = -0.3)
 # omfrow = par("mfrow")
 # par(mfrow = c(3, 3))
 # for(i in 1:9) {
 #     ind = sample(4, 4)
 #     make_plot(pos1[ind, ], smartAlign2(pos1[ind, ], range = range), range)
 # }
 # par(mfrow = omfrow)
 #
 # pos1 = rbind(c(3, 6), c(4, 7))
 # make_plot(pos1, smartAlign2(pos1, range = range), range)
 #
 # pos1 = rbind(c(1, 8), c(3, 10))
 # make_plot(pos1, smartAlign2(pos1, range = range), range)
 # make_plot(pos1, smartAlign2(pos1, range = range, range_fixed = FALSE), range)
 #
 smartAlign2 = function(start, end, range, range_fixed = TRUE) {
 
     if(missing(end)) {
         x1 = start[, 1]
         x2 = start[, 2]
     } else {
         x1 = start
         x2 = end
     }
 
     if(missing(range)) {
         range = range(c(x1, x2))
     }
 
     od = order(x1)
     rk = rank(x1, ties.method = "random")
     x1 = x1[od]
     x2 = x2[od]
     h = x2 - x1
 
     ncluster.before = -1
     ncluster = length(x1)
     i_try = 0
     while(ncluster.before != ncluster) {
         ncluster.before = ncluster
         cluster = rep(0, length(x1))
         i_cluster = 1
         cluster[1] = i_cluster
         for(i in seq_along(x1)[-1]) {
             # overlap with previous one
             if(x1[i] <= x2[i-1]) {  # this means x1 should be sorted increasingly
                 cluster[i] = i_cluster
             } else {
                 i_cluster = i_cluster + 1
                 cluster[i] = i_cluster
             }
         }
         ncluster = length(unique(cluster))
         
         if(ncluster.before == ncluster && i_try > 0) break
c6d85bdc
 
         if(i_try > 100) break
1f0f3def
         
         # tile intervals in each cluster and re-assign x1 and x2
         new_x1 = numeric(length(x1))
         new_x2 = numeric(length(x2))
         for(i_cluster in unique(cluster)) {
             index = which(cluster == i_cluster)
             total_len = sum(x2[index] - x1[index])  # sum of the height in the cluster
             mid = (min(x1[index]) + max(x2[index]))/2
             if(total_len > range[2] - range[1]) {
                 # tp = seq(range[1], range[2], length = length(index) + 1)
                 if(range_fixed) {
                     tp = cumsum(c(0, h[index]/sum(h[index])))*(range[2] - range[1]) + range[1]
                 } else {
                     tp = c(0, cumsum(h[index])) + mid - sum(h[index])/2
                 }
             } else if(mid - total_len/2 < range[1]) { # if it exceed the bottom
                 # tp = seq(range[1], range[1] + total_len, length = length(index) + 1)
                 tp = c(0, cumsum(h[index])) + range[1]
             } else if(mid + total_len/2 > range[2]) {
                 # tp = seq(range[2] - total_len, range[2], length = length(index) + 1)
                 tp = range[2] - rev(c(0, cumsum(h[index])))
             } else {
                 # tp = seq(mid - total_len/2, mid + total_len/2, length = length(index)+1)
                 tp = c(0, cumsum(h[index])) + mid - sum(h[index])/2
             }
             new_x1[index] = tp[-length(tp)]
             new_x2[index] = tp[-1]
         }
         mid = (new_x1 + new_x2)/2
         h = (x2 - x1)
         
         x1 = mid - h/2
         x2 = mid + h/2
 
         i_try = i_try + 1
     }
     
     df = data.frame(start = x1, end = x2)
     df[rk, , drop = FALSE]
 }