From 99aaf41fa03e2c17a9bf49bcbd591850fd90b58d Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Mon, 8 Jul 2019 14:57:58 +0800 Subject: [PATCH 01/32] Optimize MoMA::criterion_search --- src/moma_level1.cpp | 73 ++++++++++++------- .../test_multirank-BIC-single_problem.R | 24 ++++++ 2 files changed, 71 insertions(+), 26 deletions(-) diff --git a/src/moma_level1.cpp b/src/moma_level1.cpp index 202545c0..0a12b92c 100644 --- a/src/moma_level1.cpp +++ b/src/moma_level1.cpp @@ -52,6 +52,11 @@ Rcpp::List MoMA::criterion_search(const arma::vec &bic_au_grid, double tol = 1; int iter = 0; + int n_au = bic_au_grid.n_elem; + int n_av = bic_av_grid.n_elem; + int n_lu = bic_lu_grid.n_elem; + int n_lv = bic_lv_grid.n_elem; + // u_/v_result is a list of the following contents // Rcpp::Named("lambda") = opt_lambda_u, the chosen lambda // Rcpp::Named("alpha") = opt_alpha_u, the chosen alpha @@ -72,29 +77,42 @@ Rcpp::List MoMA::criterion_search(const arma::vec &bic_au_grid, // We conduct 2 BIC searches over 2D grids here instead // of 4 searches over 1D grids. It's consistent with // Genevera's code. - while (tol > EPS_bic && iter < max_bic_iter) + if (n_au > 1 || n_av > 1 || n_lu > 1 || n_lv > 1) { - iter++; - oldu = curu; - oldv = curv; - - // choose lambda/alpha_u - MoMALogger::debug("Start u search."); - u_result = bicsr_u.search(X * curv, curu, bic_au_grid, bic_lu_grid); - curu = Rcpp::as(u_result["vector"]); - - MoMALogger::debug("Start v search."); - v_result = bicsr_v.search(X.t() * curu, curv, bic_av_grid, bic_lv_grid); - curv = Rcpp::as(v_result["vector"]); - - double scale_u = arma::norm(oldu) == 0.0 ? 1 : arma::norm(oldu); - double scale_v = arma::norm(oldv) == 0.0 ? 1 : arma::norm(oldv); - - tol = arma::norm(oldu - u) / scale_u + arma::norm(oldv - v) / scale_v; - MoMALogger::debug("Finish nested greedy BIC search outer loop. (iter, tol) = (") - << iter << "," << tol << "), " - << "(bic_u, bic_v) = (" << (double)u_result["bic"] << "," << (double)v_result["bic"] - << ")"; + while (tol > EPS_bic && iter < max_bic_iter) + { + iter++; + oldu = curu; + oldv = curv; + + // choose lambda/alpha_u + MoMALogger::debug("Start u search."); + u_result = bicsr_u.search(X * curv, curu, bic_au_grid, bic_lu_grid); + curu = Rcpp::as(u_result["vector"]); + + MoMALogger::debug("Start v search."); + v_result = bicsr_v.search(X.t() * curu, curv, bic_av_grid, bic_lv_grid); + curv = Rcpp::as(v_result["vector"]); + + double scale_u = arma::norm(oldu) == 0.0 ? 1 : arma::norm(oldu); + double scale_v = arma::norm(oldv) == 0.0 ? 1 : arma::norm(oldv); + + tol = arma::norm(oldu - u) / scale_u + arma::norm(oldv - v) / scale_v; + MoMALogger::debug("Finish nested greedy BIC search outer loop. (iter, tol) = (") + << iter << "," << tol << "), " + << "(bic_u, bic_v) = (" << (double)u_result["bic"] << "," << (double)v_result["bic"] + << ")"; + } + } + else + { + u_result = Rcpp::List::create( + Rcpp::Named("lambda") = bic_lu_grid(0), Rcpp::Named("alpha") = bic_au_grid(0), + Rcpp::Named("vector") = initial_u, Rcpp::Named("bic") = -MOMA_INFTY); + v_result = Rcpp::List::create( + Rcpp::Named("lambda") = bic_lv_grid(0), Rcpp::Named("alpha") = bic_av_grid(0), + Rcpp::Named("vector") = initial_v, Rcpp::Named("bic") = -MOMA_INFTY); + MoMALogger::debug("Deprecated BIC grid. Skip searching."); } // A final run on the selected parameter @@ -112,6 +130,8 @@ Rcpp::List MoMA::criterion_search(const arma::vec &bic_au_grid, u_result["vector"] = u; v_result["vector"] = v; + + // NOTE: we do not update bic for the new u and v. } return Rcpp::List::create(Rcpp::Named("u_result") = u_result, @@ -213,18 +233,19 @@ Rcpp::List MoMA::grid_BIC_mix(const arma::vec &alpha_u, Rcpp::List u_result = result["u_result"]; Rcpp::List v_result = result["v_result"]; + arma::vec curu = Rcpp::as(u_result["vector"]); + arma::vec curv = Rcpp::as(v_result["vector"]); + double d = arma::as_scalar(curu.t() * X * curv); + five_d_list.insert( Rcpp::List::create(Rcpp::Named("u") = u_result, Rcpp::Named("v") = v_result, Rcpp::Named("k") = pc, - Rcpp::Named("X") = X), + Rcpp::Named("X") = X, Rcpp::Named("d") = d), i, j, k, m, pc); // Deflate the matrix if (pc < rank - 1) { - arma::vec curu = Rcpp::as(u_result["vector"]); - arma::vec curv = Rcpp::as(v_result["vector"]); - double d = arma::as_scalar(curu.t() * X * curv); deflate(d); } } diff --git a/tests/testthat/test_multirank-BIC-single_problem.R b/tests/testthat/test_multirank-BIC-single_problem.R index 86cf7819..69cca9f9 100644 --- a/tests/testthat/test_multirank-BIC-single_problem.R +++ b/tests/testthat/test_multirank-BIC-single_problem.R @@ -80,3 +80,27 @@ test_that("`cpp_multirank_BIC_grid_search` solves a single MoMA problem", { # all zeros expect_lte(zero_cnt / (n_lambda_u * n_lambda_v * n_alpha_u * n_alpha_v), 0.1) }) + +test_that("`cpp_multirank_BIC_grid_search`: a naive case", { + arglist_x_w_penalty <- modifyList( + public_arglist_wo_penalty, + list( + X = matrix(1), + Omega_u = matrix(0), + Omega_v = matrix(0), + alpha_u = 0, + alpha_v = 0, + lambda_u = 0, + lambda_v = 0 + ) + ) + + result <- do.call( + cpp_multirank_BIC_grid_search, + arglist_x_w_penalty + )[[1]] + + expect_equal(result$u$vector, as.matrix(1)) + expect_equal(result$v$vector, as.matrix(1)) + expect_equal(result$d, 1) +}) From e2c2e61485732e941e762d50bb445f7edf73188f Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Tue, 9 Jul 2019 23:09:21 +0800 Subject: [PATCH 02/32] R6 SFPCA object --- DESCRIPTION | 2 +- R/moma_sfpca.R | 262 +++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 263 insertions(+), 1 deletion(-) create mode 100644 R/moma_sfpca.R diff --git a/DESCRIPTION b/DESCRIPTION index ab94baec..10f87d3e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -11,7 +11,7 @@ Authors@R: c( Description: Unified approach to modern multivariate analysis providing sparse, smooth, and structured versions of PCA, PLS, LDA, and CCA. License: GPL (>= 2) -Imports: Rcpp +Imports: Rcpp, R6 Suggests: knitr, rmarkdown, diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R new file mode 100644 index 00000000..d14f4487 --- /dev/null +++ b/R/moma_sfpca.R @@ -0,0 +1,262 @@ +SFPCA <- R6::R6Class("SFPCA", list( + center = NULL, + scale = NULL, + grid_result = NULL, + Omega_u = NULL, + Omega_v = NULL, + u_sparsity = NULL, + v_sparsity = NULL, + rank = NULL, + alpha_u = NULL, + alpha_v = NULL, + lambda_u = NULL, + lambda_v = NULL, + selection_scheme_list = NULL, + n = NULL, + p = NULL, + X = NULL, + initialize = function(X, ..., + center = TRUE, scale = FALSE, + u_sparsity = empty(), v_sparsity = empty(), lambda_u = 0, lambda_v = 0, # lambda_u/_v is a vector or scalar + Omega_u = NULL, Omega_v = NULL, alpha_u = 0, alpha_v = 0, # so is alpha_u/_v + pg_setting = moma_pg_settings(), + selection_scheme_str = "gggg", + max_bic_iter = 5, + rank = 1) { + + # Step 1: check ALL arguments + # Step 1.1: lambdas and alphas + if (!inherits(alpha_u, c("numeric", "integer")) || + !inherits(alpha_v, c("numeric", "integer")) || + !inherits(lambda_u, c("numeric", "integer")) || + !inherits(lambda_v, c("numeric", "integer"))) { + moma_error(paste0( + "All penalty levels (", + sQuote("lambda_u"), ", ", + sQuote("lambda_v"), ", ", + sQuote("alpha_u"), ", ", + sQuote("alpha_v"), + ") must be numeric." + )) + } + self$alpha_u <- alpha_u + self$alpha_v <- alpha_v + self$lambda_u <- lambda_u + self$lambda_v <- lambda_v + + # Step 1.2: matrix + X <- as.matrix(X) + if (any(!is.finite(X))) { + moma_error("X must not have NaN, NA, or Inf.") + } + n <- dim(X)[1] + p <- dim(X)[2] + X <- scale(X, center = center, scale = scale) + + cen <- attr(X, "scaled:center") + sc <- attr(X, "scaled:scale") + if (any(sc == 0)) { + moma_error("cannot rescale a constant/zero column to unit variance") + } + + self$center <- if (is.null(cen)) FALSE else cen + self$scale <- if (is.null(sc)) FALSE else sc + self$n <- n + self$p <- p + self$X <- X + + # Step 1.3: sparsity + if (!inherits(u_sparsity, "moma_sparsity") || !inherits(v_sparsity, "moma_sparsity")) { + moma_error( + "Sparse penalty should be of class ", + sQuote("moma_sparsity"), + ". Try using, for example, `u_sparsity = lasso()`." + ) + } + self$u_sparsity <- u_sparsity + self$v_sparsity <- v_sparsity + + # Step 1.4: PG loop settings + if (!inherits(pg_setting, "moma_pg_settings")) { + moma_error( + "pg_setting penalty should be of class ", + sQuote("moma_pg_settings"), + ". Try using, for example, `pg_setting = moma_pg_settings(MAX_ITER=1e+4)`." + ) + } + + # Step 1.5: smoothness + Omega_u <- check_omega(Omega_u, alpha_u, n) + Omega_v <- check_omega(Omega_v, alpha_v, p) + self$Omega_u <- Omega_u + self$Omega_v <- Omega_v + + # Step 1.6: check selection scheme string + # "g" stands for grid search, "b" stands for BIC + if (!inherits(selection_scheme_str, "character") || nchar(selection_scheme_str) != 4 || + !all(strsplit(selection_scheme_str, split = "")[[1]] %in% c("b", "g"))) { + moma_error( + "Invalid selection_scheme_str ", selection_scheme_str, + ". It should be a four-char string containing only 'b' or 'g'." + ) + } + + # turn "b"/"g" to 1/0 + selection_scheme_list <- list( + selection_criterion_alpha_u = 0, + selection_criterion_alpha_v = 0, + selection_criterion_lambda_u = 0, + selection_criterion_lambda_v = 0 + ) + for (i in 1:4) { + selection_scheme_list[[i]] <- ifelse(strsplit(selection_scheme_str, split = "")[[1]][i] == "g", 0, 1) + } + self$selection_scheme_list <- selection_scheme_list + + # Step 1.7: check rank + if (!inherits(rank, "numeric") || rank <= 0 + || rank > min(p, n)) { + moma_error("rank should be a legit positive integer.") + } + self$rank <- rank + + # Step 2: pack all arguments in a list + algo_settings_list <- c( + list( + X = X, + lambda_u = lambda_u, + lambda_v = lambda_v, + # smoothness + alpha_u = alpha_u, + alpha_v = alpha_v, + rank = rank + ), + list( + Omega_u = check_omega(Omega_u, alpha_u, n), + Omega_v = check_omega(Omega_v, alpha_v, p), + prox_arg_list_u = add_default_prox_args(u_sparsity), + prox_arg_list_v = add_default_prox_args(v_sparsity) + ), + pg_setting, + selection_scheme_list, + list( + max_bic_iter = max_bic_iter + ) + ) + # make sure we explicitly specify ALL arguments + if (length(algo_settings_list) != length(formals(cpp_multirank_BIC_grid_search))) { + moma_error("Incomplete arguments in SFPCA::initialize.") + } + + # Step 3: call the fucntion + self$grid_result <- do.call( + cpp_multirank_BIC_grid_search, + algo_settings_list + ) + }, + + get_mat = function(alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) { + # Sanity check: if a parameter has been chosen by BIC, then + # the index for that parameter should not be specified. + parameter <- c(alpha_u, alpha_v, lambda_u, lambda_v) + if (any(parameter != 1 & self$selection_scheme_list != 0)) { + moma_error("Invalid index in SFPCA::get_mat. Do not specify indexes of parameters chosen by BIC.") + } + + + n <- self$n + p <- self$p + rank <- self$rank + + U <- matrix(0, nrow = n, ncol = rank) + V <- matrix(0, nrow = p, ncol = rank) + for (i in (1:self$rank)) { + U[, i] <- + get_5Dlist_elem(self$grid_result, + alpha_u_i = alpha_u, + lambda_u_i = lambda_u, + alpha_v_i = alpha_v, + lambda_v_i = lambda_v, rank_i = i + )[[1]]$u$vector + + V[, i] <- + get_5Dlist_elem(self$grid_result, + alpha_u_i = alpha_u, + lambda_u_i = lambda_u, + alpha_v_i = alpha_v, + lambda_v_i = lambda_v, rank_i = i + )[[1]]$v$vector + } + + coln <- if (is.null(colnames(self$X))) paste0("Xcol_", seq_len(p)) else colnames(self$X) + rown <- if (is.null(rownames(self$X))) paste0("Xrow_", seq_len(n)) else rownames(self$X) + dimnames(V) <- + list(coln, paste0("PC", seq_len(rank))) + dimnames(U) <- + list(rown, paste0("PC", seq_len(rank))) + return(list(U = U, V = V)) + }, + + print = function() { + selection_list_str <- lapply(self$selection_scheme_list, function(x) { + if (x == 0) { + return("grid search") + } + else if (x == 1) { + return("BIC search") + } + }) + + cat("An object containing solutions to the following settings\n") + cat("rank =", self$rank, "\n") + cat("Penalty and selection:\n") + + cat(paste0("alpha_u: ", selection_list_str[1]), "\n") + print(self$alpha_u) + cat(paste0("alpha_u: ", selection_list_str[2]), "\n") + print(self$alpha_v) + cat(paste0("lambda_u: ", selection_list_str[3]), "\n") + print(self$lambda_u) + cat(paste0("lambda_v: ", selection_list_str[4]), "\n") + print(self$lambda_v) + }, + + left_project = function(newX, ..., + alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) { + + # check indexes + parameter <- c(alpha_u, alpha_v, lambda_u, lambda_v) + if (any(parameter != 1 & self$selection_scheme_list != 0)) { + moma_error("Invalid index in SFPCA::get_mat. Do not specify indexes of parameters chosen by BIC.") + } + V <- self$get_mat( + alpha_u = alpha_u, + alpha_v = alpha_v, + lambda_u = lambda_u, + lambda_v = lambda_v + )$V + + + # newX should be uncencter and unscaled. + # check new X has same colnames + if (length(dim(newX)) != 2L) { + moma_error("'newX' must be a matrix or data frame") + } + + if (dim(newX)[2] != self$p) { + moma_error("`newX` is incompatible with orignal data.") + } + + PV <- solve((t(V) %*% V), t(V)) + scaled_data <- scale(newX, self$center, self$scale) + result <- scaled_data %*% t(PV) + colnames(result) <- paste0("PC", seq_len(self$rank)) + return(list( + scaled_data = scaled_data, + proj_data = result, + V = V, + PV = PV + )) + } +)) + From 789e86aa38282ef5dfd321240252dedc2b3fb4be Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Tue, 9 Jul 2019 23:09:56 +0800 Subject: [PATCH 03/32] Wrappers around SFPCA object, and special-case functions --- R/moma_sfpca.R | 205 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 205 insertions(+) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index d14f4487..a5a04f8f 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -260,3 +260,208 @@ SFPCA <- R6::R6Class("SFPCA", list( } )) +moma_sfpca <- function(X, ..., + center = TRUE, scale = FALSE, + u_sparsity = empty(), v_sparsity = empty(), lambda_u = 0, lambda_v = 0, # lambda_u/_v is a vector or scalar + Omega_u = NULL, Omega_v = NULL, alpha_u = 0, alpha_v = 0, # so is alpha_u/_v + pg_setting = moma_pg_settings(), + selection_scheme_str = "gggg", + max_bic_iter = 5, + rank = 1) { + chkDots(...) + return(SFPCA$new( + X, + center = center, scale = scale, + # sparsity + u_sparsity = u_sparsity, v_sparsity = v_sparsity, + lambda_u = lambda_u, lambda_v = lambda_v, + # smoothness + Omega_u = Omega_u, Omega_v = Omega_v, + alpha_u = alpha_u, alpha_v = alpha_v, + pg_setting = pg_setting, + selection_scheme_str = selection_scheme_str, + max_bic_iter = max_bic_iter, + rank = rank + )) +} + +moma_spca <- function(X, ..., + center = TRUE, scale = FALSE, + u_sparsity = empty(), v_sparsity = empty(), lambda_u = 0, lambda_v = 0, # lambda_u/_v is a vector or scalar + pg_setting = moma_pg_settings(), + selection_scheme_str = "g", + max_bic_iter = 5, + rank = 1) { + chkDots(...) + u_penalized <- !(missing(u_sparsity) && missing(lambda_u)) + v_penalized <- !(missing(v_sparsity) && missing(lambda_v)) + if (!u_penalized && !v_penalized) { + moma_warning("No sparsity is imposed!") + } + + if (u_penalized && v_penalized) { + moma_error("Please use `moma_twspca` if both sides are penalized.") + } + + if (nchar(selection_scheme_str) != 1 || !selection_scheme_str %in% c("g", "b")) { + moma_error("`selection_scheme_str` should be either 'g' or 'b'") + } + + # selection_scheme_str is in order of alpha_u/v, lambda_u/v + full_selection_scheme_str <- + if (u_penalized) { + paste0("gg", selection_scheme_str, "g") + } + else if (v_penalized) { + paste0("ggg", selection_scheme_str) + } + else { + "gggg" + } + + return(moma_sfpca( + X, + center = center, scale = scale, + u_sparsity = u_sparsity, v_sparsity = v_sparsity, + lambda_u = lambda_u, lambda_v = lambda_v, + # Omega_u = Omega_u, Omega_v = Omega_v, + # alpha_u = alpha_u, alpha_v = alpha_v, + pg_setting = pg_setting, + selection_scheme_str = full_selection_scheme_str, + max_bic_iter = max_bic_iter, + rank = rank + )) + # moma_error("Not implemented: SPCA") +} + +moma_twspca <- function(X, ..., + center = TRUE, scale = FALSE, + u_sparsity = empty(), v_sparsity = empty(), lambda_u = 0, lambda_v = 0, # lambda_u/_v is a vector or scalar + pg_setting = moma_pg_settings(), + selection_scheme_str = "gg", + max_bic_iter = 5, + rank = 1) { + chkDots(...) + u_penalized <- !(missing(u_sparsity) && missing(lambda_u)) + v_penalized <- !(missing(v_sparsity) && missing(lambda_v)) + if (!u_penalized && !v_penalized) { + moma_warning("No sparsity is imposed!") + } + + if (!u_penalized || !v_penalized) { + moma_warning("Please use `moma_spca` if only one side is penalized.") + } + + if (nchar(selection_scheme_str) != 2) { + moma_error("`selection_scheme_str` should be of length two.") + } + if (!all(strsplit(selection_scheme_str, split = "")[[1]] %in% c("b", "g"))) { + moma_error("`selection_scheme_str` should consist of 'g' or 'b'") + } + + # selection_scheme_str is in order of alpha_u/v, lambda_u/v + full_selection_scheme_str <- paste0("gg", selection_scheme_str) + + return(moma_sfpca( + X, + center = center, scale = scale, + u_sparsity = u_sparsity, v_sparsity = v_sparsity, + lambda_u = lambda_u, lambda_v = lambda_v, + # Omega_u = Omega_u, Omega_v = Omega_v, + # alpha_u = alpha_u, alpha_v = alpha_v, + pg_setting = pg_setting, + selection_scheme_str = full_selection_scheme_str, + max_bic_iter = max_bic_iter, + rank = rank + )) +} + +moma_fpca <- function(X, ..., + center = TRUE, scale = FALSE, + Omega_u = NULL, Omega_v = NULL, alpha_u = 0, alpha_v = 0, + pg_setting = moma_pg_settings(), + selection_scheme_str = "g", + max_bic_iter = 5, + rank = 1) { + chkDots(...) + u_penalized <- !(missing(Omega_u) && missing(alpha_u)) + v_penalized <- !(missing(Omega_v) && missing(alpha_v)) + if (!u_penalized && !v_penalized) { + moma_warning("No smoothness is imposed!") + } + + if (u_penalized && v_penalized) { + moma_error("Please use `moma_twfpca` if both sides are penalized.") + } + + if (nchar(selection_scheme_str) != 1 || !selection_scheme_str %in% c("g", "b")) { + moma_error("`selection_scheme_str` should be either 'g' or 'b'") + } + + # selection_scheme_str is in order of alpha_u/v, lambda_u/v + full_selection_scheme_str <- + if (u_penalized) { + paste0(selection_scheme_str, "ggg") + } + else if (v_penalized) { + paste0("g", selection_scheme_str, "gg") + } + else { + "gggg" + } + + return(moma_sfpca( + X, + center = center, scale = scale, + # u_sparsity = u_sparsity, v_sparsity = v_sparsity, + # lambda_u = lambda_u, lambda_v = lambda_v, + Omega_u = Omega_u, Omega_v = Omega_v, + alpha_u = alpha_u, alpha_v = alpha_v, + pg_setting = pg_setting, + selection_scheme_str = full_selection_scheme_str, + max_bic_iter = max_bic_iter, + rank = rank + )) +} + +moma_twfpca <- function(X, ..., + center = TRUE, scale = FALSE, + Omega_u = NULL, Omega_v = NULL, alpha_u = 0, alpha_v = 0, + pg_setting = moma_pg_settings(), + selection_scheme_str = "gg", + max_bic_iter = 5, + rank = 1) { + chkDots(...) + u_penalized <- !(missing(Omega_u) && missing(alpha_u)) + v_penalized <- !(missing(Omega_v) && missing(alpha_v)) + if (!u_penalized && !v_penalized) { + moma_warning("No smoothness is imposed!") + } + + if (!u_penalized || !v_penalized) { + moma_warning("Please use `moma_fpca` if only one side is penalized.") + } + + if (nchar(selection_scheme_str) != 2) { + moma_error("`selection_scheme_str` should be of length two.") + } + if (!all(strsplit(selection_scheme_str, split = "")[[1]] %in% c("b", "g"))) { + moma_error("`selection_scheme_str` should consist of 'g' or 'b'") + } + + # selection_scheme_str is in order of alpha_u/v, lambda_u/v + full_selection_scheme_str <- paste0(selection_scheme_str, "gg") + + return(moma_sfpca( + X, + center = center, scale = scale, + # u_sparsity = u_sparsity, v_sparsity = v_sparsity, + # lambda_u = lambda_u, lambda_v = lambda_v, + Omega_u = Omega_u, Omega_v = Omega_v, + alpha_u = alpha_u, alpha_v = alpha_v, + pg_setting = pg_setting, + selection_scheme_str = full_selection_scheme_str, + max_bic_iter = max_bic_iter, + rank = rank + )) +} From 64296c55b4ae6de776d8c7627079197e4073fbbd Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Tue, 9 Jul 2019 23:10:54 +0800 Subject: [PATCH 04/32] Optimize `check_omega` --- R/moma_expose.R | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/R/moma_expose.R b/R/moma_expose.R index 7ec8f371..ff3ffd5c 100644 --- a/R/moma_expose.R +++ b/R/moma_expose.R @@ -31,6 +31,13 @@ add_default_prox_args <- function(sparsity_type) { # This function checks the validity of Omega and alpha check_omega <- function(Omega, alpha, n) { + + # check if Omega is a matrix + if (!is.matrix(Omega) && !is.null(Omega)) { + moma_error("Omage_u/v is not a matrix.") + } + + # store them as sparse matrices using the package Matrix if (length(alpha) == 1 && alpha == 0) { # discard the Omega matrix specified by users Omega <- diag(n) @@ -41,6 +48,9 @@ check_omega <- function(Omega, alpha, n) { Omega <- second_diff_mat(n) } else { + # At this point, users have specified an Omega and + # non-zero penalty levels explicitly + # Check validity of Omega if users speicify both alpha and Omega if (dim(Omega)[1] != dim(Omega)[2]) { moma_error( From a71c165a440ac75f8c153d91e1ddb58ec2439ae9 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Tue, 9 Jul 2019 23:12:03 +0800 Subject: [PATCH 05/32] Tests --- tests/testthat/test_sfpca_wrapper.R | 537 ++++++++++++++++++++++++++++ 1 file changed, 537 insertions(+) create mode 100644 tests/testthat/test_sfpca_wrapper.R diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R new file mode 100644 index 00000000..b54468f1 --- /dev/null +++ b/tests/testthat/test_sfpca_wrapper.R @@ -0,0 +1,537 @@ +context("Test R6 object") +set.seed(12) +test_that("SFPCA object: a naive case, 1x1 matrix", { + a <- SFPCA$new(matrix(1), center = FALSE) + result_to_be_tested <- a$grid_result[[1]] + + expect_equal(result_to_be_tested$X, matrix(1)) + + expect_equal(result_to_be_tested$u$bic, -Inf) + expect_equal(result_to_be_tested$u$lambda, 0) + expect_equal(result_to_be_tested$u$alpha, 0) + expect_equal(result_to_be_tested$u$vector, matrix(1)) + + expect_equal(result_to_be_tested$v$bic, -Inf) + expect_equal(result_to_be_tested$v$alpha, 0) + expect_equal(result_to_be_tested$v$lambda, 0) + expect_equal(result_to_be_tested$v$vector, matrix(1)) + + expect_equal(result_to_be_tested$k, 0) + expect_equal(result_to_be_tested$d, 1) +}) + +test_that("SFPCA object: correct arguments", { + a <- SFPCA$new(matrix(runif(12), 3, 4), scale = FALSE, cen = FALSE) + expect_equal(a$Omega_u, diag(3)) + expect_equal(a$Omega_v, diag(4)) + expect_equal(a$sc, NULL) + + a <- SFPCA$new(matrix(runif(12), 3, 4), alpha_u = 1) + expect_equal(a$Omega_u, second_diff_mat(3)) + expect_equal(a$Omega_v, diag(4)) + + a <- SFPCA$new(matrix(runif(12), 3, 4), u_sparsity = lasso()) + expect_equal(a$u_sparsity, lasso()) + expect_equal(a$v_sparsity, empty()) + + expect_error( + SFPCA$new(matrix(runif(12), 3, 4), selection_scheme_str = "bbba"), + "Invalid selection_scheme_str bbba. It should be a four-char string containing only 'b' or 'g'. " + ) + + a <- SFPCA$new(matrix(runif(12), 3, 4), selection_scheme_str = "ggbb") # in order of alpha_u/v, lambda_u/v + + expect_true(all(a$selection_scheme_list == c(0, 0, 1, 1))) + + expect_equal(dim(a$grid_result), c(1, 1, 1, 1, 1)) +}) + + + +test_that("SFPCA object: as SVD", { + + # test no penalty case + X <- matrix(runif(12), 4, 3) + a <- SFPCA$new(X, rank = 3, center = FALSE, scale = FALSE) + + + mysvd <- a$get_mat(1, 1, 1, 1) + svda <- svd(X, nu = 3, nv = 3) + + expect_equal(abs(mysvd$U), abs(svda$u), check.attributes = FALSE) # use `abs` to avoid sign difference + expect_equal(abs(mysvd$V), abs(svda$v), check.attributes = FALSE) + + + # test a matrix with dimnames + rown <- paste0("row", seq(1:4)) + coln <- paste0("col", seq(1:3)) + dimnames(X) <- list(rown, coln) + + a <- SFPCA$new(X, rank = 3, center = FALSE, scale = FALSE) + mysvd <- a$get_mat(1, 1, 1, 1) + + expect_equal(rownames(mysvd$U), rown) + expect_equal(rownames(mysvd$V), coln) + + + # test default arguments + expect_equal(a$get_mat(), a$get_mat(1, 1, 1, 1)) + + # test selection with BIC seach and grid search + a <- SFPCA$new(X, + rank = 3, center = FALSE, scale = FALSE, + alpha_u = seq(0, 2, 0.2), selection_scheme_str = "bggg" + ) + expect_true(all(a$selection_scheme_list == c(1, 0, 0, 0))) + + expect_equal(dim(a$grid_result), c(1, 1, 1, 1, 3)) + expect_error( + a$get_mat(2, 1, 1, 1), + "Invalid index in SFPCA::get_mat. Do not specify indexes of parameters chosen by BIC." + ) +}) + +test_that("SFPCA object: print fucntion", { + X <- matrix(runif(12), 4, 3) + a <- SFPCA$new(X, + rank = 3, + alpha_u = c(1, 2), + lambda_u = c(2, 3), + alpha_v = c(3, 4), + lambda_v = c(6, 7), + selection_scheme_str = "bgbg" + ) + print_message <- capture.output(print(a)) + + expect_message <- c( + "An object containing solutions to the following settings", + "rank = 3 ", + "Penalty and selection:", + "alpha_u: BIC search ", + "[1] 1 2", + "alpha_u: grid search ", + "[1] 3 4", + "lambda_u: BIC search ", + "[1] 2 3", + "lambda_v: grid search ", + "[1] 6 7" + ) + expect_equal(expect_message, print_message) +}) + +test_that("SFPCA object: left-project fucntion", { + + # unamed matrix + X <- matrix(runif(17 * 8), 17, 8) + a <- SFPCA$new(X, + rank = 3, + alpha_u = c(1), + lambda_u = c(2), + alpha_v = c(3), + lambda_v = c(6), + selection_scheme_str = "bbbb" + ) + expect_error( + a$left_project(matrix(0, 4, 1)), + "`newX` is incompatible with orignal data." + ) + + new_data <- matrix(runif(24), 3, 8) + res <- a$left_project(new_data) + + V <- res$V + # we verify the the projected data + # satisfies the normal equation: + # X^T X b = X^T y + expect_equal( + t(V) %*% V %*% t(res$proj_data), + t(V) %*% t(res$scaled_data) + ) +}) + + + +test_that("SFPCA object: moma_spca", { + X <- matrix(runif(17 * 8), 17, 8) + + # test inputs + expect_error( + moma_spca(X, lambda_v = c(1, 2), lambda_u = c(2, 3)), + "Please use `moma_twspca` if both sides are penalized" + ) + expect_error( + moma_spca(X, lambda_v = c(1, 2), u_sparsity = empty()), + "Please use `moma_twspca` if both sides are penalized" + ) + expect_error( + moma_spca(X, v_sparsity = empty(), u_sparsity = empty()), + "Please use `moma_twspca` if both sides are penalized" + ) + expect_error( + moma_spca(X, + v_sparsity = empty(), lambda_v = c(1, 2), + lambda_u = c(0) + ), + "Please use `moma_twspca` if both sides are penalized" + ) + expect_error( + moma_spca(X, lambda_v = c(1, 2), lambda_u = c(0)), + "Please use `moma_twspca` if both sides are penalized" + ) + expect_warning(moma_spca(X), "No sparsity is imposed!") + + + # test selection schemes + expect_error( + moma_spca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + selection_scheme_str = "gb" + ), + "`selection_scheme_str` should be either 'g' or 'b'" + ) + + expect_error( + moma_spca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + selection_scheme_str = "c" + ), + "`selection_scheme_str` should be either 'g' or 'b'" + ) + + + a <- moma_spca(X) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + + a <- moma_spca(X, selection_scheme_str = "b") # no effects + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + + a <- moma_spca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + selection_scheme_str = "g" + ) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + + a <- moma_spca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + selection_scheme_str = "b" + ) + expect_true(all(a$selection_scheme_list == c(0, 0, 1, 0))) + + + a <- moma_spca(X, + lambda_v = seq(0, 2, 0.2), v_sparsity = lasso(), + selection_scheme_str = "g" + ) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + + a <- moma_spca(X, + lambda_v = seq(0, 2, 0.2), v_sparsity = lasso(), + selection_scheme_str = "b" + ) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 1))) +}) + + +test_that("SFPCA object: moma_twspca", { + X <- matrix(runif(17 * 8), 17, 8) + + # test inputs + expect_no_error( + moma_twspca(X, lambda_v = c(1, 2), lambda_u = c(2, 3)) + ) + expect_no_error( + moma_twspca(X, lambda_v = c(1, 2), u_sparsity = empty()) + ) + expect_no_error( + moma_twspca(X, v_sparsity = empty(), u_sparsity = empty()) + ) + + expect_no_error( + moma_twspca(X, + v_sparsity = empty(), lambda_v = c(1, 2), + lambda_u = c(0) + ) + ) + + expect_warning( + moma_twspca(X), + "No sparsity is imposed!" + ) + expect_warning( + moma_twspca(X, lambda_u = seq(0, 2, 0.2), u_sparsity = lasso()), + "Please use `moma_spca` if only one side is penalized" + ) + + + # test selection schemes + expect_error( + moma_twspca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + selection_scheme_str = "gbb" + ), + "`selection_scheme_str` should be of length two." + ) + + expect_error( + moma_twspca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + selection_scheme_str = "cc" + ), + "`selection_scheme_str` should consist of 'g' or 'b'" + ) + + expect_no_error( + moma_twspca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + selection_scheme_str = "bb" + ) + ) + + + expect_warning(a <- moma_twspca(X)) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + + expect_warning(a <- moma_twspca(X, selection_scheme_str = "bb")) + expect_true(all(a$selection_scheme_list == c(0, 0, 1, 1))) + expect_warning(a <- moma_twspca(X, selection_scheme_str = "gg")) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + expect_warning(a <- moma_twspca(X, selection_scheme_str = "bg")) + expect_true(all(a$selection_scheme_list == c(0, 0, 1, 0))) + expect_warning(a <- moma_twspca(X, selection_scheme_str = "gb")) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 1))) + + + a <- moma_twspca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + lambda_v = seq(0, 2, 0.2), v_sparsity = lasso(), + selection_scheme_str = "bb" + ) + expect_true(all(a$selection_scheme_list == c(0, 0, 1, 1))) + + a <- moma_twspca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + lambda_v = seq(0, 2, 0.2), v_sparsity = lasso(), + selection_scheme_str = "gg" + ) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + + a <- moma_twspca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + lambda_v = seq(0, 2, 0.2), v_sparsity = lasso(), + selection_scheme_str = "gb" + ) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 1))) + + a <- moma_twspca(X, + lambda_v = seq(0, 2, 0.2), v_sparsity = lasso(), + selection_scheme_str = "bg" + ) + expect_true(all(a$selection_scheme_list == c(0, 0, 1, 0))) +}) + + +test_that("SFPCA object: moma_fpca", { + X <- matrix(runif(17 * 8), 17, 8) + + # test inputs + expect_error( + moma_fpca(X, alpha_v = c(1, 2), alpha_u = c(2, 3)), + "Please use `moma_twfpca` if both sides are penalized" + ) + expect_error( + moma_fpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = lasso(), + selection_scheme_str = "g" + ), + "Omage_u/v is not a matrix." + ) + + # test when no penalty matrix is provided + a <- moma_fpca(X, alpha_v = c(1, 2)) + expect_equal(a$Omega_u, diag(17)) # Omega is set to identiy mat if u or v is unpenalzied + expect_equal(a$Omega_v, second_diff_mat(8)) # the default penalty matrix + + expect_warning( + a <- moma_fpca(X), + "No smoothness is imposed!" + ) + expect_equal(a$Omega_u, diag(17)) + expect_equal(a$Omega_v, diag(8)) + + expect_error( + moma_fpca(X, + Omega_v = 2 * diag(8), alpha_v = c(1, 2), + alpha_u = c(0) + ), + "Please use `moma_twfpca` if both sides are penalized" + ) + expect_error( + moma_fpca(X, alpha_v = c(1, 2), alpha_u = c(0)), + "Please use `moma_twfpca` if both sides are penalized" + ) + expect_warning(moma_fpca(X), "No smoothness is imposed!") + + + # test selection schemes + # error when nchar(selection_scheme_str) != 1 + expect_error( + moma_fpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = second_diff_mat(17), + selection_scheme_str = "gb" + ), + "`selection_scheme_str` should be either 'g' or 'b'" + ) + + expect_error( + moma_fpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = second_diff_mat(8), + selection_scheme_str = "c" + ), + "`selection_scheme_str` should be either 'g' or 'b'" + ) + + + a <- moma_fpca(X) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + + a <- moma_fpca(X, selection_scheme_str = "b") # no effects + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + + a <- moma_fpca(X, + alpha_u = seq(0, 2, 0.2), + selection_scheme_str = "g" + ) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + + a <- moma_fpca(X, + alpha_u = seq(0, 2, 0.2), + selection_scheme_str = "b" + ) + expect_true(all(a$selection_scheme_list == c(1, 0, 0, 0))) + expect_equal(a$Omega_u, second_diff_mat(17)) + + a <- moma_fpca(X, + alpha_v = seq(0, 2, 0.2), + selection_scheme_str = "g" + ) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + expect_equal(a$Omega_v, second_diff_mat(8)) + + a <- moma_fpca(X, + alpha_v = seq(0, 2, 0.2), + selection_scheme_str = "b" + ) + expect_true(all(a$selection_scheme_list == c(0, 1, 0, 0))) +}) + + +test_that("SFPCA object: moma_twfpca", { + X <- matrix(runif(17 * 8), 17, 8) + + # test inputs + expect_no_error( + moma_twfpca(X, alpha_v = c(1, 2), alpha_u = c(2, 3)) + ) + expect_no_error( + moma_twfpca(X, alpha_v = c(1, 2)) + ) + expect_no_error( + moma_twfpca(X, Omega_v = diag(8), Omega_u = 2.1 * second_diff_mat(17)) + ) + + # incompatible Omega + expect_error( + moma_twfpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = 2.1 * second_diff_mat(11), + selection_scheme_str = "bb" + ), + "Omega shoud be a compatible matrix. It should be of 17x17, but is actually 11x11" + ) + + expect_no_error( + moma_twfpca(X, + Omega_v = diag(8), alpha_v = c(1, 2), + alpha_u = c(0) + ) + ) + + expect_warning( + moma_twfpca(X), + "No smoothness is imposed!" + ) + expect_warning( + moma_twfpca(X, alpha_u = seq(0, 2, 0.2), Omega_u = 2.3 * second_diff_mat(17)), + "Please use `moma_fpca` if only one side is penalized" + ) + + + # test selection schemes + expect_error( + moma_twfpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = lasso(), + selection_scheme_str = "gbb" + ), + "`selection_scheme_str` should be of length two." + ) + + expect_error( + moma_twfpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = lasso(), + selection_scheme_str = "cc" + ), + "`selection_scheme_str` should consist of 'g' or 'b'" + ) + + expect_no_error( + moma_twfpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = 2.1 * second_diff_mat(17), + selection_scheme_str = "bb" + ) + ) + + + + expect_warning(a <- moma_twfpca(X)) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + + expect_warning(a <- moma_twfpca(X, selection_scheme_str = "bb")) + expect_true(all(a$selection_scheme_list == c(1, 1, 0, 0))) + expect_warning(a <- moma_twfpca(X, selection_scheme_str = "gg")) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + expect_warning(a <- moma_twfpca(X, selection_scheme_str = "bg")) + expect_true(all(a$selection_scheme_list == c(1, 0, 0, 0))) + expect_warning(a <- moma_twfpca(X, selection_scheme_str = "gb")) + expect_true(all(a$selection_scheme_list == c(0, 1, 0, 0))) + + + a <- moma_twfpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = 1.1 * second_diff_mat(17), + alpha_v = seq(0, 2, 0.2), Omega_v = 1.2 * second_diff_mat(8), + selection_scheme_str = "bb" + ) + expect_true(all(a$selection_scheme_list == c(1, 1, 0, 0))) + expect_equal(a$Omega_u, 1.1 * second_diff_mat(17)) + expect_equal(a$Omega_v, 1.2 * second_diff_mat(8)) + + a <- moma_twfpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = 1.1 * second_diff_mat(17), + alpha_v = seq(0, 2, 0.2), Omega_v = 1.2 * second_diff_mat(8), + selection_scheme_str = "gg" + ) + expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) + expect_equal(a$Omega_u, 1.1 * second_diff_mat(17)) + expect_equal(a$Omega_v, 1.2 * second_diff_mat(8)) + + a <- moma_twfpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = 1.1 * second_diff_mat(17), + alpha_v = seq(0, 2, 0.2), Omega_v = 1.2 * second_diff_mat(8), + selection_scheme_str = "gb" + ) + expect_true(all(a$selection_scheme_list == c(0, 1, 0, 0))) + + a <- moma_twfpca(X, + alpha_v = seq(0, 2, 0.2), Omega_v = 2.1 * second_diff_mat(8), + selection_scheme_str = "bg" + ) + expect_true(all(a$selection_scheme_list == c(1, 0, 0, 0))) +}) From a6e654c46f19c57d16bd4819faa5fc62589d7360 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Tue, 9 Jul 2019 23:41:38 +0800 Subject: [PATCH 06/32] Update test names --- tests/testthat/test_sfpca_wrapper.R | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index b54468f1..d4b35f62 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -151,7 +151,9 @@ test_that("SFPCA object: left-project fucntion", { -test_that("SFPCA object: moma_spca", { + + +test_that("SFPCA object wrappers: moma_spca", { X <- matrix(runif(17 * 8), 17, 8) # test inputs @@ -232,7 +234,7 @@ test_that("SFPCA object: moma_spca", { }) -test_that("SFPCA object: moma_twspca", { +test_that("SFPCA object wrappers: moma_twspca", { X <- matrix(runif(17 * 8), 17, 8) # test inputs @@ -330,7 +332,7 @@ test_that("SFPCA object: moma_twspca", { }) -test_that("SFPCA object: moma_fpca", { +test_that("SFPCA object wrappers: moma_fpca", { X <- matrix(runif(17 * 8), 17, 8) # test inputs @@ -425,7 +427,7 @@ test_that("SFPCA object: moma_fpca", { }) -test_that("SFPCA object: moma_twfpca", { +test_that("SFPCA object wrappers: moma_twfpca", { X <- matrix(runif(17 * 8), 17, 8) # test inputs From 8bda1c9fa155ce15797e2c0c9c612b10a3c804d1 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Wed, 10 Jul 2019 00:07:46 +0800 Subject: [PATCH 07/32] Add TODO --- R/moma_expose.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/moma_expose.R b/R/moma_expose.R index ff3ffd5c..995cc38e 100644 --- a/R/moma_expose.R +++ b/R/moma_expose.R @@ -37,7 +37,7 @@ check_omega <- function(Omega, alpha, n) { moma_error("Omage_u/v is not a matrix.") } - # store them as sparse matrices using the package Matrix + # TODO: store them as sparse matrices using the package Matrix if (length(alpha) == 1 && alpha == 0) { # discard the Omega matrix specified by users Omega <- diag(n) From d309ec2e10092d935bf48a0aac807dfd50d4a969 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Wed, 10 Jul 2019 00:09:54 +0800 Subject: [PATCH 08/32] Omage -> Omega --- R/moma_expose.R | 2 +- tests/testthat/test_sfpca_wrapper.R | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/moma_expose.R b/R/moma_expose.R index 995cc38e..42bbee9a 100644 --- a/R/moma_expose.R +++ b/R/moma_expose.R @@ -34,7 +34,7 @@ check_omega <- function(Omega, alpha, n) { # check if Omega is a matrix if (!is.matrix(Omega) && !is.null(Omega)) { - moma_error("Omage_u/v is not a matrix.") + moma_error("Omega_u/v is not a matrix.") } # TODO: store them as sparse matrices using the package Matrix diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index d4b35f62..304c8ea9 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -345,7 +345,7 @@ test_that("SFPCA object wrappers: moma_fpca", { alpha_u = seq(0, 2, 0.2), Omega_u = lasso(), selection_scheme_str = "g" ), - "Omage_u/v is not a matrix." + "Omega_u/v is not a matrix." ) # test when no penalty matrix is provided From 414c02ff53f25a69959b3b74bfef4a83a6ddc7c4 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Wed, 10 Jul 2019 14:25:11 +0800 Subject: [PATCH 09/32] Add a test (when empty list is passed in) --- tests/testthat/test_sfpca_wrapper.R | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 304c8ea9..33aef76e 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -201,6 +201,18 @@ test_that("SFPCA object wrappers: moma_spca", { ) + expect_error(a <- moma_spca(X, lambda_u = c()), + paste0( + "All penalty levels (", + sQuote("lambda_u"), ", ", + sQuote("lambda_v"), ", ", + sQuote("alpha_u"), ", ", + sQuote("alpha_v"), + ") must be numeric." + ), + fixed = TRUE + ) + a <- moma_spca(X) expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) From 751622844ca6bb94c73bf72c095635eefd7d7fb9 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Tue, 9 Jul 2019 15:48:42 +0800 Subject: [PATCH 10/32] set_*_gric -> construct_*_search --- src/moma_level1.cpp | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/moma_level1.cpp b/src/moma_level1.cpp index 0a12b92c..41d03e79 100644 --- a/src/moma_level1.cpp +++ b/src/moma_level1.cpp @@ -3,7 +3,7 @@ #include "moma.h" // auxiliary functions for MoMA::grid_BIC_mix -const arma::vec &set_greedy_grid(const arma::vec &grid, int want_grid) +const arma::vec &construct_grid_for_search(const arma::vec &grid, int want_grid) { if (want_grid == 1) { @@ -15,7 +15,7 @@ const arma::vec &set_greedy_grid(const arma::vec &grid, int want_grid) } } -arma::vec set_bic_grid(const arma::vec &grid, int want_bic, int i) +arma::vec construct_grid_no_search(const arma::vec &grid, int want_bic, int i) { if (want_bic == 1) { @@ -170,10 +170,10 @@ Rcpp::List MoMA::grid_BIC_mix(const arma::vec &alpha_u, // grid_au = alpha_u, bic_au_grid = [-1]. // If alpha_u is selected via nested BIC search, // then grid_au = [-1], bic_au_grid = alpha_u - const arma::vec &grid_lu = set_greedy_grid(lambda_u, !selection_criterion_lambda_u); - const arma::vec &grid_lv = set_greedy_grid(lambda_v, !selection_criterion_lambda_v); - const arma::vec &grid_au = set_greedy_grid(alpha_u, !selection_criterion_alpha_u); - const arma::vec &grid_av = set_greedy_grid(alpha_v, !selection_criterion_alpha_v); + const arma::vec &grid_lu = construct_grid_for_search(lambda_u, !selection_criterion_lambda_u); + const arma::vec &grid_lv = construct_grid_for_search(lambda_v, !selection_criterion_lambda_v); + const arma::vec &grid_au = construct_grid_for_search(alpha_u, !selection_criterion_alpha_u); + const arma::vec &grid_av = construct_grid_for_search(alpha_v, !selection_criterion_alpha_v); // Test that if a grid is set to be BIC-search grid, then // the above code should set grid_xx to the vector [-1] @@ -204,10 +204,14 @@ Rcpp::List MoMA::grid_BIC_mix(const arma::vec &alpha_u, { for (int m = 0; m < n_lambda_v; m++) { - arma::vec bic_au_grid = set_bic_grid(alpha_u, selection_criterion_alpha_u, i); - arma::vec bic_lu_grid = set_bic_grid(lambda_u, selection_criterion_lambda_u, j); - arma::vec bic_av_grid = set_bic_grid(alpha_v, selection_criterion_alpha_v, k); - arma::vec bic_lv_grid = set_bic_grid(lambda_v, selection_criterion_lambda_v, m); + arma::vec bic_au_grid = + construct_grid_no_search(alpha_u, selection_criterion_alpha_u, i); + arma::vec bic_lu_grid = + construct_grid_no_search(lambda_u, selection_criterion_lambda_u, j); + arma::vec bic_av_grid = + construct_grid_no_search(alpha_v, selection_criterion_alpha_v, k); + arma::vec bic_lv_grid = + construct_grid_no_search(lambda_v, selection_criterion_lambda_v, m); if ((selection_criterion_alpha_u == 0 && bic_au_grid.n_elem != 1) || (selection_criterion_alpha_v == 0 && bic_av_grid.n_elem != 1) || From 29acaf83a2cdeb73674610c4560509d2126a32c3 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Wed, 10 Jul 2019 23:45:15 +0800 Subject: [PATCH 11/32] Add documentations --- R/moma_arguments.R | 28 +++++++++++++++++++ R/moma_sfpca.R | 69 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 97 insertions(+) diff --git a/R/moma_arguments.R b/R/moma_arguments.R index 40981a3c..b94ca54a 100644 --- a/R/moma_arguments.R +++ b/R/moma_arguments.R @@ -1,8 +1,27 @@ +#' Sparsity-inducing penalty in \code{MoMA} +#' In the package \code{MoMA}, we support the following sparsity-inducing +#' penalty functions. +#' \itemize{ +#' \item{\code{\link{empty}}} No penalty. TODO +#' \item{\code{\link{lasso}}} TODO +#' \item{\code{\link{mcp}}} TODO +#' \item{\code{\link{scad}}} TODO +#' \item{\code{\link{slope}}} TODO +#' \item{\code{\link{grplasso}}} TODO +#' \item{\code{\link{fusedlasso}}} TODO +#' \item{\code{\link{l1tf}}} TODO +#' \item{\code{\link{spfusedlasso}}} TODO +#' \item{\code{\link{cluster}}} TODO +#' } +#' @name moma_sparsity +NULL + # Check whether `x` is a boolean value is_logical_scalar <- function(x) { return(is.logical(x) && (length(x) == 1) && !is.na(x)) } +#' @export empty <- function() { arglist <- list() class(arglist) <- "moma_sparsity" @@ -312,6 +331,15 @@ cluster <- function(..., w = NULL, ADMM = FALSE, return(arglist) } + +#' Algorithm settings for solving a penalzied SVD problem +#' +#' To find an (approximate) solution to a penalized SVD (Singular Value Decomposition) problem is to solve two +#' penalized regression problems iteratively. Each penalized regression +#' is solved using one of the three algorithms: ISTA (Iterative Shrinkage-Thresholding Algorithm), +#' FISTA (Fast Iterative Shrinkage-Thresholding Algorithm) and +#' One-step ISTA (an approximated version of ISTA). +#' @export moma_pg_settings <- function(..., EPS = 1e-10, MAX_ITER = 1000, EPS_inner = 1e-10, MAX_ITER_inner = 1e+5, solver = c("ista", "fista", "onestepista")) { diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index a5a04f8f..a307ffee 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -1,3 +1,16 @@ +#' SFPCA R6 object +#' +#' An R6 object for performing SFPCA and parameter selection +#' +#' During initialziation of an \code{SFPCA} object, \code{R} +#' calls the \code{C++}-side function, \code{cpp_multirank_BIC_grid_search}, and +#' wraps the results returned. The \code{SFPCA} object also records penalty levels +#' and selection schemes of tuning parameters. Several helper +#' methods are provivded to facilitate access to results. +#' @seealso \code{\link[MoMA:moma_sfpca]{moma_sfpca}}, +#' \code{\link[MoMA:moma_spca]{moma_spca}}, moma_fpca, moma_twspca, moma_twfpca +#' @return nothing; modifies \code{person} +#' @export SFPCA <- R6::R6Class("SFPCA", list( center = NULL, scale = NULL, @@ -260,6 +273,41 @@ SFPCA <- R6::R6Class("SFPCA", list( } )) +#' Perform two-way sparse and functional PCA +#' +#' \code{moma_sfpca} creates an \code{SFPCA} R6 object and returns. +#' @param center a logical value indicating whether the variables should be shifted to be zero centered. +#' Defaults to \code{TRUE}. +#' @param scale a logical value indicating whether the variables should be scaled to have unit variance. +#' Defaults to \code{FALSE}. +#' @param u_sparsity,v_sparsity an object of class inheriting from "\code{moma_sparsity}". Most conveniently +#' specified by functions described in \code{\link{moma_sparsity}}. It specifies the type of sparsity-inducing +#' penalty function used in the model. Note that for \code{moma_spca}, these two parameter must not be +#' specified at the same time. For \code{moma_fpca} and \code{moma_twfpca}, they must not be specified. +#' @param lambda_u,lambda_v a numeric vector or a number that specifies the penalty level of sparsity. +#' Note that for \code{moma_spca}, these two parameter must not be +#' specified at the same time. For \code{moma_fpca} and \code{moma_twfpca}, they must not be specified. +#' @param Omega_u,Omega_v a positive definite matrix that encourages smoothness. Note that for \code{moma_fpca}, these two parameter must not be +#' specified at the same time. For \code{moma_spca} and \code{moma_twspca}, they must not be specified. +#' @param alpha_u,alpha_v v a numeric vector or a number that specifies the penalty level of smoothness. +#' Note that for \code{moma_fpca}, these two parameter must not be +#' specified at the same time. For \code{moma_spca} and \code{moma_twspca}, they must not be specified. +#' @param pg_setting an object of class inheriting from "\code{moma_sparsity}". Most conviently +#' specified by functions described in \code{\link{moma_pg_settings}}. It specifies the type of algorithm +#' used to solve the problem, acceptable level of precision, and the maximum number of iterations allowed. +#' @param selection_scheme_str a one-letter, two-letter or four-letter string that specifies selection schemes for tuning +#' parameters, containing only "b" and "g". "b" stands for greedy nested BIC selection, and "g" +#' stands for exhaustive grid search. For \code{moma_sfpca}, it is a four-letter string, and selection schemes for the tuning parameters, +#' alpha_u, alpha_v, lambda_u and lambda_v are specified by the four letters of the string, respectively. For +#' \code{moma_spca} and \code{moma_fpca}, it is a one-leter string that specifies the selection scheme +#' for the paramter of interest. For \code{moma_twspca}, it is a two-letter string, and +#' selection schemes for lambda_u and lambda_v are specified by the two letters respectively. For +#' \code{moma_twfpca}, it is a two-letter string, and selection schemes for alpha_u and alpha_v +#' are specified by the two letters respectively. +#' @param max_bic_iter a positive integer. Defaults to 5. The maximum number of iterations allowed +#' in nested greedy BIC selection scheme. +#' @param rank a positive integer. Defaults to 1. The maximal rank, i.e., maximal number of principal components to be used. +#' @export moma_sfpca <- function(X, ..., center = TRUE, scale = FALSE, u_sparsity = empty(), v_sparsity = empty(), lambda_u = 0, lambda_v = 0, # lambda_u/_v is a vector or scalar @@ -285,6 +333,11 @@ moma_sfpca <- function(X, ..., )) } +#' Perform one-way sparse PCA +#' +#' \code{moma_spca} is a wrapper around R6 object \code{SFPCA} +#' @export +#' @describeIn moma_sfpca a function for one-way sparse PCA moma_spca <- function(X, ..., center = TRUE, scale = FALSE, u_sparsity = empty(), v_sparsity = empty(), lambda_u = 0, lambda_v = 0, # lambda_u/_v is a vector or scalar @@ -334,6 +387,12 @@ moma_spca <- function(X, ..., # moma_error("Not implemented: SPCA") } + +#' Perform two-way sparse PCA +#' +#' \code{moma_twspca} is a wrapper around R6 object \code{SFPCA} +#' @export +#' @describeIn moma_sfpca a function for two-way sparse PCA moma_twspca <- function(X, ..., center = TRUE, scale = FALSE, u_sparsity = empty(), v_sparsity = empty(), lambda_u = 0, lambda_v = 0, # lambda_u/_v is a vector or scalar @@ -376,6 +435,11 @@ moma_twspca <- function(X, ..., )) } +#' Perform one-way functional PCA +#' +#' \code{moma_fpca} is a wrapper around R6 object \code{SFPCA} +#' @export +#' @describeIn moma_sfpca a function for one-way functional PCA moma_fpca <- function(X, ..., center = TRUE, scale = FALSE, Omega_u = NULL, Omega_v = NULL, alpha_u = 0, alpha_v = 0, @@ -424,6 +488,11 @@ moma_fpca <- function(X, ..., )) } +#' Perform two-way functional PCA +#' +#' \code{moma_twfpca} is a wrapper around R6 object \code{SFPCA} +#' @export +#' @describeIn moma_sfpca a function for two-way functional PCA moma_twfpca <- function(X, ..., center = TRUE, scale = FALSE, Omega_u = NULL, Omega_v = NULL, alpha_u = 0, alpha_v = 0, From feb992cd6e11b251ae3abf1fc4e81fe2c5b6a90d Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Thu, 11 Jul 2019 00:28:43 +0800 Subject: [PATCH 12/32] Add documentations for R6 object SFPCA --- R/moma_sfpca.R | 41 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 36 insertions(+), 5 deletions(-) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index a307ffee..401dcf0f 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -1,15 +1,46 @@ -#' SFPCA R6 object +#' \code{SFPCA} \code{R6} object #' -#' An R6 object for performing SFPCA and parameter selection +#' An \code{R6} object for performing SFPCA and parameter selection #' #' During initialziation of an \code{SFPCA} object, \code{R} #' calls the \code{C++}-side function, \code{cpp_multirank_BIC_grid_search}, and #' wraps the results returned. The \code{SFPCA} object also records penalty levels #' and selection schemes of tuning parameters. Several helper #' methods are provivded to facilitate access to results. -#' @seealso \code{\link[MoMA:moma_sfpca]{moma_sfpca}}, -#' \code{\link[MoMA:moma_spca]{moma_spca}}, moma_fpca, moma_twspca, moma_twfpca -#' @return nothing; modifies \code{person} +#' Initialization is delegated to \code{\link{moma_sfpca}}. +#' @seealso \code{\link{moma_sfpca}}, +#' \code{\link{moma_spca}}, +#' \code{\link{moma_fpca}}, +#' \code{\link{moma_twspca}}, +#' \code{\link{moma_twfpca}} +#' +#' @section Members: +#' +#' \describe{ +#' \item{\code{center,scale}}{The attributes "\code{scaled:center}" and "\code{scaled:scale}" of function \code{scale}. +#' The numeric centering and scalings used (if any) of the data matrix.} +#' \item{\code{grid_result}}{a 5-D list containing the results evaluated on the paramter grid.} +#' \item{\code{selection_scheme_list}}{a list with elements \code{selection_criterion_alpha_u}, +#' \code{selection_criterion_alpha_v}, +#' \code{selection_criterion_lambda_u}, +#' \code{selection_criterion_lambda_v}. Each of them is either 0 or 1. 0 stands for grid search +#' and 1 stands for BIC search. TODO: descibe the mixed selection procedure.} +#' } +#' @section Methods: +#' +#' \describe{ +#' \item{\code{get_mat}}{Arguments: \code{alpha_u}, \code{alpha_v}, \code{lambda_u}, \code{lambda_v} +#' are the indices of the parameters in the paramter grid, which is specified during initialization. +#' +#' Obtain the right and left sigular penalized vectors, which are packed into matrices \code{U} and \code{V}.} +#' \item{\code{print}}{Display tuning parameters and selection schemes.} +#' \item{\code{left_project}}{Arguments: \code{newX}, a matrix (un-centered and un-scaled) of +#' the same number of columns as the original data matrix. +#' +#' Project the new data into the space spaned by the +#' penalized left singular vectors, after scaling and centering as needed.} +#' } +#' @return an \code{SFPCA} R6 object. #' @export SFPCA <- R6::R6Class("SFPCA", list( center = NULL, From f37529dd4492813c9ab80978b4bc9b8fc158cd79 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Thu, 11 Jul 2019 01:04:26 +0800 Subject: [PATCH 13/32] Add more documentations to fix the build --- R/moma_arguments.R | 10 +++++++--- R/moma_sfpca.R | 2 ++ 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/R/moma_arguments.R b/R/moma_arguments.R index b94ca54a..1e777a16 100644 --- a/R/moma_arguments.R +++ b/R/moma_arguments.R @@ -2,7 +2,6 @@ #' In the package \code{MoMA}, we support the following sparsity-inducing #' penalty functions. #' \itemize{ -#' \item{\code{\link{empty}}} No penalty. TODO #' \item{\code{\link{lasso}}} TODO #' \item{\code{\link{mcp}}} TODO #' \item{\code{\link{scad}}} TODO @@ -21,7 +20,6 @@ is_logical_scalar <- function(x) { return(is.logical(x) && (length(x) == 1) && !is.na(x)) } -#' @export empty <- function() { arglist <- list() class(arglist) <- "moma_sparsity" @@ -335,10 +333,16 @@ cluster <- function(..., w = NULL, ADMM = FALSE, #' Algorithm settings for solving a penalzied SVD problem #' #' To find an (approximate) solution to a penalized SVD (Singular Value Decomposition) problem is to solve two -#' penalized regression problems iteratively. Each penalized regression +#' penalized regression problems iteratively (outer loop). Each penalized regression (inner loop) #' is solved using one of the three algorithms: ISTA (Iterative Shrinkage-Thresholding Algorithm), #' FISTA (Fast Iterative Shrinkage-Thresholding Algorithm) and #' One-step ISTA (an approximated version of ISTA). +#' @param ... to force users to specify arguments by names +#' @param EPS precision for outer loop +#' @param MAX_ITER the maximum number of iterations for outer loop +#' @param EPS_inner precision for inner loop +#' @param MAX_ITER_inner the maximum number of iterations for inner loop +#' @param solver a string in \code{c("ista", "fista", "onestepista")}. #' @export moma_pg_settings <- function(..., EPS = 1e-10, MAX_ITER = 1000, EPS_inner = 1e-10, MAX_ITER_inner = 1e+5, diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index 401dcf0f..f34ed766 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -307,6 +307,8 @@ SFPCA <- R6::R6Class("SFPCA", list( #' Perform two-way sparse and functional PCA #' #' \code{moma_sfpca} creates an \code{SFPCA} R6 object and returns. +#' @param X data matrix. +#' @param ... force users to specify arguments by names #' @param center a logical value indicating whether the variables should be shifted to be zero centered. #' Defaults to \code{TRUE}. #' @param scale a logical value indicating whether the variables should be scaled to have unit variance. From a1ad1933ba0e1aa4cf9a9081f0b4e04da8747b65 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Thu, 11 Jul 2019 17:48:33 +0800 Subject: [PATCH 14/32] t(V) %%*% V -> crossprod(V) --- R/moma_sfpca.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index f34ed766..de117416 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -291,7 +291,7 @@ SFPCA <- R6::R6Class("SFPCA", list( moma_error("`newX` is incompatible with orignal data.") } - PV <- solve((t(V) %*% V), t(V)) + PV <- solve(crossprod(V), t(V)) scaled_data <- scale(newX, self$center, self$scale) result <- scaled_data %*% t(PV) colnames(result) <- paste0("PC", seq_len(self$rank)) From 50b07b3730fbb5562aa91dd73ba0e85f76333c49 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Thu, 11 Jul 2019 21:21:08 +0800 Subject: [PATCH 15/32] get_mat -> get_mat_by_id, and dots --- R/moma_sfpca.R | 11 +++++----- tests/testthat/test_sfpca_wrapper.R | 31 ++++++++++++++++++++++++----- 2 files changed, 32 insertions(+), 10 deletions(-) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index de117416..cc1fd2d8 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -29,7 +29,7 @@ #' @section Methods: #' #' \describe{ -#' \item{\code{get_mat}}{Arguments: \code{alpha_u}, \code{alpha_v}, \code{lambda_u}, \code{lambda_v} +#' \item{\code{get_mat_by_index}}{Arguments: \code{alpha_u}, \code{alpha_v}, \code{lambda_u}, \code{lambda_v} #' are the indices of the parameters in the paramter grid, which is specified during initialization. #' #' Obtain the right and left sigular penalized vectors, which are packed into matrices \code{U} and \code{V}.} @@ -199,12 +199,13 @@ SFPCA <- R6::R6Class("SFPCA", list( ) }, - get_mat = function(alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) { + get_mat_by_index = function(..., alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) { + chkDots(...) # Sanity check: if a parameter has been chosen by BIC, then # the index for that parameter should not be specified. parameter <- c(alpha_u, alpha_v, lambda_u, lambda_v) if (any(parameter != 1 & self$selection_scheme_list != 0)) { - moma_error("Invalid index in SFPCA::get_mat. Do not specify indexes of parameters chosen by BIC.") + moma_error("Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters chosen by BIC.") } @@ -271,9 +272,9 @@ SFPCA <- R6::R6Class("SFPCA", list( # check indexes parameter <- c(alpha_u, alpha_v, lambda_u, lambda_v) if (any(parameter != 1 & self$selection_scheme_list != 0)) { - moma_error("Invalid index in SFPCA::get_mat. Do not specify indexes of parameters chosen by BIC.") + moma_error("Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters chosen by BIC.") } - V <- self$get_mat( + V <- self$get_mat_by_index( alpha_u = alpha_u, alpha_v = alpha_v, lambda_u = lambda_u, diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 33aef76e..7794c9c4 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -55,7 +55,13 @@ test_that("SFPCA object: as SVD", { a <- SFPCA$new(X, rank = 3, center = FALSE, scale = FALSE) - mysvd <- a$get_mat(1, 1, 1, 1) + mysvd <- a$get_mat_by_index( + alpha_u = 1, + alpah_v = 1, + lambda_u = 1, + lambda_v = 1 + ) + svda <- svd(X, nu = 3, nv = 3) expect_equal(abs(mysvd$U), abs(svda$u), check.attributes = FALSE) # use `abs` to avoid sign difference @@ -68,14 +74,24 @@ test_that("SFPCA object: as SVD", { dimnames(X) <- list(rown, coln) a <- SFPCA$new(X, rank = 3, center = FALSE, scale = FALSE) - mysvd <- a$get_mat(1, 1, 1, 1) + mysvd <- a$get_mat_by_index( + alpha_u = 1, + alpah_v = 1, + lambda_u = 1, + lambda_v = 1 + ) expect_equal(rownames(mysvd$U), rown) expect_equal(rownames(mysvd$V), coln) # test default arguments - expect_equal(a$get_mat(), a$get_mat(1, 1, 1, 1)) + expect_equal(a$get_mat_by_index(), a$get_mat_by_index( + alpha_u = 1, + alpah_v = 1, + lambda_u = 1, + lambda_v = 1 + )) # test selection with BIC seach and grid search a <- SFPCA$new(X, @@ -86,8 +102,13 @@ test_that("SFPCA object: as SVD", { expect_equal(dim(a$grid_result), c(1, 1, 1, 1, 3)) expect_error( - a$get_mat(2, 1, 1, 1), - "Invalid index in SFPCA::get_mat. Do not specify indexes of parameters chosen by BIC." + a$get_mat_by_index( + alpha_u = 2, + alpah_v = 1, + lambda_u = 1, + lambda_v = 1 + ), + "Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters chosen by BIC." ) }) From 1c4b91718c6d80e000b01b8c5e600f0952baf094 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Thu, 11 Jul 2019 21:36:56 +0800 Subject: [PATCH 16/32] Add is.wholenumber(rank) and tests --- R/moma_sfpca.R | 4 +++- tests/testthat/test_sfpca_wrapper.R | 21 +++++++++++++++++++++ 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index cc1fd2d8..13ba8567 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -158,7 +158,9 @@ SFPCA <- R6::R6Class("SFPCA", list( self$selection_scheme_list <- selection_scheme_list # Step 1.7: check rank - if (!inherits(rank, "numeric") || rank <= 0 + is.wholenumber <- + function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol + if (!inherits(rank, "numeric") || !is.wholenumber(rank) || rank <= 0 || rank > min(p, n)) { moma_error("rank should be a legit positive integer.") } diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 7794c9c4..be40f3e4 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -49,6 +49,27 @@ test_that("SFPCA object: correct arguments", { test_that("SFPCA object: as SVD", { + X <- matrix(runif(12), 4, 3) + expect_error( + SFPCA$new(X, rank = "3.1", center = FALSE, scale = FALSE), + "rank should be a legit positive integer" + ) + expect_error( + SFPCA$new(X, rank = 0, center = FALSE, scale = FALSE), + "rank should be a legit positive integer" + ) + expect_error( + SFPCA$new(X, rank = -1.1, center = FALSE, scale = FALSE), + "rank should be a legit positive integer" + ) + expect_error( + SFPCA$new(X, rank = -1, center = FALSE, scale = FALSE), + "rank should be a legit positive integer" + ) + expect_error( + SFPCA$new(X, rank = 3.1, center = FALSE, scale = FALSE), + "rank should be a legit positive integer" + ) # test no penalty case X <- matrix(runif(12), 4, 3) From 041ef1bab915de007ffda293bd2642d9be104bd6 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Thu, 11 Jul 2019 21:40:50 +0800 Subject: [PATCH 17/32] strsplit -> substr --- R/moma_sfpca.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index 13ba8567..6de43011 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -153,7 +153,7 @@ SFPCA <- R6::R6Class("SFPCA", list( selection_criterion_lambda_v = 0 ) for (i in 1:4) { - selection_scheme_list[[i]] <- ifelse(strsplit(selection_scheme_str, split = "")[[1]][i] == "g", 0, 1) + selection_scheme_list[[i]] <- ifelse(substr(selection_scheme_str, i, i) == "g", 0, 1) } self$selection_scheme_list <- selection_scheme_list From c9fdb8df3f3be67ac7e410b44a52a2416cfe3529 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Fri, 12 Jul 2019 01:03:47 +0800 Subject: [PATCH 18/32] Add SFPCA::interpolate with exact --- R/moma_expose.R | 2 ++ R/moma_sfpca.R | 51 +++++++++++++++++++++++++++++ tests/testthat/test_sfpca_wrapper.R | 2 ++ 3 files changed, 55 insertions(+) diff --git a/R/moma_expose.R b/R/moma_expose.R index 42bbee9a..952cd02e 100644 --- a/R/moma_expose.R +++ b/R/moma_expose.R @@ -28,6 +28,8 @@ add_default_prox_args <- function(sparsity_type) { return(modifyList(MOMA_DEFAULT_PROX, sparsity_type)) } +is.wholenumber <- + function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol # This function checks the validity of Omega and alpha check_omega <- function(Omega, alpha, n) { diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index 6de43011..9a00112d 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -56,6 +56,7 @@ SFPCA <- R6::R6Class("SFPCA", list( lambda_u = NULL, lambda_v = NULL, selection_scheme_list = NULL, + pg_setting = NULL, n = NULL, p = NULL, X = NULL, @@ -128,6 +129,7 @@ SFPCA <- R6::R6Class("SFPCA", list( ". Try using, for example, `pg_setting = moma_pg_settings(MAX_ITER=1e+4)`." ) } + self$pg_setting <- pg_setting # Step 1.5: smoothness Omega_u <- check_omega(Omega_u, alpha_u, n) @@ -205,6 +207,21 @@ SFPCA <- R6::R6Class("SFPCA", list( chkDots(...) # Sanity check: if a parameter has been chosen by BIC, then # the index for that parameter should not be specified. + + if (!is.numeric(c(alpha_u, alpha_v, lambda_u, lambda_v)) || + any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1)) { + moma_error("Non-scalar input in SFPCA::get_mat_by_index.") + } + + if (!all( + is.wholenumber(alpha_u), + is.wholenumber(alpha_v), + is.wholenumber(lambda_u), + is.wholenumber(lambda_v) + )) { + moma_error("SFPCA::get_mat_by_index takes integer indices.") + } + parameter <- c(alpha_u, alpha_v, lambda_u, lambda_v) if (any(parameter != 1 & self$selection_scheme_list != 0)) { moma_error("Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters chosen by BIC.") @@ -244,6 +261,40 @@ SFPCA <- R6::R6Class("SFPCA", list( return(list(U = U, V = V)) }, + interpolate = function(..., alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1, exact = FALSE) { + chkDots(...) + + if (!is.numeric(c(alpha_u, alpha_v, lambda_u, lambda_v)) || + any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1)) { + moma_error("Non-scalar input in SFPCA::interpolate.") + } + + if (is.unsorted(alpha_u) || + is.unsorted(alpha_v) || + is.unsorted(lambda_u) || + is.unsorted(lambda_v)) { + moma_error("Penalty levels not sorted!") + } + + if (exact) { + # call moma_svd + a <- moma_svd( + X = self$X, + u_sparsity = self$u_sparsity, v_sparsity = self$v_sparsity, + lambda_u = lambda_u, lambda_v = lambda_v, + Omega_u = self$Omega_u, Omega_v = self$Omega_v, + pg_setting = self$pg_setting, + k = self$rank, + ) + return(list(U = a$u, V = a$v)) + } + else { + # interpolate + a <- 1 + moma_error("interpolate option not implemented.") + } + }, + print = function() { selection_list_str <- lapply(self$selection_scheme_list, function(x) { if (x == 0) { diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index be40f3e4..87769eeb 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -191,7 +191,9 @@ test_that("SFPCA object: left-project fucntion", { ) }) +test_that("SFPCA object: interpolate", { +}) From c79b37592c0862f4d44df3b0fece7b6a1d045430 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Fri, 12 Jul 2019 19:08:37 +0800 Subject: [PATCH 19/32] Add `SFPCA::fixed_list` and check indices --- R/moma_expose.R | 3 +- R/moma_sfpca.R | 661 +++++++++++++++------------- tests/testthat/test_sfpca_wrapper.R | 149 +++++-- 3 files changed, 489 insertions(+), 324 deletions(-) diff --git a/R/moma_expose.R b/R/moma_expose.R index 952cd02e..d293d18e 100644 --- a/R/moma_expose.R +++ b/R/moma_expose.R @@ -47,7 +47,7 @@ check_omega <- function(Omega, alpha, n) { else if (is.null(Omega)) { # The user wants smooth penalty # but does not specify Omega matrix - Omega <- second_diff_mat(n) + Omega <- second_diff_mat(n) # TODO: should not overwrite } else { # At this point, users have specified an Omega and @@ -67,6 +67,7 @@ check_omega <- function(Omega, alpha, n) { ", but is actually ", dim(Omega)[1], "x", dim(Omega)[1] ) } + # TODO: check definiteness and symmetry of Omega } return(Omega) } diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index 9a00112d..d8d75b38 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -42,321 +42,392 @@ #' } #' @return an \code{SFPCA} R6 object. #' @export -SFPCA <- R6::R6Class("SFPCA", list( - center = NULL, - scale = NULL, - grid_result = NULL, - Omega_u = NULL, - Omega_v = NULL, - u_sparsity = NULL, - v_sparsity = NULL, - rank = NULL, - alpha_u = NULL, - alpha_v = NULL, - lambda_u = NULL, - lambda_v = NULL, - selection_scheme_list = NULL, - pg_setting = NULL, - n = NULL, - p = NULL, - X = NULL, - initialize = function(X, ..., - center = TRUE, scale = FALSE, - u_sparsity = empty(), v_sparsity = empty(), lambda_u = 0, lambda_v = 0, # lambda_u/_v is a vector or scalar - Omega_u = NULL, Omega_v = NULL, alpha_u = 0, alpha_v = 0, # so is alpha_u/_v - pg_setting = moma_pg_settings(), - selection_scheme_str = "gggg", - max_bic_iter = 5, - rank = 1) { - - # Step 1: check ALL arguments - # Step 1.1: lambdas and alphas - if (!inherits(alpha_u, c("numeric", "integer")) || - !inherits(alpha_v, c("numeric", "integer")) || - !inherits(lambda_u, c("numeric", "integer")) || - !inherits(lambda_v, c("numeric", "integer"))) { - moma_error(paste0( - "All penalty levels (", - sQuote("lambda_u"), ", ", - sQuote("lambda_v"), ", ", - sQuote("alpha_u"), ", ", - sQuote("alpha_v"), - ") must be numeric." - )) - } - self$alpha_u <- alpha_u - self$alpha_v <- alpha_v - self$lambda_u <- lambda_u - self$lambda_v <- lambda_v - - # Step 1.2: matrix - X <- as.matrix(X) - if (any(!is.finite(X))) { - moma_error("X must not have NaN, NA, or Inf.") - } - n <- dim(X)[1] - p <- dim(X)[2] - X <- scale(X, center = center, scale = scale) - - cen <- attr(X, "scaled:center") - sc <- attr(X, "scaled:scale") - if (any(sc == 0)) { - moma_error("cannot rescale a constant/zero column to unit variance") - } - - self$center <- if (is.null(cen)) FALSE else cen - self$scale <- if (is.null(sc)) FALSE else sc - self$n <- n - self$p <- p - self$X <- X - - # Step 1.3: sparsity - if (!inherits(u_sparsity, "moma_sparsity") || !inherits(v_sparsity, "moma_sparsity")) { - moma_error( - "Sparse penalty should be of class ", - sQuote("moma_sparsity"), - ". Try using, for example, `u_sparsity = lasso()`." +SFPCA <- R6::R6Class("SFPCA", + private = list( + check_input_index = TRUE, + private_get_mat_by_index = function(alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) { + private$check_input_index <- FALSE + res <- self$get_mat_by_index( + alpha_u, alpha_v, lambda_u, lambda_v ) + private$check_input_index <- TRUE + return(res) } - self$u_sparsity <- u_sparsity - self$v_sparsity <- v_sparsity - - # Step 1.4: PG loop settings - if (!inherits(pg_setting, "moma_pg_settings")) { - moma_error( - "pg_setting penalty should be of class ", - sQuote("moma_pg_settings"), - ". Try using, for example, `pg_setting = moma_pg_settings(MAX_ITER=1e+4)`." + ), + public = list( + center = NULL, + scale = NULL, + grid_result = NULL, + Omega_u = NULL, + Omega_v = NULL, + u_sparsity = NULL, + v_sparsity = NULL, + rank = NULL, + alpha_u = NULL, + alpha_v = NULL, + lambda_u = NULL, + lambda_v = NULL, + selection_scheme_list = NULL, + pg_setting = NULL, + n = NULL, + p = NULL, + X = NULL, + fixed_list = NULL, + initialize = function(X, ..., + center = TRUE, scale = FALSE, + u_sparsity = empty(), v_sparsity = empty(), lambda_u = 0, lambda_v = 0, # lambda_u/_v is a vector or scalar + Omega_u = NULL, Omega_v = NULL, alpha_u = 0, alpha_v = 0, # so is alpha_u/_v + pg_setting = moma_pg_settings(), + selection_scheme_str = "gggg", + max_bic_iter = 5, + rank = 1) { + + # Step 1: check ALL arguments + # Step 1.1: lambdas and alphas + if (!inherits(alpha_u, c("numeric", "integer")) || + !inherits(alpha_v, c("numeric", "integer")) || + !inherits(lambda_u, c("numeric", "integer")) || + !inherits(lambda_v, c("numeric", "integer"))) { + moma_error(paste0( + "All penalty levels (", + sQuote("lambda_u"), ", ", + sQuote("lambda_v"), ", ", + sQuote("alpha_u"), ", ", + sQuote("alpha_v"), + ") must be numeric." + )) + } + self$alpha_u <- alpha_u + self$alpha_v <- alpha_v + self$lambda_u <- lambda_u + self$lambda_v <- lambda_v + + # Step 1.2: matrix + X <- as.matrix(X) + if (any(!is.finite(X))) { + moma_error("X must not have NaN, NA, or Inf.") + } + n <- dim(X)[1] + p <- dim(X)[2] + X <- scale(X, center = center, scale = scale) + + cen <- attr(X, "scaled:center") + sc <- attr(X, "scaled:scale") + if (any(sc == 0)) { + moma_error("cannot rescale a constant/zero column to unit variance") + } + + self$center <- if (is.null(cen)) FALSE else cen + self$scale <- if (is.null(sc)) FALSE else sc + self$n <- n + self$p <- p + self$X <- X + + # Step 1.3: sparsity + if (!inherits(u_sparsity, "moma_sparsity") || !inherits(v_sparsity, "moma_sparsity")) { + moma_error( + "Sparse penalty should be of class ", + sQuote("moma_sparsity"), + ". Try using, for example, `u_sparsity = lasso()`." + ) + } + self$u_sparsity <- u_sparsity + self$v_sparsity <- v_sparsity + + # Step 1.4: PG loop settings + if (!inherits(pg_setting, "moma_pg_settings")) { + moma_error( + "pg_setting penalty should be of class ", + sQuote("moma_pg_settings"), + ". Try using, for example, `pg_setting = moma_pg_settings(MAX_ITER=1e+4)`." + ) + } + self$pg_setting <- pg_setting + + # Step 1.5: smoothness + Omega_u <- check_omega(Omega_u, alpha_u, n) + Omega_v <- check_omega(Omega_v, alpha_v, p) + self$Omega_u <- Omega_u + self$Omega_v <- Omega_v + + # Step 1.6: check selection scheme string + # "g" stands for grid search, "b" stands for BIC + if (!inherits(selection_scheme_str, "character") || nchar(selection_scheme_str) != 4 || + !all(strsplit(selection_scheme_str, split = "")[[1]] %in% c("b", "g"))) { + moma_error( + "Invalid selection_scheme_str ", selection_scheme_str, + ". It should be a four-char string containing only 'b' or 'g'." + ) + } + + # turn "b"/"g" to 1/0 + selection_scheme_list <- list( + selection_criterion_alpha_u = 0, + selection_criterion_alpha_v = 0, + selection_criterion_lambda_u = 0, + selection_criterion_lambda_v = 0 ) - } - self$pg_setting <- pg_setting - - # Step 1.5: smoothness - Omega_u <- check_omega(Omega_u, alpha_u, n) - Omega_v <- check_omega(Omega_v, alpha_v, p) - self$Omega_u <- Omega_u - self$Omega_v <- Omega_v - - # Step 1.6: check selection scheme string - # "g" stands for grid search, "b" stands for BIC - if (!inherits(selection_scheme_str, "character") || nchar(selection_scheme_str) != 4 || - !all(strsplit(selection_scheme_str, split = "")[[1]] %in% c("b", "g"))) { - moma_error( - "Invalid selection_scheme_str ", selection_scheme_str, - ". It should be a four-char string containing only 'b' or 'g'." + fixed_list <- list( + # three senarios when the followings can be true + # 1. it is selected by BIC + # 2. it is a fixed scalar + is_alpha_u_fixed = FALSE, + is_alpha_v_fixed = FALSE, + is_lambda_u_fixed = FALSE, + is_lambda_v_fixed = FALSE ) - } - # turn "b"/"g" to 1/0 - selection_scheme_list <- list( - selection_criterion_alpha_u = 0, - selection_criterion_alpha_v = 0, - selection_criterion_lambda_u = 0, - selection_criterion_lambda_v = 0 - ) - for (i in 1:4) { - selection_scheme_list[[i]] <- ifelse(substr(selection_scheme_str, i, i) == "g", 0, 1) - } - self$selection_scheme_list <- selection_scheme_list - - # Step 1.7: check rank - is.wholenumber <- - function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol - if (!inherits(rank, "numeric") || !is.wholenumber(rank) || rank <= 0 - || rank > min(p, n)) { - moma_error("rank should be a legit positive integer.") - } - self$rank <- rank + parameter_list <- list( + self$alpha_u, + self$alpha_v, + self$lambda_u, + self$lambda_v + ) + for (i in 1:4) { + selection_scheme_list[[i]] <- ifelse(substr(selection_scheme_str, i, i) == "g", 0, 1) + fixed_list[[i]] <- substr(selection_scheme_str, i, i) == "b" || length(parameter_list[[i]]) == 1 + } + self$selection_scheme_list <- selection_scheme_list + self$fixed_list <- fixed_list + + # Step 1.7: check rank + is.wholenumber <- + function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol + if (!inherits(rank, "numeric") || !is.wholenumber(rank) || rank <= 0 + || rank > min(p, n)) { + moma_error("rank should be a legit positive integer.") + } + self$rank <- rank + + # Step 2: pack all arguments in a list + algo_settings_list <- c( + list( + X = X, + lambda_u = lambda_u, + lambda_v = lambda_v, + # smoothness + alpha_u = alpha_u, + alpha_v = alpha_v, + rank = rank + ), + list( + Omega_u = check_omega(Omega_u, alpha_u, n), + Omega_v = check_omega(Omega_v, alpha_v, p), + prox_arg_list_u = add_default_prox_args(u_sparsity), + prox_arg_list_v = add_default_prox_args(v_sparsity) + ), + pg_setting, + selection_scheme_list, + list( + max_bic_iter = max_bic_iter + ) + ) + # make sure we explicitly specify ALL arguments + if (length(algo_settings_list) != length(formals(cpp_multirank_BIC_grid_search))) { + moma_error("Incomplete arguments in SFPCA::initialize.") + } - # Step 2: pack all arguments in a list - algo_settings_list <- c( - list( - X = X, - lambda_u = lambda_u, - lambda_v = lambda_v, - # smoothness - alpha_u = alpha_u, - alpha_v = alpha_v, - rank = rank - ), - list( - Omega_u = check_omega(Omega_u, alpha_u, n), - Omega_v = check_omega(Omega_v, alpha_v, p), - prox_arg_list_u = add_default_prox_args(u_sparsity), - prox_arg_list_v = add_default_prox_args(v_sparsity) - ), - pg_setting, - selection_scheme_list, - list( - max_bic_iter = max_bic_iter + # Step 3: call the fucntion + self$grid_result <- do.call( + cpp_multirank_BIC_grid_search, + algo_settings_list ) - ) - # make sure we explicitly specify ALL arguments - if (length(algo_settings_list) != length(formals(cpp_multirank_BIC_grid_search))) { - moma_error("Incomplete arguments in SFPCA::initialize.") - } + }, - # Step 3: call the fucntion - self$grid_result <- do.call( - cpp_multirank_BIC_grid_search, - algo_settings_list - ) - }, - - get_mat_by_index = function(..., alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) { - chkDots(...) - # Sanity check: if a parameter has been chosen by BIC, then - # the index for that parameter should not be specified. - - if (!is.numeric(c(alpha_u, alpha_v, lambda_u, lambda_v)) || - any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1)) { - moma_error("Non-scalar input in SFPCA::get_mat_by_index.") - } + get_mat_by_index = function(..., alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) { + if (private$check_input_index) { + chkDots(...) + } + # Sanity check: if a parameter has been chosen by BIC, then + # the index for that parameter should not be specified. - if (!all( - is.wholenumber(alpha_u), - is.wholenumber(alpha_v), - is.wholenumber(lambda_u), - is.wholenumber(lambda_v) - )) { - moma_error("SFPCA::get_mat_by_index takes integer indices.") - } + if (any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1 || + !all(is.wholenumber(alpha_u), is.wholenumber(alpha_v), is.wholenumber(lambda_u), is.wholenumber(lambda_v)))) { + moma_error("Non-integer input in SFPCA::get_mat_by_index.") + } - parameter <- c(alpha_u, alpha_v, lambda_u, lambda_v) - if (any(parameter != 1 & self$selection_scheme_list != 0)) { - moma_error("Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters chosen by BIC.") - } + if (!all( + is.wholenumber(alpha_u), + is.wholenumber(alpha_v), + is.wholenumber(lambda_u), + is.wholenumber(lambda_v) + )) { + moma_error("SFPCA::get_mat_by_index takes integer indices.") + } + if (private$check_input_index) { + missing_list <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) + if (any(self$fixed_list == TRUE & missing_list != TRUE)) { + moma_error( + paste0( + "Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters ", + "i) that are chosen by BIC, or ", + "ii) that are not specified during initialization of the SFCPA object, or ", + "iii) that are scalars during initialization of the SFCPA object." + ) + ) + } + } - n <- self$n - p <- self$p - rank <- self$rank - - U <- matrix(0, nrow = n, ncol = rank) - V <- matrix(0, nrow = p, ncol = rank) - for (i in (1:self$rank)) { - U[, i] <- - get_5Dlist_elem(self$grid_result, - alpha_u_i = alpha_u, - lambda_u_i = lambda_u, - alpha_v_i = alpha_v, - lambda_v_i = lambda_v, rank_i = i - )[[1]]$u$vector - - V[, i] <- - get_5Dlist_elem(self$grid_result, - alpha_u_i = alpha_u, - lambda_u_i = lambda_u, - alpha_v_i = alpha_v, - lambda_v_i = lambda_v, rank_i = i - )[[1]]$v$vector - } + n <- self$n + p <- self$p + rank <- self$rank + + U <- matrix(0, nrow = n, ncol = rank) + V <- matrix(0, nrow = p, ncol = rank) + for (i in (1:self$rank)) { + U[, i] <- + get_5Dlist_elem(self$grid_result, + alpha_u_i = alpha_u, + lambda_u_i = lambda_u, + alpha_v_i = alpha_v, + lambda_v_i = lambda_v, rank_i = i + )[[1]]$u$vector + + V[, i] <- + get_5Dlist_elem(self$grid_result, + alpha_u_i = alpha_u, + lambda_u_i = lambda_u, + alpha_v_i = alpha_v, + lambda_v_i = lambda_v, rank_i = i + )[[1]]$v$vector + } - coln <- if (is.null(colnames(self$X))) paste0("Xcol_", seq_len(p)) else colnames(self$X) - rown <- if (is.null(rownames(self$X))) paste0("Xrow_", seq_len(n)) else rownames(self$X) - dimnames(V) <- - list(coln, paste0("PC", seq_len(rank))) - dimnames(U) <- - list(rown, paste0("PC", seq_len(rank))) - return(list(U = U, V = V)) - }, - - interpolate = function(..., alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1, exact = FALSE) { - chkDots(...) - - if (!is.numeric(c(alpha_u, alpha_v, lambda_u, lambda_v)) || - any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1)) { - moma_error("Non-scalar input in SFPCA::interpolate.") - } + coln <- if (is.null(colnames(self$X))) paste0("Xcol_", seq_len(p)) else colnames(self$X) + rown <- if (is.null(rownames(self$X))) paste0("Xrow_", seq_len(n)) else rownames(self$X) + dimnames(V) <- + list(coln, paste0("PC", seq_len(rank))) + dimnames(U) <- + list(rown, paste0("PC", seq_len(rank))) + return(list(U = U, V = V)) + }, - if (is.unsorted(alpha_u) || - is.unsorted(alpha_v) || - is.unsorted(lambda_u) || - is.unsorted(lambda_v)) { - moma_error("Penalty levels not sorted!") - } + interpolate = function(..., alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1, exact = FALSE) { + chkDots(...) - if (exact) { - # call moma_svd - a <- moma_svd( - X = self$X, - u_sparsity = self$u_sparsity, v_sparsity = self$v_sparsity, - lambda_u = lambda_u, lambda_v = lambda_v, - Omega_u = self$Omega_u, Omega_v = self$Omega_v, - pg_setting = self$pg_setting, - k = self$rank, - ) - return(list(U = a$u, V = a$v)) - } - else { - # interpolate - a <- 1 - moma_error("interpolate option not implemented.") - } - }, - print = function() { - selection_list_str <- lapply(self$selection_scheme_list, function(x) { - if (x == 0) { - return("grid search") + if (!is.numeric(c(alpha_u, alpha_v, lambda_u, lambda_v)) || + any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1)) { + moma_error("Non-scalar input in SFPCA::interpolate.") } - else if (x == 1) { - return("BIC search") + + if (is.unsorted(alpha_u) || + is.unsorted(alpha_v) || + is.unsorted(lambda_u) || + is.unsorted(lambda_v)) { + moma_error("Penalty levels not sorted!") } - }) - - cat("An object containing solutions to the following settings\n") - cat("rank =", self$rank, "\n") - cat("Penalty and selection:\n") - - cat(paste0("alpha_u: ", selection_list_str[1]), "\n") - print(self$alpha_u) - cat(paste0("alpha_u: ", selection_list_str[2]), "\n") - print(self$alpha_v) - cat(paste0("lambda_u: ", selection_list_str[3]), "\n") - print(self$lambda_u) - cat(paste0("lambda_v: ", selection_list_str[4]), "\n") - print(self$lambda_v) - }, - - left_project = function(newX, ..., - alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) { - - # check indexes - parameter <- c(alpha_u, alpha_v, lambda_u, lambda_v) - if (any(parameter != 1 & self$selection_scheme_list != 0)) { - moma_error("Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters chosen by BIC.") - } - V <- self$get_mat_by_index( - alpha_u = alpha_u, - alpha_v = alpha_v, - lambda_u = lambda_u, - lambda_v = lambda_v - )$V - - - # newX should be uncencter and unscaled. - # check new X has same colnames - if (length(dim(newX)) != 2L) { - moma_error("'newX' must be a matrix or data frame") - } - if (dim(newX)[2] != self$p) { - moma_error("`newX` is incompatible with orignal data.") - } + missing_list <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) + if (any(self$fixed_list == TRUE & missing_list != TRUE)) { + moma_error( + paste0( + "Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters ", + "i) that are chosen by BIC, or ", + "ii) that are not specified during initialization of the SFCPA object, or ", + "iii) that are scalars during initialization of the SFCPA object." + ) + ) + } - PV <- solve(crossprod(V), t(V)) - scaled_data <- scale(newX, self$center, self$scale) - result <- scaled_data %*% t(PV) - colnames(result) <- paste0("PC", seq_len(self$rank)) - return(list( - scaled_data = scaled_data, - proj_data = result, - V = V, - PV = PV - )) - } -)) + # 1 replace the defualt with internal values + # 2 for those that has ranges but is not specified, pick the first one + if (exact) { + # call moma_svd + a <- moma_svd( + X = self$X, + u_sparsity = self$u_sparsity, v_sparsity = self$v_sparsity, + lambda_u = lambda_u, lambda_v = lambda_v, + Omega_u = self$Omega_u, Omega_v = self$Omega_v, + alpha_u = self$alpha_u, alpha_v = self$alpha, + pg_setting = self$pg_setting, + k = self$rank, + ) + return(list(U = a$u, V = a$v)) + } + else { + # interpolate + a <- 1 + moma_error("interpolate option not implemented.") + } + }, + + print = function() { + selection_list_str <- lapply(self$selection_scheme_list, function(x) { + if (x == 0) { + return("grid search") + } + else if (x == 1) { + return("BIC search") + } + }) + + cat("An object containing solutions to the following settings\n") + cat("rank =", self$rank, "\n") + cat("Penalty and selection:\n") + + cat(paste0("alpha_u: ", selection_list_str[1]), "\n") + print(self$alpha_u) + cat(paste0("alpha_u: ", selection_list_str[2]), "\n") + print(self$alpha_v) + cat(paste0("lambda_u: ", selection_list_str[3]), "\n") + print(self$lambda_u) + cat(paste0("lambda_v: ", selection_list_str[4]), "\n") + print(self$lambda_v) + }, + + left_project = function(newX, ..., + alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) { + + # check indexes + missing_list <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) + if (any(self$fixed_list == TRUE & missing_list != TRUE)) { + moma_error( + paste0( + "Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters ", + "i) that are chosen by BIC, or ", + "ii) that are not specified during initialization of the SFCPA object, or ", + "iii) that are scalars during initialization of the SFCPA object." + ) + ) + } + + if (any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1 || + !all(is.wholenumber(alpha_u), is.wholenumber(alpha_v), is.wholenumber(lambda_u), is.wholenumber(lambda_v)))) { + moma_error("Non-integer input in SFPCA::left_project.") + } + + V <- private$private_get_mat_by_index( + alpha_u = alpha_u, + alpha_v = alpha_v, + lambda_u = lambda_u, + lambda_v = lambda_v + )$V + + + # newX should be uncencter and unscaled. + # check new X has same colnames + if (length(dim(newX)) != 2L) { + moma_error("'newX' must be a matrix or data frame") + } + + if (dim(newX)[2] != self$p) { + moma_error("`newX` is incompatible with orignal data.") + } + + PV <- solve(crossprod(V), t(V)) + scaled_data <- scale(newX, self$center, self$scale) + result <- scaled_data %*% t(PV) + colnames(result) <- paste0("PC", seq_len(self$rank)) + return(list( + scaled_data = scaled_data, + proj_data = result, + V = V, + PV = PV + )) + } + ) +) #' Perform two-way sparse and functional PCA #' diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 87769eeb..154079cb 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -76,12 +76,7 @@ test_that("SFPCA object: as SVD", { a <- SFPCA$new(X, rank = 3, center = FALSE, scale = FALSE) - mysvd <- a$get_mat_by_index( - alpha_u = 1, - alpah_v = 1, - lambda_u = 1, - lambda_v = 1 - ) + mysvd <- a$get_mat_by_index() svda <- svd(X, nu = 3, nv = 3) @@ -95,24 +90,14 @@ test_that("SFPCA object: as SVD", { dimnames(X) <- list(rown, coln) a <- SFPCA$new(X, rank = 3, center = FALSE, scale = FALSE) - mysvd <- a$get_mat_by_index( - alpha_u = 1, - alpah_v = 1, - lambda_u = 1, - lambda_v = 1 - ) + mysvd <- a$get_mat_by_index() expect_equal(rownames(mysvd$U), rown) expect_equal(rownames(mysvd$V), coln) # test default arguments - expect_equal(a$get_mat_by_index(), a$get_mat_by_index( - alpha_u = 1, - alpah_v = 1, - lambda_u = 1, - lambda_v = 1 - )) + expect_equal(a$get_mat_by_index(), a$get_mat_by_index()) # test selection with BIC seach and grid search a <- SFPCA$new(X, @@ -122,15 +107,6 @@ test_that("SFPCA object: as SVD", { expect_true(all(a$selection_scheme_list == c(1, 0, 0, 0))) expect_equal(dim(a$grid_result), c(1, 1, 1, 1, 3)) - expect_error( - a$get_mat_by_index( - alpha_u = 2, - alpah_v = 1, - lambda_u = 1, - lambda_v = 1 - ), - "Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters chosen by BIC." - ) }) test_that("SFPCA object: print fucntion", { @@ -173,6 +149,7 @@ test_that("SFPCA object: left-project fucntion", { lambda_v = c(6), selection_scheme_str = "bbbb" ) + SFPCA$debug("left_project") expect_error( a$left_project(matrix(0, 4, 1)), "`newX` is incompatible with orignal data." @@ -191,8 +168,98 @@ test_that("SFPCA object: left-project fucntion", { ) }) -test_that("SFPCA object: interpolate", { +test_that("SFPCA object: `fixed_list` functions as expected", { + set.seed(113) + X <- matrix(runif(17 * 8), 17, 8) * 10 + + invaid_indices_error <- + paste0( + "Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters ", + "i) that are chosen by BIC, or ", + "ii) that are not specified during initialization of the SFCPA object, or ", + "iii) that are scalars during initialization of the SFCPA object." + ) + + # case 1: + # parameters that did not appear in initialization + # should not appear in SFPCA::get_mat_by_index at all + a <- moma_sfpca(X) + expect_no_error( + a$get_mat_by_index() + ) + expect_error( + a$get_mat_by_index(alpha_u = 1), + invaid_indices_error + ) + expect_warning( + a$get_mat_by_index(alphau = 1), + paste0("extra argument ", sQuote("alphau"), " will be disregarded") + ) + + # case 2: + # parameters that were scalars in initialization + # should not appear in SFPCA::get_mat_by_index either + a <- moma_sfpca(X, alpha_u = 1) + expect_no_error( + a$get_mat_by_index() + ) + expect_error( + a$get_mat_by_index(alpha_u = 1), + invaid_indices_error + ) + expect_error( + a$get_mat_by_index(alpha_v = 2), + invaid_indices_error + ) + expect_error( + a$get_mat_by_index(lambda_u = 1), + invaid_indices_error + ) + + + # case 2: + # parameters that were selected by BIC + # should not appear in SFPCA::get_mat_by_index either + a <- moma_sfpca(X, + alpha_u = c(1, 2), + alpha_v = c(1, 2, 3), # selected by BIC + selection_scheme_str = "gbgg" + ) + expect_no_error( + a$get_mat_by_index() + ) + expect_no_error( + a$get_mat_by_index(alpha_u = 2) + ) + expect_error( + a$get_mat_by_index(alpha_v = 0), + invaid_indices_error + ) + expect_error( + a$get_mat_by_index(lambda_u = 0), + invaid_indices_error + ) + + a <- moma_sfpca(X, alpha_u = c(1, 2)) + expect_no_error( + a$get_mat_by_index() + ) + expect_no_error( + a$get_mat_by_index(alpha_u = 2) + ) + expect_error( + a$get_mat_by_index(alpha_u = 2, alpha_v = 0), + invaid_indices_error + ) + expect_error( + a$get_mat_by_index(alpha_v = 0), + invaid_indices_error + ) + expect_error( + a$get_mat_by_index(lambda_u = 0), + invaid_indices_error + ) }) @@ -593,3 +660,29 @@ test_that("SFPCA object wrappers: moma_twfpca", { ) expect_true(all(a$selection_scheme_list == c(1, 0, 0, 0))) }) + +test_that("SFPCA object: get_mat_by_index and left_project takes non-ingeters", { + set.seed(113) + X <- matrix(runif(17 * 8), 17, 8) * 10 + + a <- moma_sfpca(X, + alpha_u = c(1, 2), alpha_v = c(1.1, 2.1), + lambda_u = c(1.3, 2.3), lambda_v = c(1.4, 2.4) + ) + + expect_no_error( + a$get_mat_by_index(alpha_u = 2) + ) + expect_error( + a$get_mat_by_index(alpha_u = 2.1), + "Non-integer input in SFPCA::get_mat_by_index" + ) + + expect_no_error( + a$left_project(X) + ) + expect_error( + a$left_project(X, alpha_u = 2.1), + "Non-integer input in SFPCA::left_project" + ) +}) From 1ca99e8598efac25fe10bf9ce2973b60c98dfcb7 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Fri, 12 Jul 2019 20:59:09 +0800 Subject: [PATCH 20/32] Add exact mode and tests --- R/moma_sfpca.R | 56 +++++++++++++++++------------ tests/testthat/test_sfpca_wrapper.R | 47 ++++++++++++++++++++++++ 2 files changed, 80 insertions(+), 23 deletions(-) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index d8d75b38..4212ee90 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -302,23 +302,47 @@ SFPCA <- R6::R6Class("SFPCA", return(list(U = U, V = V)) }, - interpolate = function(..., alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1, exact = FALSE) { + interpolate = function(..., alpha_u = 0, alpha_v = 0, lambda_u = 0, lambda_v = 0, exact = FALSE) { chkDots(...) + if (any(self$selection_scheme_list != 0)) { + moma_error("R6 ojbect SFPCA do not support interpolation when BIC selection scheme has been used.") + } + + missing_list <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) + + if (exact) { + # TODO: for moma_spca, moma_twspca etc., this can be relaxed + if (any(missing_list == TRUE)) { + moma_error("SFPCA::interpolate does not support exact mode unless all parameters are specified.") + } + + a <- moma_svd( + X = self$X, + # smoothness + u_sparsity = self$u_sparsity, v_sparsity = self$v_sparsity, + lambda_u = lambda_u, lambda_v = lambda_v, + # sparsity + Omega_u = self$Omega_u, Omega_v = self$Omega_v, + alpha_u = alpha_u, alpha_v = alpha_v, + pg_setting = self$pg_setting, + k = self$rank, + ) + return(list(U = a$u, V = a$v)) + } + if (!is.numeric(c(alpha_u, alpha_v, lambda_u, lambda_v)) || any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1)) { moma_error("Non-scalar input in SFPCA::interpolate.") } - if (is.unsorted(alpha_u) || - is.unsorted(alpha_v) || - is.unsorted(lambda_u) || - is.unsorted(lambda_v)) { - moma_error("Penalty levels not sorted!") - } - missing_list <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) + # If users specify a "fixed" parameter, give errors. + # "Fixed" parameters are those + # i) that are chosen by BIC, or ", + # ii) that are not specified during initialization of the SFCPA object, or " + # iii) that are scalars during initialization of the SFCPA object." if (any(self$fixed_list == TRUE & missing_list != TRUE)) { moma_error( paste0( @@ -330,21 +354,7 @@ SFPCA <- R6::R6Class("SFPCA", ) } - # 1 replace the defualt with internal values - # 2 for those that has ranges but is not specified, pick the first one - if (exact) { - # call moma_svd - a <- moma_svd( - X = self$X, - u_sparsity = self$u_sparsity, v_sparsity = self$v_sparsity, - lambda_u = lambda_u, lambda_v = lambda_v, - Omega_u = self$Omega_u, Omega_v = self$Omega_v, - alpha_u = self$alpha_u, alpha_v = self$alpha, - pg_setting = self$pg_setting, - k = self$rank, - ) - return(list(U = a$u, V = a$v)) - } + else { # interpolate a <- 1 diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 154079cb..3f41ed0e 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -686,3 +686,50 @@ test_that("SFPCA object: get_mat_by_index and left_project takes non-ingeters", "Non-integer input in SFPCA::left_project" ) }) + +test_that("SFPCA object: interpolate", { + set.seed(113) + X <- matrix(runif(17 * 8), 17, 8) * 10 + + a <- moma_sfpca(X, + u_sparsity = lasso(), v_sparsity = lasso(), + lambda_u = c(1, 2), lambda_v = c(1, 2), + alpha_u = c(1, 2), alpha_v = c(1, 2), + Omega_u = second_diff_mat(17), Omega_v = second_diff_mat(8), + selection_scheme_str = "gggb" + ) + expect_error( + a$interpolate(), + "R6 ojbect SFPCA do not support interpolation when BIC selection scheme has been used." + ) + + a <- moma_sfpca(X, + u_sparsity = lasso(), v_sparsity = lasso(), + lambda_u = 1.2, lambda_v = 1, + alpha_u = 0.7, alpha_v = 0.3, + Omega_u = second_diff_mat(17), Omega_v = second_diff_mat(8) + ) + expect_error( + a$interpolate( + alpha_u = 1, exact = TRUE + ), + "SFPCA::interpolate does not support exact mode unless all parameters are specified." + ) + + SFPCA$debug("interpolate") + internal_call <- a$interpolate( + alpha_u = 1, alpha_v = 1, + lambda_v = 1, lambda_u = 1, + exact = TRUE + ) + + dir_call <- moma_svd( + scale(X, center = a$center, scale = a$scale), + u_sparsity = lasso(), v_sparsity = lasso(), + alpha_u = 1, alpha_v = 1, + lambda_v = 1, lambda_u = 1, + Omega_u = second_diff_mat(17), Omega_v = second_diff_mat(8) + ) + expect_equal(internal_call$U, dir_call$u) + expect_equal(internal_call$V, dir_call$v) +}) From 0cb125b17fdf8d078c725b98db927c43823fdefc Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Sat, 13 Jul 2019 00:56:33 +0800 Subject: [PATCH 21/32] Add interpolate (inexact) and some tests --- R/moma_sfpca.R | 137 ++++++++++++++++++++++++---- tests/testthat/test_sfpca_wrapper.R | 61 ++++++++++++- 2 files changed, 178 insertions(+), 20 deletions(-) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index 4212ee90..869d703b 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -168,9 +168,10 @@ SFPCA <- R6::R6Class("SFPCA", selection_criterion_lambda_v = 0 ) fixed_list <- list( - # three senarios when the followings can be true - # 1. it is selected by BIC - # 2. it is a fixed scalar + # "Fixed" parameters are those + # i) that are chosen by BIC, or ", + # ii) that are not specified during initialization of the SFCPA object, or " + # iii) that are scalars during initialization of the SFCPA object." is_alpha_u_fixed = FALSE, is_alpha_v_fixed = FALSE, is_lambda_u_fixed = FALSE, @@ -310,6 +311,11 @@ SFPCA <- R6::R6Class("SFPCA", moma_error("R6 ojbect SFPCA do not support interpolation when BIC selection scheme has been used.") } + if (!is.numeric(c(alpha_u, alpha_v, lambda_u, lambda_v)) || + any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1)) { + moma_error("Non-scalar input in SFPCA::interpolate.") + } + missing_list <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) if (exact) { @@ -327,26 +333,15 @@ SFPCA <- R6::R6Class("SFPCA", Omega_u = self$Omega_u, Omega_v = self$Omega_v, alpha_u = alpha_u, alpha_v = alpha_v, pg_setting = self$pg_setting, - k = self$rank, + k = self$rank ) return(list(U = a$u, V = a$v)) } - if (!is.numeric(c(alpha_u, alpha_v, lambda_u, lambda_v)) || - any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1)) { - moma_error("Non-scalar input in SFPCA::interpolate.") - } - - - # If users specify a "fixed" parameter, give errors. - # "Fixed" parameters are those - # i) that are chosen by BIC, or ", - # ii) that are not specified during initialization of the SFCPA object, or " - # iii) that are scalars during initialization of the SFCPA object." if (any(self$fixed_list == TRUE & missing_list != TRUE)) { moma_error( paste0( - "Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters ", + "Invalid index in SFPCA::interpolate. Do not specify indexes of parameters ", "i) that are chosen by BIC, or ", "ii) that are not specified during initialization of the SFCPA object, or ", "iii) that are scalars during initialization of the SFCPA object." @@ -354,11 +349,115 @@ SFPCA <- R6::R6Class("SFPCA", ) } + if (is.unsorted(self$alpha_u) || + is.unsorted(self$alpha_v) || + is.unsorted(self$lambda_u) || + is.unsorted(self$lambda_v)) { + moma_error("Penalty levels not sorted!") + } + + # a bool: want to interpolate on u side? + inter_u <- !missing(alpha_u) && !missing(lambda_u) && + !self$fixed_list$is_alpha_u_fixed && + !self$fixed_list$is_lambda_u_fixed && + self$fixed_list$is_alpha_v_fixed && + self$fixed_list$is_lambda_v_fixed + + inter_v <- !missing(alpha_v) && !missing(lambda_v) && + self$fixed_list$is_alpha_u_fixed && + self$fixed_list$is_lambda_u_fixed && + !self$fixed_list$is_alpha_v_fixed && + !self$fixed_list$is_lambda_v_fixed + + if (!xor(inter_u, inter_v)) { + moma_error("SFPCA::interpolate only supports one-sided interpolation.") + } + + if (inter_v) { + + # test that it is in the known range + if (alpha_v >= max(self$alpha_v) || alpha_v <= min(self$alpha_v)) { + moma_error("Invalid range: alpha_v.") + } + # find the cloest alpha + alpha_v_i <- which.min(abs(self$alpha_v - alpha_v)) + + + # find the bin where lambda lies in + if (lambda_v >= max(self$lambda_v) || lambda_v <= min(self$lambda_v)) { + moma_error("Invalid range: lambda_v.") + } + lambda_v_i_lo <- findInterval(lambda_v, self$lambda_v) + lambda_v_i_hi <- lambda_v_i_lo + 1 + + + if (lambda_v_i_hi > length(self$lambda_v)) { + moma_error("SFPCA::interpolate, error in findInterval") + } + + result_lo <- private$private_get_mat_by_index( + alpha_u = 1, + alpha_v = alpha_v_i, + lambda_u = 1, + lambda_v = lambda_v_i_lo + ) + + result_hi <- private$private_get_mat_by_index( + alpha_u = 1, + alpha_v = alpha_v_i, + lambda_u = 1, + lambda_v = lambda_v_i_hi + ) + + U <- 0.5 * (result_lo$U + result_hi$U) + V <- 0.5 * (result_lo$V + result_hi$V) + + # newSv <- diag(self$p) + alpha_v * self$Omega_v + # newSu <- diag(self$n) + alpha_u * self$Omega_u + + # LvT <- chol(newSv) + # LuT <- chol(newSu) # L^T L = newS, L is a lower triagular matrix + + return(list(U = U, V = V)) + } + else if (inter_u) { + if (alpha_u >= max(self$alpha_u) || alpha_u <= min(self$alpha_u)) { + moma_error("Invalid range: alpha_u.") + } + # find the cloest alpha + alpha_u_i <- which.min(abs(self$alpha_u - alpha_u)) + + + # find the bin where lambda lies in + if (lambda_u >= max(self$lambda_u) || lambda_u <= min(self$lambda_u)) { + moma_error("Invalid range: lambda_u.") + } + lambda_u_i_lo <- findInterval(lambda_u, self$lambda_u) + lambda_u_i_hi <- lambda_u_i_lo + 1 + + if (lambda_u_i_hi > length(self$lambda_u)) { + moma_error("SFPCA::interpolate, error in findInterval.") + } + result_lo <- private$private_get_mat_by_index( + alpha_v = 1, + alpha_u = alpha_u_i, + lambda_v = 1, + lambda_u = lambda_u_i_lo + ) + + result_hi <- private$private_get_mat_by_index( + alpha_v = 1, + alpha_u = alpha_u_i, + lambda_v = 1, + lambda_u = lambda_u_i_hi + ) + + U <- 0.5 * (result_lo$U + result_hi$U) + V <- 0.5 * (result_lo$V + result_hi$V) + } else { - # interpolate - a <- 1 - moma_error("interpolate option not implemented.") + moma_error("UNKNOWN.") } }, diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 3f41ed0e..4a94af05 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -687,7 +687,7 @@ test_that("SFPCA object: get_mat_by_index and left_project takes non-ingeters", ) }) -test_that("SFPCA object: interpolate", { +test_that("SFPCA object: interpolate, exact mode", { set.seed(113) X <- matrix(runif(17 * 8), 17, 8) * 10 @@ -733,3 +733,62 @@ test_that("SFPCA object: interpolate", { expect_equal(internal_call$U, dir_call$u) expect_equal(internal_call$V, dir_call$v) }) + + +test_that("SFPCA object: interpolate, inexact mode", { + set.seed(113) + X <- matrix(runif(17 * 8), 17, 8) * 10 + + # interpolation on v side + a <- moma_sfpca(X, + u_sparsity = lasso(), v_sparsity = lasso(), + alpha_v = seq(0.1, 1, 0.15), alpha_u = 0.32, + Omega_v = second_diff_mat(8), Omega_u = second_diff_mat(17), + lambda_v = seq(0.1, 1, 0.2), lambda_u = 2.1 + ) + expect_no_error( + a$interpolate(alpha_v = 0.23, lambda_v = 0.121) + ) + + # error becasue both alpha_v and lambda_v should be specified + expect_error( + a$interpolate(alpha_v = 0.23), + "SFPCA::interpolate only supports one-sided interpolation." + ) + + # error becasue alpha_u is a scalar during initialization + expect_error( + a$interpolate(alpha_u = 0.2323), + paste0( + "Invalid index in SFPCA::interpolate. Do not specify indexes of parameters ", + "i) that are chosen by BIC, or ", + "ii) that are not specified during initialization of the SFCPA object, or ", + "iii) that are scalars during initialization of the SFCPA object." + ) + ) + + # error becasue both alpha_u should not be specified + expect_error( + a$interpolate(alpha_v = 0.23, lambda_v = 0.121, alpha_u = 1), + paste0( + "Invalid index in SFPCA::interpolate. Do not specify indexes of parameters ", + "i) that are chosen by BIC, or ", + "ii) that are not specified during initialization of the SFCPA object, or ", + "iii) that are scalars during initialization of the SFCPA object." + ) + ) + + # unsorted alpha_v + a <- moma_sfpca(X, + u_sparsity = lasso(), v_sparsity = lasso(), + alpha_v = c(seq(0.1, 1, 0.15), 0.02), + alpha_u = 0.32, + Omega_v = second_diff_mat(8), Omega_u = second_diff_mat(17), + lambda_v = seq(0.1, 1, 0.2), lambda_u = 2.1 + ) + + expect_error( + a$interpolate(alpha_v = 0.23, lambda_v = 0.121), + "Penalty levels not sorted" + ) +}) From 09dff04031723c62fc6739e263496b8b1883a900 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Sat, 13 Jul 2019 01:20:30 +0800 Subject: [PATCH 22/32] More comments --- R/moma_sfpca.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index 869d703b..a02c4ab2 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -242,11 +242,13 @@ SFPCA <- R6::R6Class("SFPCA", # Sanity check: if a parameter has been chosen by BIC, then # the index for that parameter should not be specified. + # they should be of length 1 if (any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1 || !all(is.wholenumber(alpha_u), is.wholenumber(alpha_v), is.wholenumber(lambda_u), is.wholenumber(lambda_v)))) { moma_error("Non-integer input in SFPCA::get_mat_by_index.") } + # indices should be integers if (!all( is.wholenumber(alpha_u), is.wholenumber(alpha_v), From 8242f10378bd2a67fceb0cf0ad78658477337bb0 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Sat, 13 Jul 2019 15:18:25 +0800 Subject: [PATCH 23/32] Imporved UX in SFPCA::interpolate --- R/moma_sfpca.R | 94 +++++++++++++-------- tests/testthat/test_sfpca_wrapper.R | 123 +++++++++++++++++++++------- 2 files changed, 152 insertions(+), 65 deletions(-) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index a02c4ab2..47d5640b 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -46,6 +46,8 @@ SFPCA <- R6::R6Class("SFPCA", private = list( check_input_index = TRUE, private_get_mat_by_index = function(alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) { + # private functions can be called only by + # internal functions private$check_input_index <- FALSE res <- self$get_mat_by_index( alpha_u, alpha_v, lambda_u, lambda_v @@ -169,9 +171,9 @@ SFPCA <- R6::R6Class("SFPCA", ) fixed_list <- list( # "Fixed" parameters are those - # i) that are chosen by BIC, or ", - # ii) that are not specified during initialization of the SFCPA object, or " - # iii) that are scalars during initialization of the SFCPA object." + # i) that are chosen by BIC, or + # ii) that are not specified during initialization of the SFCPA object, or + # iii) that are scalars as opposed to vectors during initialization of the SFCPA object. is_alpha_u_fixed = FALSE, is_alpha_v_fixed = FALSE, is_lambda_u_fixed = FALSE, @@ -239,8 +241,6 @@ SFPCA <- R6::R6Class("SFPCA", if (private$check_input_index) { chkDots(...) } - # Sanity check: if a parameter has been chosen by BIC, then - # the index for that parameter should not be specified. # they should be of length 1 if (any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1 || @@ -258,9 +258,16 @@ SFPCA <- R6::R6Class("SFPCA", moma_error("SFPCA::get_mat_by_index takes integer indices.") } + # A "fixed" parameter should not be specified + # at all (this is a bit stringent, can be improved later). + # "Fixed" parameters are those + # i) that are chosen by BIC, or + # ii) that are not specified during initialization of the SFCPA object, or + # iii) that are scalars as opposed to vectors during initialization of the SFCPA object. if (private$check_input_index) { - missing_list <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) - if (any(self$fixed_list == TRUE & missing_list != TRUE)) { + is_missing <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) + is_fixed <- self$fixed_list + if (any(is_fixed == TRUE & is_missing == FALSE)) { moma_error( paste0( "Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters ", @@ -308,23 +315,52 @@ SFPCA <- R6::R6Class("SFPCA", interpolate = function(..., alpha_u = 0, alpha_v = 0, lambda_u = 0, lambda_v = 0, exact = FALSE) { chkDots(...) - + # If BIC scheme has been used for any parameters, exit. if (any(self$selection_scheme_list != 0)) { moma_error("R6 ojbect SFPCA do not support interpolation when BIC selection scheme has been used.") } + # Reject inputs like alpha_u = "1" or "alpha_u" = c(1,2,3) if (!is.numeric(c(alpha_u, alpha_v, lambda_u, lambda_v)) || any(c(length(alpha_u), length(alpha_v), length(lambda_u), length(lambda_v)) > 1)) { moma_error("Non-scalar input in SFPCA::interpolate.") } - missing_list <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) + + # Parameters that are specified explictly is not "fixed". + # Parameters that are "fixed" must not be specified. + is_missing <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) + is_fixed <- self$fixed_list == TRUE + param_str_list <- c("alpha_u", "alpha_v", "lambda_u", "lambda_v") + if (any(is_missing == FALSE & is_fixed == TRUE)) { + output_para <- is_missing == FALSE & is_fixed == TRUE + moma_error( + paste0( + "Invalid index in SFPCA::interpolate: ", + paste(param_str_list[output_para], collapse = ","), + ". Do not specify indexes of parameters ", + "i) that are chosen by BIC, or ", + "ii) that are not specified during initialization of the SFCPA object, or ", + "iii) that are scalars during initialization of the SFCPA object." + ) + ) + } + if (any(is_missing == TRUE & is_fixed == FALSE)) { + output_para <- is_missing == TRUE & is_fixed == FALSE + moma_error( + paste0( + "Please spesify the following argument(s): ", + paste(param_str_list[output_para], collapse = ","), + "." + ) + ) + } if (exact) { - # TODO: for moma_spca, moma_twspca etc., this can be relaxed - if (any(missing_list == TRUE)) { - moma_error("SFPCA::interpolate does not support exact mode unless all parameters are specified.") - } + alpha_u <- ifelse(self$fixed_list$is_alpha_u_fixed, self$alpha_u, alpha_u) + alpha_v <- ifelse(self$fixed_list$is_alpha_v_fixed, self$alpha_v, alpha_v) + lambda_u <- ifelse(self$fixed_list$is_lambda_u_fixed, self$lambda_u, lambda_u) + lambda_v <- ifelse(self$fixed_list$is_lambda_v_fixed, self$lambda_v, lambda_v) a <- moma_svd( X = self$X, @@ -340,17 +376,7 @@ SFPCA <- R6::R6Class("SFPCA", return(list(U = a$u, V = a$v)) } - if (any(self$fixed_list == TRUE & missing_list != TRUE)) { - moma_error( - paste0( - "Invalid index in SFPCA::interpolate. Do not specify indexes of parameters ", - "i) that are chosen by BIC, or ", - "ii) that are not specified during initialization of the SFCPA object, or ", - "iii) that are scalars during initialization of the SFCPA object." - ) - ) - } - + # Function `findInterval` requires sorted breakpoints if (is.unsorted(self$alpha_u) || is.unsorted(self$alpha_v) || is.unsorted(self$lambda_u) || @@ -371,7 +397,8 @@ SFPCA <- R6::R6Class("SFPCA", !self$fixed_list$is_alpha_v_fixed && !self$fixed_list$is_lambda_v_fixed - if (!xor(inter_u, inter_v)) { + # If both of them are ture or false at the same time, exit. + if (inter_u == inter_v) { moma_error("SFPCA::interpolate only supports one-sided interpolation.") } @@ -381,19 +408,17 @@ SFPCA <- R6::R6Class("SFPCA", if (alpha_v >= max(self$alpha_v) || alpha_v <= min(self$alpha_v)) { moma_error("Invalid range: alpha_v.") } + # find the cloest alpha alpha_v_i <- which.min(abs(self$alpha_v - alpha_v)) - # find the bin where lambda lies in if (lambda_v >= max(self$lambda_v) || lambda_v <= min(self$lambda_v)) { moma_error("Invalid range: lambda_v.") } lambda_v_i_lo <- findInterval(lambda_v, self$lambda_v) lambda_v_i_hi <- lambda_v_i_lo + 1 - - - if (lambda_v_i_hi > length(self$lambda_v)) { + if (lambda_v_i_hi > length(self$lambda_v) || lambda_v_i_lo <= 0) { moma_error("SFPCA::interpolate, error in findInterval") } @@ -403,7 +428,6 @@ SFPCA <- R6::R6Class("SFPCA", lambda_u = 1, lambda_v = lambda_v_i_lo ) - result_hi <- private$private_get_mat_by_index( alpha_u = 1, alpha_v = alpha_v_i, @@ -436,8 +460,7 @@ SFPCA <- R6::R6Class("SFPCA", } lambda_u_i_lo <- findInterval(lambda_u, self$lambda_u) lambda_u_i_hi <- lambda_u_i_lo + 1 - - if (lambda_u_i_hi > length(self$lambda_u)) { + if (lambda_u_i_hi > length(self$lambda_u) || lambda_u_i_lo <= 0) { moma_error("SFPCA::interpolate, error in findInterval.") } @@ -447,7 +470,6 @@ SFPCA <- R6::R6Class("SFPCA", lambda_v = 1, lambda_u = lambda_u_i_lo ) - result_hi <- private$private_get_mat_by_index( alpha_v = 1, alpha_u = alpha_u_i, @@ -457,6 +479,7 @@ SFPCA <- R6::R6Class("SFPCA", U <- 0.5 * (result_lo$U + result_hi$U) V <- 0.5 * (result_lo$V + result_hi$V) + return(list(U = U, V = V)) } else { moma_error("UNKNOWN.") @@ -491,8 +514,9 @@ SFPCA <- R6::R6Class("SFPCA", alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) { # check indexes - missing_list <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) - if (any(self$fixed_list == TRUE & missing_list != TRUE)) { + is_missing <- list(missing(alpha_u), missing(alpha_v), missing(lambda_u), missing(lambda_v)) + is_fixed <- self$fixed_list + if (any(is_fixed == TRUE & is_missing == FALSE)) { moma_error( paste0( "Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters ", diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 4a94af05..762e2dfe 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -149,7 +149,6 @@ test_that("SFPCA object: left-project fucntion", { lambda_v = c(6), selection_scheme_str = "bbbb" ) - SFPCA$debug("left_project") expect_error( a$left_project(matrix(0, 4, 1)), "`newX` is incompatible with orignal data." @@ -703,6 +702,8 @@ test_that("SFPCA object: interpolate, exact mode", { "R6 ojbect SFPCA do not support interpolation when BIC selection scheme has been used." ) + # case 1: interpolate cannot be used if all + # parameters are specified as scalars a <- moma_sfpca(X, u_sparsity = lasso(), v_sparsity = lasso(), lambda_u = 1.2, lambda_v = 1, @@ -713,25 +714,87 @@ test_that("SFPCA object: interpolate, exact mode", { a$interpolate( alpha_u = 1, exact = TRUE ), - "SFPCA::interpolate does not support exact mode unless all parameters are specified." + "Invalid index in SFPCA::interpolate: alpha_u." + ) + expect_error( + a$interpolate( + alpha_v = 1, exact = TRUE + ), + "Invalid index in SFPCA::interpolate: alpha_v." + ) + expect_error( + a$interpolate( + lambda_u = 1, exact = TRUE + ), + "Invalid index in SFPCA::interpolate: lambda_u." + ) + expect_error( + a$interpolate( + lambda_v = 1, exact = TRUE + ), + "Invalid index in SFPCA::interpolate: lambda_v." ) - SFPCA$debug("interpolate") - internal_call <- a$interpolate( - alpha_u = 1, alpha_v = 1, - lambda_v = 1, lambda_u = 1, - exact = TRUE + + # case 2: interpolation on v side, both alpha and lambda are vectors + a <- moma_sfpca(X, + u_sparsity = lasso(), v_sparsity = lasso(), + alpha_v = seq(0.1, 1, 0.15), alpha_u = 0.32, + Omega_v = second_diff_mat(8), Omega_u = second_diff_mat(17), + lambda_v = seq(0.1, 1, 0.2), lambda_u = 2.1 + ) + # incomplete arguments + expect_error( + a$interpolate( + lambda_v = 1, exact = TRUE + ), + "Please spesify the following argument(s): alpha_v.", + fixed = TRUE + ) + # extra arguments + expect_error( + a$interpolate( + lambda_v = 1, alpha_v = 1, alpha_u = 1, + exact = TRUE + ), + "Invalid index in SFPCA::interpolate: alpha_u." + ) + # alpha_v too large + expect_error( + a$interpolate( + lambda_v = 0.21, alpha_v = 1 + ), + "Invalid range: alpha_v." ) - dir_call <- moma_svd( - scale(X, center = a$center, scale = a$scale), + # case 3: interpolation on v side, only lambda is a vector + a <- moma_sfpca(X, u_sparsity = lasso(), v_sparsity = lasso(), - alpha_u = 1, alpha_v = 1, - lambda_v = 1, lambda_u = 1, - Omega_u = second_diff_mat(17), Omega_v = second_diff_mat(8) + alpha_v = 0.12, alpha_u = 0.32, + Omega_v = second_diff_mat(8), Omega_u = second_diff_mat(17), + lambda_v = seq(0.1, 1, 0.2), lambda_u = 2.1 + ) + # correct arguments + expect_no_error( + a$interpolate( + lambda_v = 1, exact = TRUE + ) + ) + # extra arguments + expect_error( + a$interpolate( + lambda_v = 1, alpha_v = 1, alpha_u = 1, + exact = TRUE + ), + "Invalid index in SFPCA::interpolate: alpha_u." + ) + # alpha_v must not be specified + expect_error( + a$interpolate( + lambda_v = 0.21, alpha_v = 0.09 + ), + "Invalid index in SFPCA::interpolate: alpha_v." ) - expect_equal(internal_call$U, dir_call$u) - expect_equal(internal_call$V, dir_call$v) }) @@ -750,32 +813,32 @@ test_that("SFPCA object: interpolate, inexact mode", { a$interpolate(alpha_v = 0.23, lambda_v = 0.121) ) - # error becasue both alpha_v and lambda_v should be specified + # error because both alpha_v and lambda_v should be specified expect_error( a$interpolate(alpha_v = 0.23), - "SFPCA::interpolate only supports one-sided interpolation." + "Please spesify the following argument(s): lambda_v.", + fixed = TRUE + ) + expect_error( + a$interpolate(), + "lease spesify the following argument(s): alpha_v,lambda_v.", + fixed = TRUE ) - # error becasue alpha_u is a scalar during initialization + # error because alpha_u is a scalar during initialization expect_error( a$interpolate(alpha_u = 0.2323), - paste0( - "Invalid index in SFPCA::interpolate. Do not specify indexes of parameters ", - "i) that are chosen by BIC, or ", - "ii) that are not specified during initialization of the SFCPA object, or ", - "iii) that are scalars during initialization of the SFCPA object." - ) + "Invalid index in SFPCA::interpolate: alpha_u" ) - # error becasue both alpha_u should not be specified + # error because alpha_u should not be specified expect_error( a$interpolate(alpha_v = 0.23, lambda_v = 0.121, alpha_u = 1), - paste0( - "Invalid index in SFPCA::interpolate. Do not specify indexes of parameters ", - "i) that are chosen by BIC, or ", - "ii) that are not specified during initialization of the SFCPA object, or ", - "iii) that are scalars during initialization of the SFCPA object." - ) + "Invalid index in SFPCA::interpolate: alpha_u" + ) + expect_error( + a$interpolate(alpha_v = 0.23, lambda_v = 0.121, alpha_u = 1, lambda_u = 1.3), + "Invalid index in SFPCA::interpolate: alpha_u,lambda_u." ) # unsorted alpha_v From c696f0f3abb0f2b9436a9b49b6ccb3348bb9e5ee Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Sun, 14 Jul 2019 15:54:58 +0800 Subject: [PATCH 24/32] Fix typo and more comments --- R/moma_sfpca.R | 17 +++-- tests/testthat/test_sfpca_wrapper.R | 98 +++++++++++++++++++---------- 2 files changed, 77 insertions(+), 38 deletions(-) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index 47d5640b..53b56aca 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -198,7 +198,7 @@ SFPCA <- R6::R6Class("SFPCA", function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol if (!inherits(rank, "numeric") || !is.wholenumber(rank) || rank <= 0 || rank > min(p, n)) { - moma_error("rank should be a legit positive integer.") + moma_error("`rank` should be a positive integer smaller than the rank of the data matrix.") } self$rank <- rank @@ -309,7 +309,7 @@ SFPCA <- R6::R6Class("SFPCA", list(coln, paste0("PC", seq_len(rank))) dimnames(U) <- list(rown, paste0("PC", seq_len(rank))) - return(list(U = U, V = V)) + return(list(U = U, V = V, d = diag(t(U) %*% self$X %*% V))) }, interpolate = function(..., alpha_u = 0, alpha_v = 0, lambda_u = 0, lambda_v = 0, exact = FALSE) { @@ -317,7 +317,7 @@ SFPCA <- R6::R6Class("SFPCA", # If BIC scheme has been used for any parameters, exit. if (any(self$selection_scheme_list != 0)) { - moma_error("R6 ojbect SFPCA do not support interpolation when BIC selection scheme has been used.") + moma_error("R6 object SFPCA do not support interpolation when BIC selection scheme has been used.") } # Reject inputs like alpha_u = "1" or "alpha_u" = c(1,2,3) @@ -337,7 +337,7 @@ SFPCA <- R6::R6Class("SFPCA", moma_error( paste0( "Invalid index in SFPCA::interpolate: ", - paste(param_str_list[output_para], collapse = ","), + paste(param_str_list[output_para], collapse = ", "), ". Do not specify indexes of parameters ", "i) that are chosen by BIC, or ", "ii) that are not specified during initialization of the SFCPA object, or ", @@ -350,7 +350,7 @@ SFPCA <- R6::R6Class("SFPCA", moma_error( paste0( "Please spesify the following argument(s): ", - paste(param_str_list[output_para], collapse = ","), + paste(param_str_list[output_para], collapse = ", "), "." ) ) @@ -547,7 +547,12 @@ SFPCA <- R6::R6Class("SFPCA", } if (dim(newX)[2] != self$p) { - moma_error("`newX` is incompatible with orignal data.") + moma_error( + paste0( + "`newX` is incompatible with orignal data. ", + "It must have ", self$p, " columns." + ) + ) } PV <- solve(crossprod(V), t(V)) diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 762e2dfe..afe1d83e 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -52,23 +52,23 @@ test_that("SFPCA object: as SVD", { X <- matrix(runif(12), 4, 3) expect_error( SFPCA$new(X, rank = "3.1", center = FALSE, scale = FALSE), - "rank should be a legit positive integer" + "`rank` should be a positive integer smaller than the rank of the data matrix." ) expect_error( SFPCA$new(X, rank = 0, center = FALSE, scale = FALSE), - "rank should be a legit positive integer" + "`rank` should be a positive integer smaller than the rank of the data matrix." ) expect_error( SFPCA$new(X, rank = -1.1, center = FALSE, scale = FALSE), - "rank should be a legit positive integer" + "`rank` should be a positive integer smaller than the rank of the data matrix." ) expect_error( SFPCA$new(X, rank = -1, center = FALSE, scale = FALSE), - "rank should be a legit positive integer" + "`rank` should be a positive integer smaller than the rank of the data matrix." ) expect_error( SFPCA$new(X, rank = 3.1, center = FALSE, scale = FALSE), - "rank should be a legit positive integer" + "`rank` should be a positive integer smaller than the rank of the data matrix." ) # test no penalty case @@ -82,7 +82,7 @@ test_that("SFPCA object: as SVD", { expect_equal(abs(mysvd$U), abs(svda$u), check.attributes = FALSE) # use `abs` to avoid sign difference expect_equal(abs(mysvd$V), abs(svda$v), check.attributes = FALSE) - + expect_equal(mysvd$d, svda$d, check.attributes = FALSE) # test a matrix with dimnames rown <- paste0("row", seq(1:4)) @@ -96,8 +96,31 @@ test_that("SFPCA object: as SVD", { expect_equal(rownames(mysvd$V), coln) - # test default arguments - expect_equal(a$get_mat_by_index(), a$get_mat_by_index()) + # `a` is specified without + # any parameters (all four parameters default + # to 0 actually), so users must not specify any + # parameters when calling `get_mat_by_index`. + # A bit stringent though. + expect_error( + a$get_mat_by_index(alpha_u = 1), + "Invalid index in SFPCA::get_mat_by_index." + ) + expect_error( + a$get_mat_by_index(alpha_v = 1), + "Invalid index in SFPCA::get_mat_by_index." + ) + expect_error( + a$get_mat_by_index(lambda_u = 1), + "Invalid index in SFPCA::get_mat_by_index." + ) + expect_error( + a$get_mat_by_index(lambda_v = 1), + "Invalid index in SFPCA::get_mat_by_index." + ) + expect_no_error( + a$get_mat_by_index(), + "Invalid index in SFPCA::get_mat_by_index." + ) # test selection with BIC seach and grid search a <- SFPCA$new(X, @@ -121,7 +144,7 @@ test_that("SFPCA object: print fucntion", { ) print_message <- capture.output(print(a)) - expect_message <- c( + expected_message <- c( "An object containing solutions to the following settings", "rank = 3 ", "Penalty and selection:", @@ -134,12 +157,10 @@ test_that("SFPCA object: print fucntion", { "lambda_v: grid search ", "[1] 6 7" ) - expect_equal(expect_message, print_message) + expect_equal(expected_message, print_message) }) test_that("SFPCA object: left-project fucntion", { - - # unamed matrix X <- matrix(runif(17 * 8), 17, 8) a <- SFPCA$new(X, rank = 3, @@ -151,16 +172,16 @@ test_that("SFPCA object: left-project fucntion", { ) expect_error( a$left_project(matrix(0, 4, 1)), - "`newX` is incompatible with orignal data." + "`newX` is incompatible with orignal data. It must have 8 columns." ) new_data <- matrix(runif(24), 3, 8) res <- a$left_project(new_data) V <- res$V - # we verify the the projected data + # We verify the the projected data # satisfies the normal equation: - # X^T X b = X^T y + # X^T X b = X^T y. expect_equal( t(V) %*% V %*% t(res$proj_data), t(V) %*% t(res$scaled_data) @@ -171,7 +192,7 @@ test_that("SFPCA object: `fixed_list` functions as expected", { set.seed(113) X <- matrix(runif(17 * 8), 17, 8) * 10 - invaid_indices_error <- + invalid_indices_error <- paste0( "Invalid index in SFPCA::get_mat_by_index. Do not specify indexes of parameters ", "i) that are chosen by BIC, or ", @@ -181,15 +202,16 @@ test_that("SFPCA object: `fixed_list` functions as expected", { # case 1: # parameters that did not appear in initialization - # should not appear in SFPCA::get_mat_by_index at all + # should not appear in `SFPCA::get_mat_by_index` at all a <- moma_sfpca(X) expect_no_error( a$get_mat_by_index() ) expect_error( a$get_mat_by_index(alpha_u = 1), - invaid_indices_error + invalid_indices_error ) + # when an unused argument is given, or a typo. expect_warning( a$get_mat_by_index(alphau = 1), paste0("extra argument ", sQuote("alphau"), " will be disregarded") @@ -205,19 +227,19 @@ test_that("SFPCA object: `fixed_list` functions as expected", { ) expect_error( a$get_mat_by_index(alpha_u = 1), - invaid_indices_error + invalid_indices_error ) expect_error( a$get_mat_by_index(alpha_v = 2), - invaid_indices_error + invalid_indices_error ) expect_error( a$get_mat_by_index(lambda_u = 1), - invaid_indices_error + invalid_indices_error ) - # case 2: + # case 3: # parameters that were selected by BIC # should not appear in SFPCA::get_mat_by_index either a <- moma_sfpca(X, @@ -233,11 +255,11 @@ test_that("SFPCA object: `fixed_list` functions as expected", { ) expect_error( a$get_mat_by_index(alpha_v = 0), - invaid_indices_error + invalid_indices_error ) expect_error( a$get_mat_by_index(lambda_u = 0), - invaid_indices_error + invalid_indices_error ) a <- moma_sfpca(X, alpha_u = c(1, 2)) @@ -249,15 +271,15 @@ test_that("SFPCA object: `fixed_list` functions as expected", { ) expect_error( a$get_mat_by_index(alpha_u = 2, alpha_v = 0), - invaid_indices_error + invalid_indices_error ) expect_error( a$get_mat_by_index(alpha_v = 0), - invaid_indices_error + invalid_indices_error ) expect_error( a$get_mat_by_index(lambda_u = 0), - invaid_indices_error + invalid_indices_error ) }) @@ -699,10 +721,10 @@ test_that("SFPCA object: interpolate, exact mode", { ) expect_error( a$interpolate(), - "R6 ojbect SFPCA do not support interpolation when BIC selection scheme has been used." + "R6 object SFPCA do not support interpolation when BIC selection scheme has been used." ) - # case 1: interpolate cannot be used if all + # case 1: `SFPCA::interpolate` cannot be used at all if all # parameters are specified as scalars a <- moma_sfpca(X, u_sparsity = lasso(), v_sparsity = lasso(), @@ -710,6 +732,18 @@ test_that("SFPCA object: interpolate, exact mode", { alpha_u = 0.7, alpha_v = 0.3, Omega_u = second_diff_mat(17), Omega_v = second_diff_mat(8) ) + # this returns `a$get_mat_by_id()` + expect_no_error( + a$interpolate( + exact = TRUE + ) + ) + expect_error( + a$interpolate( + exact = FALSE + ), + "SFPCA::interpolate only supports one-sided interpolation" + ) expect_error( a$interpolate( alpha_u = 1, exact = TRUE @@ -786,7 +820,7 @@ test_that("SFPCA object: interpolate, exact mode", { lambda_v = 1, alpha_v = 1, alpha_u = 1, exact = TRUE ), - "Invalid index in SFPCA::interpolate: alpha_u." + "Invalid index in SFPCA::interpolate: alpha_u, alpha_v." ) # alpha_v must not be specified expect_error( @@ -821,7 +855,7 @@ test_that("SFPCA object: interpolate, inexact mode", { ) expect_error( a$interpolate(), - "lease spesify the following argument(s): alpha_v,lambda_v.", + "Please spesify the following argument(s): alpha_v, lambda_v.", fixed = TRUE ) @@ -838,7 +872,7 @@ test_that("SFPCA object: interpolate, inexact mode", { ) expect_error( a$interpolate(alpha_v = 0.23, lambda_v = 0.121, alpha_u = 1, lambda_u = 1.3), - "Invalid index in SFPCA::interpolate: alpha_u,lambda_u." + "Invalid index in SFPCA::interpolate: alpha_u, lambda_u." ) # unsorted alpha_v From bec0c766bd29c57acd8bfec044d3f80f3dd0918c Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Sun, 14 Jul 2019 15:59:18 +0800 Subject: [PATCH 25/32] Add `set.seed` --- tests/testthat/test_sfpca_wrapper.R | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index afe1d83e..1c324f9b 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -1,6 +1,7 @@ context("Test R6 object") -set.seed(12) test_that("SFPCA object: a naive case, 1x1 matrix", { + set.seed(12) + a <- SFPCA$new(matrix(1), center = FALSE) result_to_be_tested <- a$grid_result[[1]] @@ -21,6 +22,8 @@ test_that("SFPCA object: a naive case, 1x1 matrix", { }) test_that("SFPCA object: correct arguments", { + set.seed(12) + a <- SFPCA$new(matrix(runif(12), 3, 4), scale = FALSE, cen = FALSE) expect_equal(a$Omega_u, diag(3)) expect_equal(a$Omega_v, diag(4)) @@ -49,6 +52,8 @@ test_that("SFPCA object: correct arguments", { test_that("SFPCA object: as SVD", { + set.seed(12) + X <- matrix(runif(12), 4, 3) expect_error( SFPCA$new(X, rank = "3.1", center = FALSE, scale = FALSE), @@ -133,6 +138,8 @@ test_that("SFPCA object: as SVD", { }) test_that("SFPCA object: print fucntion", { + set.seed(12) + X <- matrix(runif(12), 4, 3) a <- SFPCA$new(X, rank = 3, @@ -161,6 +168,8 @@ test_that("SFPCA object: print fucntion", { }) test_that("SFPCA object: left-project fucntion", { + set.seed(12) + X <- matrix(runif(17 * 8), 17, 8) a <- SFPCA$new(X, rank = 3, @@ -286,6 +295,8 @@ test_that("SFPCA object: `fixed_list` functions as expected", { test_that("SFPCA object wrappers: moma_spca", { + set.seed(12) + X <- matrix(runif(17 * 8), 17, 8) # test inputs @@ -379,6 +390,8 @@ test_that("SFPCA object wrappers: moma_spca", { test_that("SFPCA object wrappers: moma_twspca", { + set.seed(12) + X <- matrix(runif(17 * 8), 17, 8) # test inputs @@ -477,6 +490,8 @@ test_that("SFPCA object wrappers: moma_twspca", { test_that("SFPCA object wrappers: moma_fpca", { + set.seed(12) + X <- matrix(runif(17 * 8), 17, 8) # test inputs @@ -572,6 +587,8 @@ test_that("SFPCA object wrappers: moma_fpca", { test_that("SFPCA object wrappers: moma_twfpca", { + set.seed(12) + X <- matrix(runif(17 * 8), 17, 8) # test inputs From 12d67f0498e6263df457e8c34e1afb1034f36850 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Sun, 14 Jul 2019 16:46:08 +0800 Subject: [PATCH 26/32] More comments for tests for interpolate --- tests/testthat/test_sfpca_wrapper.R | 31 +++++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 1c324f9b..600f6125 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -23,7 +23,7 @@ test_that("SFPCA object: a naive case, 1x1 matrix", { test_that("SFPCA object: correct arguments", { set.seed(12) - + a <- SFPCA$new(matrix(runif(12), 3, 4), scale = FALSE, cen = FALSE) expect_equal(a$Omega_u, diag(3)) expect_equal(a$Omega_v, diag(4)) @@ -729,10 +729,13 @@ test_that("SFPCA object: interpolate, exact mode", { set.seed(113) X <- matrix(runif(17 * 8), 17, 8) * 10 + # Once BIC is used, `SFPCA::interpolate` + # cannot be called at all. a <- moma_sfpca(X, u_sparsity = lasso(), v_sparsity = lasso(), lambda_u = c(1, 2), lambda_v = c(1, 2), - alpha_u = c(1, 2), alpha_v = c(1, 2), + alpha_u = c(1, 2), + alpha_v = c(1, 2), # selected by BIC Omega_u = second_diff_mat(17), Omega_v = second_diff_mat(8), selection_scheme_str = "gggb" ) @@ -740,6 +743,13 @@ test_that("SFPCA object: interpolate, exact mode", { a$interpolate(), "R6 object SFPCA do not support interpolation when BIC selection scheme has been used." ) + expect_error( + a$interpolate( + lambda_u = 1.1, lambda_v = 1.3, + alpha_u = 1.4 + ), + "R6 object SFPCA do not support interpolation when BIC selection scheme has been used." + ) # case 1: `SFPCA::interpolate` cannot be used at all if all # parameters are specified as scalars @@ -905,4 +915,21 @@ test_that("SFPCA object: interpolate, inexact mode", { a$interpolate(alpha_v = 0.23, lambda_v = 0.121), "Penalty levels not sorted" ) + + # interpolate on both sides is not supported + a <- moma_sfpca(X, + u_sparsity = lasso(), v_sparsity = lasso(), + alpha_v = seq(0.1, 1, 0.15), alpha_u = seq(0.1, 1, 0.15), + Omega_v = second_diff_mat(8), Omega_u = second_diff_mat(17), + lambda_v = seq(0.1, 1.4, 0.4), lambda_u = seq(0.1, 1.3, 0.45) + ) + # trying to do two-sided interpolation + # but not supported now. + expect_error( + a$interpolate( + alpha_v = 0.23, lambda_v = 0.121, + alpha_u = 0.1, lambda_u = 0.3 + ), + "SFPCA::interpolate only supports one-sided interpolation" + ) }) From ef9bb20acd71b5c43eed13d70d1a0fecf70037e7 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Mon, 15 Jul 2019 13:53:17 +0800 Subject: [PATCH 27/32] Add test "SFPCA object: X contains string" --- tests/testthat/test_sfpca_wrapper.R | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 600f6125..7006f781 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -21,6 +21,14 @@ test_that("SFPCA object: a naive case, 1x1 matrix", { expect_equal(result_to_be_tested$d, 1) }) +test_that("SFPCA object: X contains string", { + set.seed(12) + SFPCA$debug("initialize") + X <- matrix(seq(1:12), 4, 3) + X[1, 1] <- "abc" + expect_error(SFPCA$new(X), "X must not have NaN, NA, or Inf.") +}) + test_that("SFPCA object: correct arguments", { set.seed(12) From 948feaedb958d465dd2859d2566c4617b1a6cea2 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Wed, 17 Jul 2019 14:29:16 +0800 Subject: [PATCH 28/32] Replace is.null with %||% --- R/moma_sfpca.R | 8 ++++---- R/sfpca.R | 4 ++-- R/util.R | 7 +++++++ tests/testthat/test_sfpca_wrapper.R | 1 - 4 files changed, 13 insertions(+), 7 deletions(-) create mode 100644 R/util.R diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index 53b56aca..cde8031c 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -119,8 +119,8 @@ SFPCA <- R6::R6Class("SFPCA", moma_error("cannot rescale a constant/zero column to unit variance") } - self$center <- if (is.null(cen)) FALSE else cen - self$scale <- if (is.null(sc)) FALSE else sc + self$center <- cen %||% FALSE + self$scale <- sc %||% FALSE self$n <- n self$p <- p self$X <- X @@ -303,8 +303,8 @@ SFPCA <- R6::R6Class("SFPCA", )[[1]]$v$vector } - coln <- if (is.null(colnames(self$X))) paste0("Xcol_", seq_len(p)) else colnames(self$X) - rown <- if (is.null(rownames(self$X))) paste0("Xrow_", seq_len(n)) else rownames(self$X) + coln <- colnames(self$X) %||% paste0("Xcol_", seq_len(p)) + rown <- rownames(self$X) %||% paste0("Xrow_", seq_len(n)) dimnames(V) <- list(coln, paste0("PC", seq_len(rank))) dimnames(U) <- diff --git a/R/sfpca.R b/R/sfpca.R index 6e236e0a..dd34a100 100644 --- a/R/sfpca.R +++ b/R/sfpca.R @@ -55,8 +55,8 @@ sfpca <- function(X, lambda_u <- as.vector(lambda_u) lambda_v <- as.vector(lambda_v) - Omega_u <- if (is.null(Omega_u)) diag(dim(X)[1]) else Omega_u - Omega_v <- if (is.null(Omega_v)) diag(dim(X)[2]) else Omega_v + Omega_u <- Omega_u %||% diag(dim(X)[1]) + Omega_v <- Omega_v %||% diag(dim(X)[2]) prox_arg_list_u <- list( w = w_u, diff --git a/R/util.R b/R/util.R new file mode 100644 index 00000000..c1c07f71 --- /dev/null +++ b/R/util.R @@ -0,0 +1,7 @@ +`%||%` <- function(lhs, rhs) { + if (!is.null(lhs)) { + lhs + } else { + rhs + } +} diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 7006f781..8f736c49 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -23,7 +23,6 @@ test_that("SFPCA object: a naive case, 1x1 matrix", { test_that("SFPCA object: X contains string", { set.seed(12) - SFPCA$debug("initialize") X <- matrix(seq(1:12), 4, 3) X[1, 1] <- "abc" expect_error(SFPCA$new(X), "X must not have NaN, NA, or Inf.") From 70f35b04f05de0001582ef3f95eb269beaa8f75f Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Mon, 15 Jul 2019 19:13:33 +0800 Subject: [PATCH 29/32] Expect warning in `test_sfpca_wrapper.R ` --- tests/testthat/test_sfpca_wrapper.R | 158 ++++++++++++++++++---------- 1 file changed, 101 insertions(+), 57 deletions(-) diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index 8f736c49..ed9e77f0 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -177,7 +177,7 @@ test_that("SFPCA object: print fucntion", { test_that("SFPCA object: left-project fucntion", { set.seed(12) - X <- matrix(runif(17 * 8), 17, 8) + X <- matrix(runif(17 * 8), 17, 8) * 10 a <- SFPCA$new(X, rank = 3, alpha_u = c(1), @@ -304,7 +304,7 @@ test_that("SFPCA object: `fixed_list` functions as expected", { test_that("SFPCA object wrappers: moma_spca", { set.seed(12) - X <- matrix(runif(17 * 8), 17, 8) + X <- matrix(runif(17 * 8), 17, 8) * 10 # test inputs expect_error( @@ -363,10 +363,16 @@ test_that("SFPCA object wrappers: moma_spca", { fixed = TRUE ) - a <- moma_spca(X) + expect_warning( + a <- moma_spca(X), + "No sparsity is imposed!" + ) expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) - a <- moma_spca(X, selection_scheme_str = "b") # no effects + expect_warning( + a <- moma_spca(X, selection_scheme_str = "b"), # no effects + "No sparsity is imposed!" + ) expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) a <- moma_spca(X, @@ -399,7 +405,7 @@ test_that("SFPCA object wrappers: moma_spca", { test_that("SFPCA object wrappers: moma_twspca", { set.seed(12) - X <- matrix(runif(17 * 8), 17, 8) + X <- matrix(runif(17 * 8), 17, 8) * 10 # test inputs expect_no_error( @@ -430,30 +436,40 @@ test_that("SFPCA object wrappers: moma_twspca", { # test selection schemes - expect_error( - moma_twspca(X, - lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), - selection_scheme_str = "gbb" + expect_warning( + expect_error( + moma_twspca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + selection_scheme_str = "gbb" + ), + "`selection_scheme_str` should be of length two." ), - "`selection_scheme_str` should be of length two." + "Please use `moma_spca` if only one side is penalized." ) - expect_error( - moma_twspca(X, - lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), - selection_scheme_str = "cc" + expect_warning( + expect_error( + moma_twspca(X, + lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), + selection_scheme_str = "cc" + ), + "`selection_scheme_str` should consist of 'g' or 'b'" ), - "`selection_scheme_str` should consist of 'g' or 'b'" + "Please use `moma_spca` if only one side is penalized" ) - expect_no_error( - moma_twspca(X, - lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), - selection_scheme_str = "bb" - ) + expect_warning( + expect_no_error( + moma_twspca(X, + lambda_u = seq(0, 1, 0.1), u_sparsity = lasso(), + selection_scheme_str = "bb" + ) + ), + "Please use `moma_spca` if only one side is penalized" ) + expect_warning(a <- moma_twspca(X)) expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) @@ -475,9 +491,10 @@ test_that("SFPCA object wrappers: moma_twspca", { expect_true(all(a$selection_scheme_list == c(0, 0, 1, 1))) a <- moma_twspca(X, - lambda_u = seq(0, 2, 0.2), u_sparsity = lasso(), - lambda_v = seq(0, 2, 0.2), v_sparsity = lasso(), - selection_scheme_str = "gg" + lambda_u = seq(0, 1, 0.1), u_sparsity = lasso(), + lambda_v = seq(0, 1, 0.1), v_sparsity = lasso(), + selection_scheme_str = "gg", + pg_setting = moma_pg_settings(MAX_ITER = 2000) ) expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) @@ -488,9 +505,12 @@ test_that("SFPCA object wrappers: moma_twspca", { ) expect_true(all(a$selection_scheme_list == c(0, 0, 0, 1))) - a <- moma_twspca(X, - lambda_v = seq(0, 2, 0.2), v_sparsity = lasso(), - selection_scheme_str = "bg" + expect_warning( + a <- moma_twspca(X, + lambda_v = seq(0, 2, 0.2), v_sparsity = lasso(), + selection_scheme_str = "bg" + ), + "Please use `moma_spca` if only one side is penalized" ) expect_true(all(a$selection_scheme_list == c(0, 0, 1, 0))) }) @@ -499,7 +519,7 @@ test_that("SFPCA object wrappers: moma_twspca", { test_that("SFPCA object wrappers: moma_fpca", { set.seed(12) - X <- matrix(runif(17 * 8), 17, 8) + X <- matrix(runif(17 * 8), 17, 8) * 10 # test inputs expect_error( @@ -559,10 +579,16 @@ test_that("SFPCA object wrappers: moma_fpca", { ) - a <- moma_fpca(X) + expect_warning( + a <- moma_fpca(X), + "No smoothness is imposed!" + ) expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) - a <- moma_fpca(X, selection_scheme_str = "b") # no effects + expect_warning( + a <- moma_fpca(X, selection_scheme_str = "b"), # no effects + "No smoothness is imposed!" + ) expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) a <- moma_fpca(X, @@ -596,28 +622,36 @@ test_that("SFPCA object wrappers: moma_fpca", { test_that("SFPCA object wrappers: moma_twfpca", { set.seed(12) - X <- matrix(runif(17 * 8), 17, 8) + X <- matrix(runif(17 * 8), 17, 8) * 10 # test inputs expect_no_error( moma_twfpca(X, alpha_v = c(1, 2), alpha_u = c(2, 3)) ) - expect_no_error( - moma_twfpca(X, alpha_v = c(1, 2)) + + expect_warning( + expect_no_error( + moma_twfpca(X, alpha_v = c(1, 2)) + ), + "Please use `moma_fpca` if only one side is penalized" ) expect_no_error( moma_twfpca(X, Omega_v = diag(8), Omega_u = 2.1 * second_diff_mat(17)) ) # incompatible Omega - expect_error( - moma_twfpca(X, - alpha_u = seq(0, 2, 0.2), Omega_u = 2.1 * second_diff_mat(11), - selection_scheme_str = "bb" + expect_warning( + expect_error( + moma_twfpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = 2.1 * second_diff_mat(11), + selection_scheme_str = "bb" + ), + "Omega shoud be a compatible matrix. It should be of 17x17, but is actually 11x11" ), - "Omega shoud be a compatible matrix. It should be of 17x17, but is actually 11x11" + "Please use `moma_fpca` if only one side is penalized" ) + expect_no_error( moma_twfpca(X, Omega_v = diag(8), alpha_v = c(1, 2), @@ -636,31 +670,38 @@ test_that("SFPCA object wrappers: moma_twfpca", { # test selection schemes - expect_error( - moma_twfpca(X, - alpha_u = seq(0, 2, 0.2), Omega_u = lasso(), - selection_scheme_str = "gbb" + expect_warning( + expect_error( + moma_twfpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = lasso(), + selection_scheme_str = "gbb" + ), + "`selection_scheme_str` should be of length two." ), - "`selection_scheme_str` should be of length two." + "Please use `moma_fpca` if only one side is penalized" ) - expect_error( - moma_twfpca(X, - alpha_u = seq(0, 2, 0.2), Omega_u = lasso(), - selection_scheme_str = "cc" + expect_warning( + expect_error( + moma_twfpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = lasso(), + selection_scheme_str = "cc" + ), + "`selection_scheme_str` should consist of 'g' or 'b'" ), - "`selection_scheme_str` should consist of 'g' or 'b'" + "Please use `moma_fpca` if only one side is penalized" ) - expect_no_error( - moma_twfpca(X, - alpha_u = seq(0, 2, 0.2), Omega_u = 2.1 * second_diff_mat(17), - selection_scheme_str = "bb" - ) + expect_warning( + expect_no_error( + moma_twfpca(X, + alpha_u = seq(0, 2, 0.2), Omega_u = 2.1 * second_diff_mat(17), + selection_scheme_str = "bb" + ) + ), + "Please use `moma_fpca` if only one side is penalized" ) - - expect_warning(a <- moma_twfpca(X)) expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0))) @@ -699,9 +740,12 @@ test_that("SFPCA object wrappers: moma_twfpca", { ) expect_true(all(a$selection_scheme_list == c(0, 1, 0, 0))) - a <- moma_twfpca(X, - alpha_v = seq(0, 2, 0.2), Omega_v = 2.1 * second_diff_mat(8), - selection_scheme_str = "bg" + expect_warning( + a <- moma_twfpca(X, + alpha_v = seq(0, 2, 0.2), Omega_v = 2.1 * second_diff_mat(8), + selection_scheme_str = "bg" + ), + "Please use `moma_fpca` if only one side is penalized" ) expect_true(all(a$selection_scheme_list == c(1, 0, 0, 0))) }) From 0dedb2a404b8c5650be4855b46a16699ede56674 Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Wed, 17 Jul 2019 19:42:14 +0800 Subject: [PATCH 30/32] Add `d` in SFPCA$get_mat_by_id --- R/moma_sfpca.R | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index cde8031c..44ba7f57 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -285,22 +285,17 @@ SFPCA <- R6::R6Class("SFPCA", U <- matrix(0, nrow = n, ncol = rank) V <- matrix(0, nrow = p, ncol = rank) + d <- vector(mode = "numeric", length = rank) for (i in (1:self$rank)) { - U[, i] <- - get_5Dlist_elem(self$grid_result, - alpha_u_i = alpha_u, - lambda_u_i = lambda_u, - alpha_v_i = alpha_v, - lambda_v_i = lambda_v, rank_i = i - )[[1]]$u$vector - - V[, i] <- - get_5Dlist_elem(self$grid_result, - alpha_u_i = alpha_u, - lambda_u_i = lambda_u, - alpha_v_i = alpha_v, - lambda_v_i = lambda_v, rank_i = i - )[[1]]$v$vector + rank_i_result <- get_5Dlist_elem(self$grid_result, + alpha_u_i = alpha_u, + lambda_u_i = lambda_u, + alpha_v_i = alpha_v, + lambda_v_i = lambda_v, rank_i = i + )[[1]] + U[, i] <- rank_i_result$u$vector + V[, i] <- rank_i_result$v$vector + d[i] <- rank_i_result$d } coln <- colnames(self$X) %||% paste0("Xcol_", seq_len(p)) @@ -309,7 +304,7 @@ SFPCA <- R6::R6Class("SFPCA", list(coln, paste0("PC", seq_len(rank))) dimnames(U) <- list(rown, paste0("PC", seq_len(rank))) - return(list(U = U, V = V, d = diag(t(U) %*% self$X %*% V))) + return(list(U = U, V = V, d = d)) }, interpolate = function(..., alpha_u = 0, alpha_v = 0, lambda_u = 0, lambda_v = 0, exact = FALSE) { From 8ca700047e7951b1086c44a13c1d1b0e1d8a20fe Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Fri, 19 Jul 2019 14:07:25 +0800 Subject: [PATCH 31/32] Add a comment --- R/moma_sfpca.R | 1 + 1 file changed, 1 insertion(+) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index 44ba7f57..eef8bb8e 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -194,6 +194,7 @@ SFPCA <- R6::R6Class("SFPCA", self$fixed_list <- fixed_list # Step 1.7: check rank + # TODO: check that `rank` < min(rank(X), rank(Y)) is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol if (!inherits(rank, "numeric") || !is.wholenumber(rank) || rank <= 0 From 8cbbfa879c4fc5c74639af37aa3557257e64235d Mon Sep 17 00:00:00 2001 From: Luofeng Liao Date: Mon, 5 Aug 2019 15:23:23 +0800 Subject: [PATCH 32/32] Add `is.double(X)` --- R/moma_sfpca.R | 3 +++ tests/testthat/test_sfpca_wrapper.R | 8 ++++++-- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/R/moma_sfpca.R b/R/moma_sfpca.R index eef8bb8e..192c0cdc 100644 --- a/R/moma_sfpca.R +++ b/R/moma_sfpca.R @@ -106,6 +106,9 @@ SFPCA <- R6::R6Class("SFPCA", # Step 1.2: matrix X <- as.matrix(X) + if (!is.double(X)) { + moma_error("X must contain numbers only.") + } if (any(!is.finite(X))) { moma_error("X must not have NaN, NA, or Inf.") } diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R index ed9e77f0..44dfac7e 100644 --- a/tests/testthat/test_sfpca_wrapper.R +++ b/tests/testthat/test_sfpca_wrapper.R @@ -21,11 +21,15 @@ test_that("SFPCA object: a naive case, 1x1 matrix", { expect_equal(result_to_be_tested$d, 1) }) -test_that("SFPCA object: X contains string", { +test_that("SFPCA object: X contains strings", { set.seed(12) X <- matrix(seq(1:12), 4, 3) X[1, 1] <- "abc" - expect_error(SFPCA$new(X), "X must not have NaN, NA, or Inf.") + expect_error(SFPCA$new(X), "X must contain numbers only") + + X <- matrix(seq(1:12), 4, 3) + X[1, 1] <- Inf + expect_error(SFPCA$new(X), "X must not have NaN, NA, or Inf") }) test_that("SFPCA object: correct arguments", {