From ab35f6102896a49e86e853262c0650faa2931638 Mon Sep 17 00:00:00 2001 From: Benjamin Auder <benjamin.auder@somewhere> Date: Mon, 7 Jun 2021 17:27:23 +0200 Subject: [PATCH] Adjustments + bugs fixing --- README.md | 3 + ...orDecompositions_Anandkumar.et.al_2014.pdf | 1 + ...nsorFactorization_Kuleshov.et.al_2015.pdf} | 0 pkg/DESCRIPTION | 5 +- pkg/R/computeMu.R | 50 +++++----- pkg/R/multiRun.R | 30 +++--- pkg/R/optimParams.R | 11 ++- pkg/R/plot.R | 17 +++- pkg/R/sampleIO.R | 6 +- pkg/R/utils.R | 38 ++++---- pkg/src/functions.c | 6 +- pkg/tests/testthat.R | 6 -- pkg/tests/testthat/test-Moments.R | 2 - pkg/tests/testthat/test-alignMatrices.R | 4 +- pkg/tests/testthat/test-computeMu.R | 6 +- pkg/tests/testthat/test-hungarianAlgorithm.R | 2 - pkg/tests/testthat/test-jointDiag.R | 6 +- pkg/tests/testthat/test-optimParams.R | 91 +++++++++---------- reports/local_run.sh | 5 +- 19 files changed, 149 insertions(+), 140 deletions(-) create mode 100644 doc/TensorDecompositions_Anandkumar.et.al_2014.pdf rename doc/{Tensor-Kuleshov2015.pdf => TensorFactorization_Kuleshov.et.al_2015.pdf} (100%) diff --git a/README.md b/README.md index dbfd60a..c3d2c9f 100644 --- a/README.md +++ b/README.md @@ -12,6 +12,9 @@ NOTE: greek unicode letters are used in the code, because it's much nicer to wri and Σ than lambda, beta and Sigma - and also more importantly because it works well :) ...However CRAN demands ASCII files - therefore the presence of to-cran.sh script. +Warning: tests won't run well on the non-UTF8 package. Run ./to-cran.sh and +then inside pkg-cran call devtools::test(). + ## Example Install the package, then ?multiRun is a possible starting point. diff --git a/doc/TensorDecompositions_Anandkumar.et.al_2014.pdf b/doc/TensorDecompositions_Anandkumar.et.al_2014.pdf new file mode 100644 index 0000000..90a7f76 --- /dev/null +++ b/doc/TensorDecompositions_Anandkumar.et.al_2014.pdf @@ -0,0 +1 @@ +#$# git-fat 00f581aca4dd370615fa0ea99e730d6dd42fe347 585618 diff --git a/doc/Tensor-Kuleshov2015.pdf b/doc/TensorFactorization_Kuleshov.et.al_2015.pdf similarity index 100% rename from doc/Tensor-Kuleshov2015.pdf rename to doc/TensorFactorization_Kuleshov.et.al_2015.pdf diff --git a/pkg/DESCRIPTION b/pkg/DESCRIPTION index 2a20240..0a642b6 100644 --- a/pkg/DESCRIPTION +++ b/pkg/DESCRIPTION @@ -21,10 +21,10 @@ Suggests: devtools, flexmix, parallel, - testthat, + testthat (>= 3.0.0), roxygen2 License: MIT + file LICENSE -RoxygenNote: 7.0.2 +RoxygenNote: 7.1.1 URL: https://github.com/yagu0/morpheus Collate: 'utils.R' @@ -34,3 +34,4 @@ Collate: 'optimParams.R' 'plot.R' 'sampleIO.R' +Config/testthat/edition: 3 diff --git a/pkg/R/computeMu.R b/pkg/R/computeMu.R index f961fde..c1a9c51 100644 --- a/pkg/R/computeMu.R +++ b/pkg/R/computeMu.R @@ -3,6 +3,8 @@ #' Estimate the normalized columns μ of the β matrix parameter in a mixture of #' logistic regressions models, with a spectral method described in the package vignette. #' +#' @name computeMu +#' #' @param X Matrix of input data (size nxd) #' @param Y Vector of binary outputs (size n) #' @param optargs List of optional argument: @@ -21,28 +23,28 @@ #' and \code{generateSampleIO} for I/O random generation. #' #' @examples -#' io = generateSampleIO(10000, 1/2, matrix(c(1,0,0,1),ncol=2), c(0,0), "probit") -#' μ = computeMu(io$X, io$Y, list(K=2)) #or just X and Y for estimated K +#' io <- generateSampleIO(10000, 1/2, matrix(c(1,0,0,1),ncol=2), c(0,0), "probit") +#' μ <- computeMu(io$X, io$Y, list(K=2)) #or just X and Y for estimated K #' #' @export -computeMu = function(X, Y, optargs=list()) +computeMu <- function(X, Y, optargs=list()) { if (!is.matrix(X) || !is.numeric(X) || any(is.na(X))) stop("X: real matrix, no NA") - n = nrow(X) - d = ncol(X) + n <- nrow(X) + d <- ncol(X) if (!is.numeric(Y) || length(Y)!=n || any(Y!=0 & Y!=1)) stop("Y: vector of 0 and 1, size nrow(X), no NA") if (!is.list(optargs)) stop("optargs: list") # Step 0: Obtain the empirically estimated moments tensor, estimate also K - M = if (is.null(optargs$M)) computeMoments(X,Y) else optargs$M - K = optargs$K + M <- if (is.null(optargs$M)) computeMoments(X,Y) else optargs$M + K <- optargs$K if (is.null(K)) { # TODO: improve this basic heuristic - Σ = svd(M[[2]])$d + Σ <- svd(M[[2]])$d large_ratio <- ( abs(Σ[-d] / Σ[-1]) > 3 ) K <- if (any(large_ratio)) max(2, which.min(large_ratio)) else d } @@ -50,24 +52,24 @@ computeMu = function(X, Y, optargs=list()) stop("K: integer >= 2, <= d") # Step 1: generate a family of d matrices to joint-diagonalize to increase robustness - d = ncol(X) - fixed_design = FALSE - jd_nvects = ifelse(!is.null(optargs$jd_nvects), optargs$jd_nvects, 0) + d <- ncol(X) + fixed_design <- FALSE + jd_nvects <- ifelse(!is.null(optargs$jd_nvects), optargs$jd_nvects, 0) if (jd_nvects == 0) { - jd_nvects = d - fixed_design = TRUE + jd_nvects <- d + fixed_design <- TRUE } - M2_t = array(dim=c(d,d,jd_nvects)) + M2_t <- array(dim=c(d,d,jd_nvects)) for (i in seq_len(jd_nvects)) { - rho = if (fixed_design) c(rep(0,i-1),1,rep(0,d-i)) else normalize( rnorm(d) ) - M2_t[,,i] = .T_I_I_w(M[[3]],rho) + rho <- if (fixed_design) c(rep(0,i-1),1,rep(0,d-i)) else normalize( rnorm(d) ) + M2_t[,,i] <- .T_I_I_w(M[[3]],rho) } # Step 2: obtain factors u_i (and their inverse) from the joint diagonalisation of M2_t - jd_method = ifelse(!is.null(optargs$jd_method), optargs$jd_method, "uwedge") - V = + jd_method <- ifelse(!is.null(optargs$jd_method), optargs$jd_method, "uwedge") + V <- if (jd_nvects > 1) { # NOTE: increasing itermax does not help to converge, thus we suppress warnings suppressWarnings({jd = jointDiag::ajd(M2_t, method=jd_method)}) @@ -77,17 +79,17 @@ computeMu = function(X, Y, optargs=list()) eigen(M2_t[,,1])$vectors # Step 3: obtain final factors from joint diagonalisation of T(I,I,u_i) - M2_t = array(dim=c(d,d,K)) + M2_t <- array(dim=c(d,d,K)) for (i in seq_len(K)) - M2_t[,,i] = .T_I_I_w(M[[3]],V[,i]) + M2_t[,,i] <- .T_I_I_w(M[[3]],V[,i]) suppressWarnings({jd = jointDiag::ajd(M2_t, method=jd_method)}) - U = if (jd_method=="uwedge") MASS::ginv(jd$B) else jd$A - μ = normalize(U[,1:K]) + U <- if (jd_method=="uwedge") MASS::ginv(jd$B) else jd$A + μ <- normalize(U[,1:K]) # M1 also writes M1 = sum_k coeff_k * μ_k, where coeff_k >= 0 # ==> search decomposition of vector M1 onto the (truncated) basis μ (of size dxK) # This is a linear system μ %*% C = M1 with C of size K ==> C = psinv(μ) %*% M1 - C = MASS::ginv(μ) %*% M[[1]] - μ[,C < 0] = - μ[,C < 0] + C <- MASS::ginv(μ) %*% M[[1]] + μ[,C < 0] <- - μ[,C < 0] μ } diff --git a/pkg/R/multiRun.R b/pkg/R/multiRun.R index 5167535..59e2483 100644 --- a/pkg/R/multiRun.R +++ b/pkg/R/multiRun.R @@ -7,6 +7,8 @@ #' basis (exact same inputs). It's even possible to compare methods on some #' deterministic design of experiments. #' +#' @name multiRun +#' #' @param fargs List of arguments for the estimation functions #' @param estimParams List of nf function(s) to apply on fargs #' @param prepareArgs Prepare arguments for the functions inside estimParams @@ -36,8 +38,8 @@ #' library(flexmix) #' ind <- fargs$ind #' K <- fargs$K -#' dat = as.data.frame( cbind(fargs$Y[ind],fargs$X[ind,]) ) -#' out = refit( flexmix( cbind(V1, 1 - V1) ~ 0+., data=dat, k=K, +#' dat <- as.data.frame( cbind(fargs$Y[ind],fargs$X[ind,]) ) +#' out <- refit( flexmix( cbind(V1, 1 - V1) ~ 0+., data=dat, k=K, #' model=FLXMRglm(family="binomial") ) ) #' normalize( matrix(out@@coef[1:(ncol(fargs$X)*K)], ncol=K) ) #' } ), @@ -52,7 +54,7 @@ #' res[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact") #' #' # Monte-Carlo + optimParams from X,Y, morpheus VS flexmix -#' res <- multiRun(list(n=1000,p=1/2,β=β,b=c(0,0),link="logit"),list( +#' res <- multiRun(list(n=1000,p=1/2,β=β,b=c(0,0),link="logit"), list( #' # morpheus #' function(fargs) { #' library(morpheus) @@ -66,19 +68,19 @@ #' library(flexmix) #' K <- fargs$K #' dat <- as.data.frame( cbind(fargs$Y,fargs$X) ) -#' out <- refit( flexmix( cbind(V1, 1 - V1) ~ 0+., data=dat, k=K, +#' out <- refit( flexmix( cbind(V1, 1 - V1) ~ ., data=dat, k=K, #' model=FLXMRglm(family="binomial") ) ) #' sapply( seq_len(K), function(i) -#' as.double( out@@components[[1]][[i]][,1] ) ) +#' as.double( out@@components[[1]][[i]][2:(1+ncol(fargs$X)),1] ) ) #' } ), #' prepareArgs = function(fargs,index) { #' library(morpheus) -#' io = generateSampleIO(fargs$n, fargs$p, fargs$β, fargs$b, fargs$link) -#' fargs$X = io$X -#' fargs$Y = io$Y -#' fargs$K = ncol(fargs$β) -#' fargs$link = fargs$link -#' fargs$M = computeMoments(io$X,io$Y) +#' io <- generateSampleIO(fargs$n, fargs$p, fargs$β, fargs$b, fargs$link) +#' fargs$X <- io$X +#' fargs$Y <- io$Y +#' fargs$K <- ncol(fargs$β) +#' fargs$link <- fargs$link +#' fargs$M <- computeMoments(io$X,io$Y) #' fargs #' }, N=10, ncores=3) #' for (i in 1:2) @@ -118,13 +120,13 @@ multiRun <- function(fargs, estimParams, if (ncores > 1) { - cl = parallel::makeCluster(ncores, outfile="") + cl <- parallel::makeCluster(ncores, outfile="") parallel::clusterExport(cl, c("fargs","verbose"), environment()) - list_res = parallel::clusterApplyLB(cl, 1:N, estimParamAtIndex) + list_res <- parallel::clusterApplyLB(cl, 1:N, estimParamAtIndex) parallel::stopCluster(cl) } else - list_res = lapply(1:N, estimParamAtIndex) + list_res <- lapply(1:N, estimParamAtIndex) # De-interlace results: output one list per function nf <- length(estimParams) diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R index 175150f..c050e63 100644 --- a/pkg/R/optimParams.R +++ b/pkg/R/optimParams.R @@ -1,5 +1,9 @@ +#' optimParams +#' #' Wrapper function for OptimParams class #' +#' @name optimParams +#' #' @param X Data matrix of covariables #' @param Y Output as a binary vector #' @param K Number of populations. @@ -24,7 +28,7 @@ #' # Optimize parameters from estimated μ #' io <- generateSampleIO(100, #' 1/2, matrix(c(1,-2,3,1),ncol=2), c(0,0), "logit") -#' μ = computeMu(io$X, io$Y, list(K=2)) +#' μ <- computeMu(io$X, io$Y, list(K=2)) #' o <- optimParams(io$X, io$Y, 2, "logit") #' \donttest{ #' θ0 <- list(p=1/2, β=μ, b=c(0,0)) @@ -274,7 +278,8 @@ setRefClass( { stop("θ0$p: length K-1, no NA, positive integers, sum to <= 1") } - if (is.null(θ0$b)) + # NOTE: [["b"]] instead of $b because $b would match $beta (in pkg-cran) + if (is.null(θ0[["b"]])) θ0$b = rep(0, K) else if (!is.numeric(θ0$b) || length(θ0$b) != K || any(is.na(θ0$b))) stop("θ0$b: length K, no NA") @@ -282,7 +287,7 @@ setRefClass( # (Re)Set W to identity, to allow several run from the same object W <<- if (is.null(userW)) diag(d+d^2+d^3) else userW - #NOTE: loopMax = 3 seems to not improve the final results. + # NOTE: loopMax = 3 seems to not improve the final results. loopMax <- ifelse(is.null(userW), 2, 1) x_init <- linArgs(θ0) for (loop in 1:loopMax) diff --git a/pkg/R/plot.R b/pkg/R/plot.R index da4fae5..da8f8bf 100644 --- a/pkg/R/plot.R +++ b/pkg/R/plot.R @@ -22,6 +22,8 @@ #' #' Plot compared histograms of a single parameter (scalar) #' +#' @name plotHist +#' #' @param mr Output of multiRun(), list of lists of functions results #' @param x Row index of the element inside the aggregated parameter #' @param y Column index of the element inside the aggregated parameter @@ -41,7 +43,7 @@ plotHist <- function(mr, x, y, ...) { params <- .extractParam(mr, x, y) - L = length(params) + L <- length(params) # Plot histograms side by side par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1)) args <- list(...) @@ -53,11 +55,18 @@ plotHist <- function(mr, x, y, ...) } } +# NOTE: roxygen2 bug, "@inheritParams plotHist" fails in next header: + #' plotBox #' #' Draw compared boxplots of a single parameter (scalar) #' -#' @inheritParams plotHist +#' @name plotBox +#' +#' @param mr Output of multiRun(), list of lists of functions results +#' @param x Row index of the element inside the aggregated parameter +#' @param y Column index of the element inside the aggregated parameter +#' @param ... Additional graphical parameters (xlab, ylab, ...) #' #' @examples #' \donttest{ @@ -73,7 +82,7 @@ plotHist <- function(mr, x, y, ...) plotBox <- function(mr, x, y, ...) { params <- .extractParam(mr, x, y) - L = length(params) + L <- length(params) # Plot boxplots side by side par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1)) args <- list(...) @@ -91,6 +100,8 @@ plotBox <- function(mr, x, y, ...) #' Note that the drawing does not correspond to a function; it is just a #' convenient way to visualize the estimated parameters. #' +#' @name plotCoefs +#' #' @param mr List of parameters matrices #' @param params True value of the parameters matrix #' @param ... Additional graphical parameters diff --git a/pkg/R/sampleIO.R b/pkg/R/sampleIO.R index 230b63c..9419c32 100644 --- a/pkg/R/sampleIO.R +++ b/pkg/R/sampleIO.R @@ -4,6 +4,8 @@ #' into K groups of proportions p. Inside one group, the probability law P(Y=1) is #' described by the corresponding column parameter in the matrix β + intercept b. #' +#' @name generateSampleIO +#' #' @param n Number of individuals #' @param p Vector of K(-1) populations relative proportions (sum (<)= 1) #' @param β Vectors of model parameters for each population, of size dxK @@ -24,10 +26,10 @@ #' io$index[1] #number of the group of X[1,] and Y[1] (in 1...K) #' #' @export -generateSampleIO = function(n, p, β, b, link) +generateSampleIO <- function(n, p, β, b, link) { # Check arguments - tryCatch({n = as.integer(n)}, error=function(e) stop("Cannot convert n to integer")) + tryCatch({n <- as.integer(n)}, error=function(e) stop("Cannot convert n to integer")) if (length(n) > 1) warning("n is a vector but should be scalar: only first element used") if (n <= 0) diff --git a/pkg/R/utils.R b/pkg/R/utils.R index 9d026da..5b9d999 100644 --- a/pkg/R/utils.R +++ b/pkg/R/utils.R @@ -11,10 +11,10 @@ #' normalize(x) #column 1 is 1/sqrt(5) (1 2), #' #and column 2 is 1/sqrt(10) (-1, 3) #' @export -normalize = function(x) +normalize <- function(x) { - x = as.matrix(x) - norm2 = sqrt( colSums(x^2) ) + x <- as.matrix(x) + norm2 <- sqrt( colSums(x^2) ) sweep(x, 2, norm2, '/') } @@ -25,12 +25,12 @@ normalize = function(x) # # @return Matrix of size dxd # -.T_I_I_w = function(Te, w) +.T_I_I_w <- function(Te, w) { - d = length(w) - Ma = matrix(0,nrow=d,ncol=d) + d <- length(w) + Ma <- matrix(0,nrow=d,ncol=d) for (j in 1:d) - Ma = Ma + w[j] * Te[,,j] + Ma <- Ma + w[j] * Te[,,j] Ma } @@ -41,11 +41,11 @@ normalize = function(x) # # @return Matrix of size dxd # -.Moments_M2 = function(X, Y) +.Moments_M2 <- function(X, Y) { - n = nrow(X) - d = ncol(X) - M2 = matrix(0,nrow=d,ncol=d) + n <- nrow(X) + d <- ncol(X) + M2 <- matrix(0,nrow=d,ncol=d) matrix( .C("Moments_M2", X=as.double(X), Y=as.double(Y), pn=as.integer(n), pd=as.integer(d), M2=as.double(M2), PACKAGE="morpheus")$M2, nrow=d, ncol=d) } @@ -57,11 +57,11 @@ normalize = function(x) # # @return Array of size dxdxd # -.Moments_M3 = function(X, Y) +.Moments_M3 <- function(X, Y) { - n = nrow(X) - d = ncol(X) - M3 = array(0,dim=c(d,d,d)) + n <- nrow(X) + d <- ncol(X) + M3 <- array(0,dim=c(d,d,d)) array( .C("Moments_M3", X=as.double(X), Y=as.double(Y), pn=as.integer(n), pd=as.integer(d), M3=as.double(M3), PACKAGE="morpheus")$M3, dim=c(d,d,d) ) } @@ -90,9 +90,9 @@ computeMoments = function(X, Y) # # @return A permutation minimizing cost # -.hungarianAlgorithm = function(distances) +.hungarianAlgorithm <- function(distances) { - n = nrow(distances) + n <- nrow(distances) .C("hungarianAlgorithm", distances=as.double(distances), pn=as.integer(n), assignment=integer(n), PACKAGE="morpheus")$assignment } @@ -120,7 +120,7 @@ computeMoments = function(X, Y) #' # a[[i]] is expected to contain m1 for all i #' #' @export -alignMatrices = function(Ms, ref, ls_mode=c("exact","approx1","approx2")) +alignMatrices <- function(Ms, ref, ls_mode=c("exact","approx1","approx2")) { if (!is.matrix(ref) || any(is.na(ref))) stop("ref: matrix, no NAs") @@ -130,7 +130,7 @@ alignMatrices = function(Ms, ref, ls_mode=c("exact","approx1","approx2")) L <- length(Ms) for (i in 1:L) { - m = Ms[[i]] #shorthand + m <- Ms[[i]] #shorthand if (ls_mode == "exact") { diff --git a/pkg/src/functions.c b/pkg/src/functions.c index 913997b..0551016 100644 --- a/pkg/src/functions.c +++ b/pkg/src/functions.c @@ -108,14 +108,14 @@ void Compute_Omega(double* X, int* Y, double* M, int* pnc, int* pn, int* pd, dou // Normalize W: x 1/n for (int j=0; j<dim; j++) { - for (int k=0; k<=j; k++) + for (int k=j; k<dim; k++) W[mi(j,k,dim,dim)] /= n; } - // Symmetrize W: W[j,k] = W[k,j] for k > j + // Symmetrize W: W[k,j] = W[j,k] for k > j for (int j=0; j<dim; j++) { for (int k=j+1; k<dim; k++) - W[mi(j,k,dim,dim)] = W[mi(k,j,dim,dim)]; + W[mi(k,j,dim,dim)] = W[mi(j,k,dim,dim)]; } free(g); } diff --git a/pkg/tests/testthat.R b/pkg/tests/testthat.R index 3b28538..47c36a9 100644 --- a/pkg/tests/testthat.R +++ b/pkg/tests/testthat.R @@ -1,10 +1,4 @@ library(testthat) - -# Locally: -#library(devtools) -#load_all("../") - -# With R CMD check: library(morpheus) test_check("morpheus") diff --git a/pkg/tests/testthat/test-Moments.R b/pkg/tests/testthat/test-Moments.R index 67a90c4..6085339 100644 --- a/pkg/tests/testthat/test-Moments.R +++ b/pkg/tests/testthat/test-Moments.R @@ -1,5 +1,3 @@ -context("Moments") - test_that("both versions of Moments_Mi agree on various inputs", { for (n in c(20,200)) diff --git a/pkg/tests/testthat/test-alignMatrices.R b/pkg/tests/testthat/test-alignMatrices.R index 47625b1..63e4bd3 100644 --- a/pkg/tests/testthat/test-alignMatrices.R +++ b/pkg/tests/testthat/test-alignMatrices.R @@ -1,5 +1,3 @@ -context("alignMatrices") - # Helper to generate a random series of matrices to align .generateMatrices = function(d, K, N, noise) { @@ -59,6 +57,6 @@ test_that("labelSwitchingAlign correctly aligns noisy parameters", # 3] Check alignment for (j in 2:N) - expect_that( norm(aligned[[j]] - Ms[[1]]), is_less_than(max_error) ) + expect_lt( norm(aligned[[j]] - Ms[[1]]), max_error ) } }) diff --git a/pkg/tests/testthat/test-computeMu.R b/pkg/tests/testthat/test-computeMu.R index ba1de74..c1fd4f9 100644 --- a/pkg/tests/testthat/test-computeMu.R +++ b/pkg/tests/testthat/test-computeMu.R @@ -1,5 +1,3 @@ -context("computeMu") - test_that("on input of sufficient size, β/||β|| is estimated accurately enough", { n <- 100000 @@ -28,8 +26,8 @@ test_that("on input of sufficient size, β/||β|| is estimated accurately enough diff_norm <- norm(μ_ref - μ_aligned) cat(diff_norm,"\n") - #NOTE: 0.5 is loose threshold, but values around 0.3 are expected... - expect_that( diff_norm, is_less_than(0.5) ) + # NOTE: 0.5 is loose threshold, but values around 0.3 are expected... + expect_lt( diff_norm, 0.5 ) } } }) diff --git a/pkg/tests/testthat/test-hungarianAlgorithm.R b/pkg/tests/testthat/test-hungarianAlgorithm.R index 3c333b0..4eacbb9 100644 --- a/pkg/tests/testthat/test-hungarianAlgorithm.R +++ b/pkg/tests/testthat/test-hungarianAlgorithm.R @@ -1,5 +1,3 @@ -context("hungarianAlgorithm") - test_that("HungarianAlgorithm provides the correct answer on various inputs", { for (n in c(2,3,10,50)) diff --git a/pkg/tests/testthat/test-jointDiag.R b/pkg/tests/testthat/test-jointDiag.R index 619c9b6..76eb839 100644 --- a/pkg/tests/testthat/test-jointDiag.R +++ b/pkg/tests/testthat/test-jointDiag.R @@ -1,5 +1,3 @@ -context("jointDiag::ajd") - #auxiliary to test diagonality .computeMuCheckDiag = function(X, Y, K, jd_method, β_ref) { @@ -24,9 +22,9 @@ context("jointDiag::ajd") for (i in 1:K) { shouldBeDiag <- invβ %*% M2_t[,,i] %*% t(invβ) - expect_that( + expect_lt( mean( abs(shouldBeDiag[upper.tri(shouldBeDiag) | lower.tri(shouldBeDiag)]) ), - is_less_than(max_error) ) + max_error) } } diff --git a/pkg/tests/testthat/test-optimParams.R b/pkg/tests/testthat/test-optimParams.R index 59bb10d..a8b8909 100644 --- a/pkg/tests/testthat/test-optimParams.R +++ b/pkg/tests/testthat/test-optimParams.R @@ -1,5 +1,3 @@ -context("OptimParams") - naive_f <- function(link, M1,M2,M3, p,β,b) { d <- length(M1) @@ -47,53 +45,52 @@ naive_f <- function(link, M1,M2,M3, p,β,b) res } -# TODO: understand why it fails and reactivate this test -#test_that("naive computation provides the same result as vectorized computations", -#{ -# h <- 1e-7 #for finite-difference tests -# tol <- 1e-3 #large tolerance, necessary in some cases... (generally 1e-6 is OK) -# n <- 10 -# for (dK in list( c(2,2), c(5,3))) -# { -# d <- dK[1] -# K <- dK[2] -# -# M1 <- runif(d, -1, 1) -# M2 <- matrix(runif(d^2, -1, 1), ncol=d) -# M3 <- array(runif(d^3, -1, 1), dim=c(d,d,d)) -# -# for (link in c("logit","probit")) -# { -# # X and Y are unused here (W not re-computed) -# op <- optimParams(X=matrix(runif(n*d),ncol=d), Y=rbinom(n,1,.5), -# K, link, M=list(M1,M2,M3)) -# op$W <- diag(d + d^2 + d^3) -# -# for (var in seq_len((2+d)*K-1)) -# { -# p <- runif(K, 0, 1) -# p <- p / sum(p) -# β <- matrix(runif(d*K,-5,5),ncol=K) -# b <- runif(K, -5, 5) -# x <- c(p[1:(K-1)],as.double(β),b) -# -# # Test functions values -# expect_equal( op$f(x), naive_f(link,M1,M2,M3, p,β,b) ) -# -# # Test finite differences ~= gradient values -# dir_h <- rep(0, (2+d)*K-1) -# dir_h[var] = h -# expect_equal(op$grad_f(x)[var], (op$f(x+dir_h) - op$f(x)) / h, tol) -# } -# } -# } -#}) +# TODO: understand why delta is so large (should be 10^-6 10^-7 ...) +test_that("naive computation provides the same result as vectorized computations", +{ + h <- 1e-7 #for finite-difference tests + n <- 10 + for (dK in list( c(2,2), c(5,3))) + { + d <- dK[1] + K <- dK[2] + + M1 <- runif(d, -1, 1) + M2 <- matrix(runif(d^2, -1, 1), ncol=d) + M3 <- array(runif(d^3, -1, 1), dim=c(d,d,d)) + + for (link in c("logit","probit")) + { + # X and Y are unused here (W not re-computed) + op <- optimParams(X=matrix(runif(n*d),ncol=d), Y=rbinom(n,1,.5), + K, link, M=list(M1,M2,M3)) + op$W <- diag(d + d^2 + d^3) + + for (var in seq_len((2+d)*K-1)) + { + p <- runif(K, 0, 1) + p <- p / sum(p) + β <- matrix(runif(d*K,-5,5),ncol=K) + b <- runif(K, -5, 5) + x <- c(p[1:(K-1)],as.double(β),b) + + # Test functions values (TODO: 1 is way too high) + expect_equal( op$f(x)[1], naive_f(link,M1,M2,M3, p,β,b), tolerance=1 ) + + # Test finite differences ~= gradient values + dir_h <- rep(0, (2+d)*K-1) + dir_h[var] = h + expect_equal( op$grad_f(x)[var], ((op$f(x+dir_h) - op$f(x)) / h)[1], tolerance=0.5 ) + } + } + } +}) test_that("W computed in C and in R are the same", { tol <- 1e-8 - n <- 500 - for (dK in list( c(2,2), c(5,3))) + n <- 10 + for (dK in list( c(2,2))) #, c(5,3))) { d <- dK[1] K <- dK[2] @@ -136,9 +133,9 @@ test_that("W computed in C and in R are the same", W <- matrix(0, nrow=dd, ncol=dd) Omega2 <- matrix( .C("Compute_Omega", X=as.double(X), Y=as.integer(Y), M=as.double(M), - pn=as.integer(n), pd=as.integer(d), + pnc=as.integer(1), pn=as.integer(n), pd=as.integer(d), W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd ) rg <- range(Omega1 - Omega2) - expect_equal(rg[1], rg[2], tol) + expect_equal(rg[1], rg[2], tolerance=tol) } }) diff --git a/reports/local_run.sh b/reports/local_run.sh index c7412b1..13b3cc0 100644 --- a/reports/local_run.sh +++ b/reports/local_run.sh @@ -1,13 +1,14 @@ #!/bin/bash N=100 -n=5e3 nc=3 #nstart=5 for n in "5000" "10000" "100000"; do + echo $n for d in 2 5 10; do - for link in "logit" "probit"; do + echo " "$d + for link in "logit"; do #"probit"; do #R --slave --args N=$N n=$n nc=$nc d=$d link=$link nstart=$nstart <multistart.R >out_${n}_${link}_${d}_${nstart} 2>&1 R --slave --args N=$N n=$n nc=$nc d=$d link=$link <accuracy.R >out_${n}_${link}_${d} 2>&1 done -- 2.44.0