R/tSNE_internal.R

Defines functions whiten x2p Hbeta

Hbeta <- function(D, beta){
  P = exp(-D * beta)
  sumP = sum(P, na.rm =T )
  if (sumP == 0) {
    H = 0
    P = D * 0
  }else {
    H = log(sumP) + beta * sum(D %*% P, na.rm=T)/sumP
    P = P/sumP
  }
  r = {
  }
  r$H = H
  r$P = P
  return(r)
}
# Get Probability distribution (P) from data (X)
x2p <- function(X, perplexity = 15, tol = 1e-05){
  if (any(class(X) %in% "dist")) {
    D = X
    n = attr(D, "Size")
  }else {
    D = dist(X)
    n = attr(D, "Size")
  }
  D = as.matrix(D)
  P = matrix(0, n, n)
  beta = rep(1, n)
  logU = log(perplexity)
  for (i in 1:n) {
    betamin = -Inf
    betamax = Inf
    Di = D[i, -i]
    hbeta = Hbeta(Di, beta[i])
    H = hbeta$H
    thisP = hbeta$P
    Hdiff = H - logU
    tries = 0
    while (abs(Hdiff) > tol && tries < 50) {
      if (Hdiff > 0) {
        betamin = beta[i]
        if (is.infinite(betamax))
          beta[i] = beta[i] * 2
        else beta[i] = (beta[i] + betamax)/2
      }
      else {
        betamax = beta[i]
        if (is.infinite(betamin))
          beta[i] = beta[i]/2
        else beta[i] = (beta[i] + betamin)/2
      }
      hbeta = Hbeta(Di, beta[i])
      H = hbeta$H
      thisP = hbeta$P
      Hdiff = H - logU
      tries = tries + 1
    }
    P[i, -i] = thisP
  }
  r = {
  }
  r$P = P
  r$beta = beta
  sigma = sqrt(1/beta)
  message("sigma summary: ", paste(names(summary(sigma)),
                                   ":", summary(sigma), "|", collapse = ""))
  return(r)
}

# PCA -- Whiten
whiten <- function(X, row.norm = FALSE, verbose = FALSE, n.comp = ncol(X)){
  n.comp
  if (verbose)
    message("Centering")
  n = nrow(X)
  p = ncol(X)
  X <- scale(X, scale = FALSE)
  X <- if (row.norm)
    t(scale(X, scale = row.norm))
  else t(X)
  if (verbose)
    message("Whitening")
  V <- X %*% t(X)/n
  s <- La.svd(V)
  D <- diag(c(1/sqrt(s$d)))
  K <- D %*% t(s$u)
  K <- matrix(K[1:n.comp, ], n.comp, p)
  X = t(K %*% X)
  return(X)
}
theorod93/multiSNE documentation built on April 22, 2024, 11:16 a.m.