From: Benjamin Auder Date: Mon, 16 Dec 2019 11:35:10 +0000 (+0100) Subject: Fix OptimParams + clean test file X-Git-Url: https://git.auder.net/?p=morpheus.git;a=commitdiff_plain;h=44559add0e38058d9ce539c4b91246e4a088f67a Fix OptimParams + clean test file --- diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R index b7d901c..c42e6c5 100644 --- a/pkg/R/optimParams.R +++ b/pkg/R/optimParams.R @@ -104,7 +104,7 @@ setRefClass( list( # p: dimension K-1, need to be completed "p" = c(v[1:(K-1)], 1-sum(v[1:(K-1)])), - "β" = matrix(v[K:(K+d*K-1)], ncol=K), + "β" = t(matrix(v[K:(K+d*K-1)], ncol=d)), "b" = v[(K+d*K):(K+(d+1)*K-1)]) }, @@ -112,12 +112,12 @@ setRefClass( { "Linearize vectors+matrices from list L into a vector" - c(L$p[1:(K-1)], as.double(L$β), L$b) + # β linearized row by row, to match derivatives order + c(L$p[1:(K-1)], as.double(t(L$β)), L$b) }, computeW = function(θ) { - #return (diag(c(rep(6,d), rep(3, d^2), rep(1,d^3)))) require(MASS) dd <- d + d^2 + d^3 M <- Moments(θ) diff --git a/pkg/tests/testthat/test-optimParams.R b/pkg/tests/testthat/test-optimParams.R index 993015f..305c36f 100644 --- a/pkg/tests/testthat/test-optimParams.R +++ b/pkg/tests/testthat/test-optimParams.R @@ -86,130 +86,38 @@ test_that("naive computation provides the same result as vectorized computations } }) -# TODO: test computeW -# computeW = function(θ) -# { -# require(MASS) -# dd <- d + d^2 + d^3 -# M <- Moments(θ) -# Id <- as.double(diag(d)) -# E <- diag(d) -# v1 <- Y * X -# v2 <- Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) ) -# v3 <- Y * t( apply(X, 1, function(Xi) { return (Xi %o% Xi %o% Xi -# - Reduce('+', lapply(1:d, function(j) as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d)) -# - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d)) -# - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) ) -# Wtmp <- matrix(0, nrow=dd, ncol=dd) -# -# -#g <- matrix(nrow=n, ncol=dd); for (i in 1:n) g[i,] = c(v1[i,], v2[i,], v3[i,]) - M -# -# -# -# -# -# -# p <- θ$p -# β <- θ$β -# b <- θ$b -# -# -# -# -## # Random generation of the size of each population in X~Y (unordered) -## classes <- rmultinom(1, n, p) -## -## #d <- nrow(β) -## zero_mean <- rep(0,d) -## id_sigma <- diag(rep(1,d)) -## X <- matrix(nrow=0, ncol=d) -## Y <- c() -## for (i in 1:ncol(β)) #K = ncol(β) -## { -## newXblock <- MASS::mvrnorm(classes[i], zero_mean, id_sigma) -## arg_link <- newXblock %*% β[,i] + b[i] -## probas <- -## if (li == "logit") -## { -## e_arg_link = exp(arg_link) -## e_arg_link / (1 + e_arg_link) -## } -## else #"probit" -## pnorm(arg_link) -## probas[is.nan(probas)] <- 1 #overflow of exp(x) -## X <- rbind(X, newXblock) -## Y <- c( Y, vapply(probas, function(p) (rbinom(1,1,p)), 1) ) -## } -# -# -# -# -# -# -# -# -# Mhatt <- c( -# colMeans(Y * X), -# colMeans(Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) )), -# colMeans(Y * t( apply(X, 1, function(Xi) { return (Xi %o% Xi %o% Xi -# - Reduce('+', lapply(1:d, function(j) as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d)) -# - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d)) -# - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) ) )) -# λ <- sqrt(colSums(β^2)) -# β2 <- apply(β, 2, function(col) col %o% col) -# β3 <- apply(β, 2, function(col) col %o% col %o% col) -# M <- c( -# β %*% (p * .G(li,1,λ,b)), -# β2 %*% (p * .G(li,2,λ,b)), -# β3 %*% (p * .G(li,3,λ,b)) ) -# print(sum(abs(Mhatt - M))) -# -#save(list=c("X", "Y"), file="v2.RData") -# -# -# -# -#browser() -# for (i in 1:n) -# { -# gi <- t(as.matrix(c(v1[i,], v2[i,], v3[i,]) - M)) -# Wtmp <- Wtmp + t(gi) %*% gi / n -# } -# Wtmp -# #MASS::ginv(Wtmp) -# }, -# -# #TODO: compare with R version? -# computeW_orig = function(θ) -# { -# require(MASS) -# dd <- d + d^2 + d^3 -# M <- Moments(θ) -# Omega <- matrix( .C("Compute_Omega", -# X=as.double(X), Y=as.double(Y), M=as.double(M), -# pn=as.integer(n), pd=as.integer(d), -# W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd ) -# Omega -# #MASS::ginv(Omega) #, tol=1e-4) -# }, -# -# Moments = function(θ) -# { -# "Vector of moments, of size d+d^2+d^3" -# -# p <- θ$p -# β <- θ$β -# λ <- sqrt(colSums(β^2)) -# b <- θ$b -# -# # Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1 -# β2 <- apply(β, 2, function(col) col %o% col) -# β3 <- apply(β, 2, function(col) col %o% col %o% col) -# -# c( -# β %*% (p * .G(li,1,λ,b)), -# β2 %*% (p * .G(li,2,λ,b)), -# β3 %*% (p * .G(li,3,λ,b))) -# }, -# +test_that("W computed in C and in R are the same", +{ + # TODO: provide data X,Y + parameters theta + dd <- d + d^2 + d^3 + p <- θ$p + β <- θ$β + λ <- sqrt(colSums(β^2)) + b <- θ$b + β2 <- apply(β, 2, function(col) col %o% col) + β3 <- apply(β, 2, function(col) col %o% col %o% col) + M <- c( + β %*% (p * .G(li,1,λ,b)), + β2 %*% (p * .G(li,2,λ,b)), + β3 %*% (p * .G(li,3,λ,b))) + Id <- as.double(diag(d)) + E <- diag(d) + v1 <- Y * X + v2 <- Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) ) + v3 <- Y * t( apply(X, 1, function(Xi) { return (Xi %o% Xi %o% Xi + - Reduce('+', lapply(1:d, function(j) as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d)) + - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d)) + - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) ) + Omega1 <- matrix(0, nrow=dd, ncol=dd) + for (i in 1:n) + { + gi <- t(as.matrix(c(v1[i,], v2[i,], v3[i,]) - M)) + Omega1 <- Omega1 + t(gi) %*% gi / n + } + Omega2 <- matrix( .C("Compute_Omega", + X=as.double(X), Y=as.double(Y), M=as.double(M), + pn=as.integer(n), pd=as.integer(d), + W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd ) + rg <- range(Omega1 - Omega2) + expect_that(rg[2] - rg[1] <= 1e8) +})