|
| 1 | +#' Perform MDS on landmarks and project other samples to the same space |
| 2 | +#' |
| 3 | +#' @param dist_2lm Distance matrix between the landmarks and all the samples in original dataset |
| 4 | +#' @param ndim The number of dimensions |
| 5 | +#' @param rescale Whether or not to rescale the final dimensionality reduction (recommended) |
| 6 | +#' @param ... Extra params to pass to [irlba::irlba()] |
| 7 | +#' |
| 8 | +#' @importFrom irlba partial_eigen |
| 9 | +#' @importFrom dynutils scale_uniform |
| 10 | +#' |
| 11 | +#' @export |
| 12 | +#' |
| 13 | +#' @examples |
| 14 | +#' library(Matrix) |
| 15 | +#' x <- Matrix::rsparsematrix(10000, 1000, .01) |
| 16 | +#' dist_2lm <- select_landmarks(x) |
| 17 | +#' cmdscale_landmarks(dist_2lm) |
| 18 | +cmdscale_landmarks <- function(dist_2lm, ndim = 3, rescale = TRUE, ...) { |
| 19 | + assert_that( |
| 20 | + !is.null(attr(dist_2lm, "landmark_ix")) || (!is.null(rownames(dist_2lm)) && !is.null(colnames(dist_2lm))), |
| 21 | + is.numeric(ndim), length(ndim) == 1, ndim >= 1, |
| 22 | + is.logical(rescale), length(rescale) == 1, !is.na(rescale) |
| 23 | + ) |
| 24 | + |
| 25 | + ix_lm <- attr(dist_2lm, "landmark_ix") |
| 26 | + if (is.null(ix_lm)) { |
| 27 | + ix_lm <- match(rownames(dist_2lm), colnames(dist_2lm)) |
| 28 | + } |
| 29 | + |
| 30 | + # short hand notations |
| 31 | + x <- dist_2lm[, ix_lm, drop = FALSE]^2 |
| 32 | + n <- as.integer(nrow(x)) |
| 33 | + N <- as.integer(ncol(dist_2lm)) |
| 34 | + |
| 35 | + # double center data |
| 36 | + mu_n <- rowMeans(x) |
| 37 | + mu <- mean(x) |
| 38 | + x_dc <- |
| 39 | + sweep( |
| 40 | + sweep(x, 1, mu_n, "-"), |
| 41 | + 2, mu_n, "-" |
| 42 | + ) + mu |
| 43 | + |
| 44 | + # classical MDS on landmarks |
| 45 | + e <- irlba::partial_eigen(-x_dc / 2, symmetric = TRUE, n = ndim, ...) |
| 46 | + ev <- e$values |
| 47 | + evec <- e$vectors |
| 48 | + ndim1 <- sum(ev > 0) |
| 49 | + if (ndim1 < ndim) { |
| 50 | + warning(gettextf("only %d of the first %d eigenvalues are > 0", ndim1, ndim), domain = NA) |
| 51 | + evec <- evec[, ev > 0, drop = FALSE] |
| 52 | + ev <- ev[ev > 0] |
| 53 | + ndim <- ndim1 |
| 54 | + } |
| 55 | + |
| 56 | + # distance-based triangulation |
| 57 | + points_inv <- evec / rep(sqrt(ev), each = n) |
| 58 | + dimred <- (-t(dist_2lm - rep(mu_n, each = N)) / 2) %*% points_inv |
| 59 | + |
| 60 | + if (rescale) { |
| 61 | + dimred <- dynutils::scale_uniform(dimred) |
| 62 | + } |
| 63 | + |
| 64 | + dimred |
| 65 | +} |
0 commit comments