pkg/R/utils.R

   1 #' normalize
   2 #'
   3 #' Normalize a vector or a matrix (by columns), using euclidian norm
   4 #'
   5 #' @param x Vector or matrix to be normalized
   6 #'
   7 #' @return The normalized matrix (1 column if x is a vector)
   8 #'
   9 #' @examples
  10 #' x <- matrix(c(1,2,-1,3), ncol=2)
  11 #' normalize(x) #column 1 is 1/sqrt(5) (1 2),
  12 #'              #and column 2 is 1/sqrt(10) (-1, 3)
  13 #' @export
  14 normalize <- function(x)
  15 {
  16   x <- as.matrix(x)
  17   norm2 <- sqrt( colSums(x^2) )
  18   sweep(x, 2, norm2, '/')
  19 }
  20
  21 # Computes a tensor-vector product
  22 #
  23 # @param Te third-order tensor (size dxdxd)
  24 # @param w vector of size d
  25 #
  26 # @return Matrix of size dxd
  27 #
  28 .T_I_I_w <- function(Te, w)
  29 {
  30   d <- length(w)
  31   Ma <- matrix(0,nrow=d,ncol=d)
  32   for (j in 1:d)
  33     Ma <- Ma + w[j] * Te[,,j]
  34   Ma
  35 }
  36
  37 # Computes the second-order empirical moment between input X and output Y
  38 #
  39 # @param X matrix of covariates (of size n*d)
  40 # @param Y vector of responses (of size n)
  41 #
  42 # @return Matrix of size dxd
  43 #
  44 .Moments_M2 <- function(X, Y)
  45 {
  46   n <- nrow(X)
  47   d <- ncol(X)
  48   M2 <- matrix(0,nrow=d,ncol=d)
  49   matrix( .C("Moments_M2", X=as.double(X), Y=as.double(Y), pn=as.integer(n),
  50     pd=as.integer(d), M2=as.double(M2), PACKAGE="morpheus")$M2, nrow=d, ncol=d)
  51 }
  52
  53 # Computes the third-order empirical moment between input X and output Y
  54 #
  55 # @param X matrix of covariates (of size n*d)
  56 # @param Y vector of responses (of size n)
  57 #
  58 # @return Array of size dxdxd
  59 #
  60 .Moments_M3 <- function(X, Y)
  61 {
  62   n <- nrow(X)
  63   d <- ncol(X)
  64   M3 <- array(0,dim=c(d,d,d))
  65   array( .C("Moments_M3", X=as.double(X), Y=as.double(Y), pn=as.integer(n),
  66     pd=as.integer(d), M3=as.double(M3), PACKAGE="morpheus")$M3, dim=c(d,d,d) )
  67 }
  68
  69 #' computeMoments
  70 #'
  71 #' Compute cross-moments of order 1,2,3 from X,Y
  72 #'
  73 #' @inheritParams computeMu
  74 #'
  75 #' @return A list L where L[[i]] is the i-th cross-moment
  76 #'
  77 #' @examples
  78 #' X <- matrix(rnorm(100), ncol=2)
  79 #' Y <- rbinom(100, 1, .5)
  80 #' M <- computeMoments(X, Y)
  81 #'
  82 #' @export
  83 computeMoments = function(X, Y)
  84   list( colMeans(Y * X), .Moments_M2(X,Y), .Moments_M3(X,Y) )
  85
  86 # Find the optimal assignment (permutation) between two sets (minimize cost)
  87 #
  88 # @param distances The distances matrix, in columns
  89 #   (distances[i,j] is distance between i and j)
  90 #
  91 # @return A permutation minimizing cost
  92 #
  93 .hungarianAlgorithm <- function(distances)
  94 {
  95   n <- nrow(distances)
  96   .C("hungarianAlgorithm", distances=as.double(distances), pn=as.integer(n),
  97     assignment=integer(n), PACKAGE="morpheus")$assignment
  98 }
  99
 100 #' alignMatrices
 101 #'
 102 #' Align a set of parameters matrices, with potential permutations.
 103 #'
 104 #' @param Ms A list of matrices, all of same size DxK
 105 #' @param ref A reference matrix to align other matrices with
 106 #' @param ls_mode How to compute the labels assignment: "exact" for exact algorithm
 107 #'   (default, but might be time-consuming, complexity is O(K^3) ), or "approx1", or
 108 #'   "approx2" to apply a greedy matching algorithm (heuristic) which for each column in
 109 #'   reference (resp. in current row) compare to all unassigned columns in current row
 110 #'   (resp. in reference)
 111 #'
 112 #' @return The aligned list (of matrices), of same size as Ms
 113 #'
 114 #' @examples
 115 #' m1 <- matrix(c(1,1,0,0),ncol=2)
 116 #' m2 <- matrix(c(0,0,1,1),ncol=2)
 117 #' ref <- m1
 118 #' Ms <- list(m1, m2, m1, m2)
 119 #' a <- alignMatrices(Ms, ref, "exact")
 120 #' # a[[i]] is expected to contain m1 for all i
 121 #'
 122 #' @export
 123 alignMatrices <- function(Ms, ref, ls_mode=c("exact","approx1","approx2"))
 124 {
 125   if (!is.matrix(ref) || any(is.na(ref)))
 126     stop("ref: matrix, no NAs")
 127   ls_mode <- match.arg(ls_mode)
 128
 129   K <- ncol(Ms[[1]])
 130   L <- length(Ms)
 131   for (i in 1:L)
 132   {
 133     m <- Ms[[i]] #shorthand
 134
 135     if (ls_mode == "exact")
 136     {
 137       #distances[i,j] = distance between m column i and ref column j
 138       distances = apply( ref, 2, function(col) ( sqrt(colSums((m-col)^2)) ) )
 139       assignment = .hungarianAlgorithm(distances)
 140       col <- m[,assignment]
 141       Ms[[i]] <- col
 142     }
 143     else
 144     {
 145       # Greedy matching:
 146       #   approx1: li[[i]][,j] is assigned to m[,k] minimizing dist(li[[i]][,j],m[,k'])
 147       #   approx2: m[,j] is assigned to li[[i]][,k] minimizing dist(m[,j],li[[i]][,k'])
 148       available_indices = 1:K
 149       for (j in 1:K)
 150       {
 151         distances =
 152           if (ls_mode == "approx1")
 153           {
 154             apply(as.matrix(m[,available_indices]), 2,
 155               function(col) ( sqrt(sum((col - ref[,j])^2)) ) )
 156           }
 157           else #approx2
 158           {
 159             apply(as.matrix(ref[,available_indices]), 2,
 160               function(col) ( sqrt(sum((col - m[,j])^2)) ) )
 161           }
 162         indMin = which.min(distances)
 163         if (ls_mode == "approx1")
 164         {
 165           col <- m[ , available_indices[indMin] ]
 166           Ms[[i]][,j] <- col
 167         }
 168         else #approx2
 169         {
 170           col <- available_indices[indMin]
 171           Ms[[i]][,col] <- m[,j]
 172         }
 173         available_indices = available_indices[-indMin]
 174       }
 175     }
 176   }
 177   Ms
 178 }