From: Benjamin Auder Date: Mon, 7 Jun 2021 15:27:23 +0000 (+0200) Subject: Adjustments + bugs fixing X-Git-Url: https://git.auder.net/doc/html/scripts/%7B%7B%20asset%28%27mixstore/css/DESCRIPTION?a=commitdiff_plain;h=ab35f6102896a49e86e853262c0650faa2931638;p=morpheus.git Adjustments + bugs fixing --- 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 j + // Symmetrize W: W[k,j] = W[j,k] for k > j for (int j=0; jout_${n}_${link}_${d}_${nstart} 2>&1 R --slave --args N=$N n=$n nc=$nc d=$d link=$link out_${n}_${link}_${d} 2>&1 done