update
[morpheus.git] / pkg / R / utils.R
CommitLineData
cbd88fe5
BA
1#' normalize
2#'
3#' Normalize a vector or a matrix (by columns), using euclidian norm
4#'
2b3a6af5 5#' @param x Vector or matrix to be normalized
cbd88fe5 6#'
2b3a6af5 7#' @return The normalized matrix (1 column if x is a vector)
cbd88fe5 8#'
2b3a6af5
BA
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)
cbd88fe5 13#' @export
ab35f610 14normalize <- function(x)
cbd88fe5 15{
ab35f610
BA
16 x <- as.matrix(x)
17 norm2 <- sqrt( colSums(x^2) )
2b3a6af5 18 sweep(x, 2, norm2, '/')
cbd88fe5
BA
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#
ab35f610 28.T_I_I_w <- function(Te, w)
cbd88fe5 29{
ab35f610
BA
30 d <- length(w)
31 Ma <- matrix(0,nrow=d,ncol=d)
6dd5c2ac 32 for (j in 1:d)
ab35f610 33 Ma <- Ma + w[j] * Te[,,j]
6dd5c2ac 34 Ma
cbd88fe5
BA
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#
ab35f610 44.Moments_M2 <- function(X, Y)
cbd88fe5 45{
ab35f610
BA
46 n <- nrow(X)
47 d <- ncol(X)
48 M2 <- matrix(0,nrow=d,ncol=d)
6dd5c2ac
BA
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)
cbd88fe5
BA
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#
ab35f610 60.Moments_M3 <- function(X, Y)
cbd88fe5 61{
ab35f610
BA
62 n <- nrow(X)
63 d <- ncol(X)
64 M3 <- array(0,dim=c(d,d,d))
6dd5c2ac
BA
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) )
cbd88fe5
BA
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#'
2b3a6af5
BA
77#' @examples
78#' X <- matrix(rnorm(100), ncol=2)
79#' Y <- rbinom(100, 1, .5)
80#' M <- computeMoments(X, Y)
81#'
cbd88fe5
BA
82#' @export
83computeMoments = function(X, Y)
6dd5c2ac 84 list( colMeans(Y * X), .Moments_M2(X,Y), .Moments_M3(X,Y) )
cbd88fe5
BA
85
86# Find the optimal assignment (permutation) between two sets (minimize cost)
87#
2b3a6af5
BA
88# @param distances The distances matrix, in columns
89# (distances[i,j] is distance between i and j)
cbd88fe5
BA
90#
91# @return A permutation minimizing cost
92#
ab35f610 93.hungarianAlgorithm <- function(distances)
cbd88fe5 94{
ab35f610 95 n <- nrow(distances)
6dd5c2ac
BA
96 .C("hungarianAlgorithm", distances=as.double(distances), pn=as.integer(n),
97 assignment=integer(n), PACKAGE="morpheus")$assignment
cbd88fe5
BA
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
2b3a6af5 105#' @param ref A reference matrix to align other matrices with
cbd88fe5
BA
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#'
2b3a6af5
BA
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#'
cbd88fe5 122#' @export
ab35f610 123alignMatrices <- function(Ms, ref, ls_mode=c("exact","approx1","approx2"))
cbd88fe5 124{
2b3a6af5
BA
125 if (!is.matrix(ref) || any(is.na(ref)))
126 stop("ref: matrix, no NAs")
127 ls_mode <- match.arg(ls_mode)
cbd88fe5 128
6dd5c2ac 129 K <- ncol(Ms[[1]])
6dd5c2ac 130 L <- length(Ms)
2b3a6af5 131 for (i in 1:L)
6dd5c2ac 132 {
ab35f610 133 m <- Ms[[i]] #shorthand
cbd88fe5 134
6dd5c2ac
BA
135 if (ls_mode == "exact")
136 {
137 #distances[i,j] = distance between m column i and ref column j
2b3a6af5 138 distances = apply( ref, 2, function(col) ( sqrt(colSums((m-col)^2)) ) )
6dd5c2ac
BA
139 assignment = .hungarianAlgorithm(distances)
140 col <- m[,assignment]
2b3a6af5 141 Ms[[i]] <- col
6dd5c2ac
BA
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,
2b3a6af5 155 function(col) ( sqrt(sum((col - ref[,j])^2)) ) )
6dd5c2ac
BA
156 }
157 else #approx2
158 {
2b3a6af5 159 apply(as.matrix(ref[,available_indices]), 2,
6dd5c2ac
BA
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] ]
2b3a6af5 166 Ms[[i]][,j] <- col
6dd5c2ac
BA
167 }
168 else #approx2
169 {
170 col <- available_indices[indMin]
2b3a6af5 171 Ms[[i]][,col] <- m[,j]
6dd5c2ac
BA
172 }
173 available_indices = available_indices[-indMin]
174 }
175 }
6dd5c2ac
BA
176 }
177 Ms
cbd88fe5 178}