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_arguments.R b/R/moma_arguments.R
index 40981a3c..1e777a16 100644
--- a/R/moma_arguments.R
+++ b/R/moma_arguments.R
@@ -1,3 +1,20 @@
+#' Sparsity-inducing penalty in \code{MoMA}
+#' In the package \code{MoMA}, we support the following sparsity-inducing
+#' penalty functions.
+#' \itemize{
+#'   \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))
@@ -312,6 +329,21 @@ 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 (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,
                              solver = c("ista", "fista", "onestepista")) {
diff --git a/R/moma_expose.R b/R/moma_expose.R
index 7ec8f371..d293d18e 100644
--- a/R/moma_expose.R
+++ b/R/moma_expose.R
@@ -28,9 +28,18 @@ 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) {
+
+    # check if Omega is a matrix
+    if (!is.matrix(Omega) && !is.null(Omega)) {
+        moma_error("Omega_u/v is not a 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)
@@ -38,9 +47,12 @@ 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
+        # 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(
@@ -55,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
new file mode 100644
index 00000000..192c0cdc
--- /dev/null
+++ b/R/moma_sfpca.R
@@ -0,0 +1,833 @@
+#' \code{SFPCA} \code{R6} object
+#'
+#' 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.
+#' 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_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}.}
+#'   \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",
+    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
+            )
+            private$check_input_index <- TRUE
+            return(res)
+        }
+    ),
+    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 (!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.")
+            }
+            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 <- cen %||% FALSE
+            self$scale <- sc %||% FALSE
+            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
+            )
+            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 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,
+                is_lambda_v_fixed = FALSE
+            )
+
+            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
+            # 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
+            || rank > min(p, n)) {
+                moma_error("`rank` should be a positive integer smaller than the rank of the data matrix.")
+            }
+            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_by_index = function(..., alpha_u = 1, alpha_v = 1, lambda_u = 1, lambda_v = 1) {
+            if (private$check_input_index) {
+                chkDots(...)
+            }
+
+            # 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),
+                is.wholenumber(lambda_u),
+                is.wholenumber(lambda_v)
+            )) {
+                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) {
+                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 ",
+                            "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)
+            d <- vector(mode = "numeric", length = rank)
+            for (i in (1:self$rank)) {
+                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))
+            rown <- rownames(self$X) %||% paste0("Xrow_", seq_len(n))
+            dimnames(V) <-
+                list(coln, paste0("PC", seq_len(rank)))
+            dimnames(U) <-
+                list(rown, paste0("PC", seq_len(rank)))
+            return(list(U = U, V = V, d = d))
+        },
+
+        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 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)
+            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.")
+            }
+
+
+            # 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) {
+                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,
+                    # 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))
+            }
+
+            # Function `findInterval` requires sorted breakpoints
+            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 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.")
+            }
+
+            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) || lambda_v_i_lo <= 0) {
+                    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) || lambda_u_i_lo <= 0) {
+                    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)
+                return(list(U = U, V = V))
+            }
+            else {
+                moma_error("UNKNOWN.")
+            }
+        },
+
+        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 <SFPCA> 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
+            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 ",
+                        "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(
+                    paste0(
+                        "`newX` is incompatible with orignal data. ",
+                        "It must have ", self$p, " columns."
+                    )
+                )
+            }
+
+            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
+#'
+#' \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.
+#' 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
+                       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
+    ))
+}
+
+#' 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
+                      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")
+}
+
+
+#' 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
+                        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
+    ))
+}
+
+#' 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,
+                      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
+    ))
+}
+
+#' 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,
+                        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
+    ))
+}
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/src/moma_level1.cpp b/src/moma_level1.cpp
index b9876c18..41d03e79 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<Rcpp::NumericVector>(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<Rcpp::NumericVector>(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<Rcpp::NumericVector>(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<Rcpp::NumericVector>(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,
@@ -217,18 +237,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<Rcpp::NumericVector>(u_result["vector"]);
+                        arma::vec curv = Rcpp::as<Rcpp::NumericVector>(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<Rcpp::NumericVector>(u_result["vector"]);
-                            arma::vec curv = Rcpp::as<Rcpp::NumericVector>(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)
+})
diff --git a/tests/testthat/test_sfpca_wrapper.R b/tests/testthat/test_sfpca_wrapper.R
new file mode 100644
index 00000000..44dfac7e
--- /dev/null
+++ b/tests/testthat/test_sfpca_wrapper.R
@@ -0,0 +1,990 @@
+context("Test R6 object")
+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]]
+
+    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: X contains strings", {
+    set.seed(12)
+    X <- matrix(seq(1:12), 4, 3)
+    X[1, 1] <- "abc"
+    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", {
+    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))
+    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", {
+    set.seed(12)
+
+    X <- matrix(runif(12), 4, 3)
+    expect_error(
+        SFPCA$new(X, rank = "3.1", center = FALSE, scale = FALSE),
+        "`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 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 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 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 positive integer smaller than the rank of the data matrix."
+    )
+
+    # test no penalty case
+    X <- matrix(runif(12), 4, 3)
+    a <- SFPCA$new(X, rank = 3, center = FALSE, scale = FALSE)
+
+
+    mysvd <- a$get_mat_by_index()
+
+    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)
+    expect_equal(mysvd$d, svda$d, 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_by_index()
+
+    expect_equal(rownames(mysvd$U), rown)
+    expect_equal(rownames(mysvd$V), coln)
+
+
+    # `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,
+        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))
+})
+
+test_that("SFPCA object: print fucntion", {
+    set.seed(12)
+
+    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))
+
+    expected_message <- c(
+        "An <SFPCA> 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(expected_message, print_message)
+})
+
+test_that("SFPCA object: left-project fucntion", {
+    set.seed(12)
+
+    X <- matrix(runif(17 * 8), 17, 8) * 10
+    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. 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
+    # 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: `fixed_list` functions as expected", {
+    set.seed(113)
+    X <- matrix(runif(17 * 8), 17, 8) * 10
+
+    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 ",
+            "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),
+        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")
+    )
+
+
+    # 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),
+        invalid_indices_error
+    )
+    expect_error(
+        a$get_mat_by_index(alpha_v = 2),
+        invalid_indices_error
+    )
+    expect_error(
+        a$get_mat_by_index(lambda_u = 1),
+        invalid_indices_error
+    )
+
+
+    # case 3:
+    # 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),
+        invalid_indices_error
+    )
+    expect_error(
+        a$get_mat_by_index(lambda_u = 0),
+        invalid_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),
+        invalid_indices_error
+    )
+    expect_error(
+        a$get_mat_by_index(alpha_v = 0),
+        invalid_indices_error
+    )
+    expect_error(
+        a$get_mat_by_index(lambda_u = 0),
+        invalid_indices_error
+    )
+})
+
+
+
+test_that("SFPCA object wrappers: moma_spca", {
+    set.seed(12)
+
+    X <- matrix(runif(17 * 8), 17, 8) * 10
+
+    # 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'"
+    )
+
+
+    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
+    )
+
+    expect_warning(
+        a <- moma_spca(X),
+        "No sparsity is imposed!"
+    )
+    expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0)))
+
+    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,
+        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 wrappers: moma_twspca", {
+    set.seed(12)
+
+    X <- matrix(runif(17 * 8), 17, 8) * 10
+
+    # 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_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."
+        ),
+        "Please use `moma_spca` if only one side is penalized."
+    )
+
+    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'"
+        ),
+        "Please use `moma_spca` if only one side is penalized"
+    )
+
+    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)))
+
+    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, 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)))
+
+    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)))
+
+    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)))
+})
+
+
+test_that("SFPCA object wrappers: moma_fpca", {
+    set.seed(12)
+
+    X <- matrix(runif(17 * 8), 17, 8) * 10
+
+    # 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"
+        ),
+        "Omega_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'"
+    )
+
+
+    expect_warning(
+        a <- moma_fpca(X),
+        "No smoothness is imposed!"
+    )
+    expect_true(all(a$selection_scheme_list == c(0, 0, 0, 0)))
+
+    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,
+        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 wrappers: moma_twfpca", {
+    set.seed(12)
+
+    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_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_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"
+        ),
+        "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),
+            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_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."
+        ),
+        "Please use `moma_fpca` if only one side is penalized"
+    )
+
+    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'"
+        ),
+        "Please use `moma_fpca` if only one side is penalized"
+    )
+
+    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)))
+
+    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)))
+
+    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)))
+})
+
+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"
+    )
+})
+
+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), # selected by BIC
+        Omega_u = second_diff_mat(17), Omega_v = second_diff_mat(8),
+        selection_scheme_str = "gggb"
+    )
+    expect_error(
+        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
+    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)
+    )
+    # 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
+        ),
+        "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."
+    )
+
+
+    # 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."
+    )
+
+    # case 3: interpolation on v side, only lambda is a vector
+    a <- moma_sfpca(X,
+        u_sparsity = lasso(), v_sparsity = lasso(),
+        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."
+    )
+    # 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."
+    )
+})
+
+
+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 because both alpha_v and lambda_v should be specified
+    expect_error(
+        a$interpolate(alpha_v = 0.23),
+        "Please spesify the following argument(s): lambda_v.",
+        fixed = TRUE
+    )
+    expect_error(
+        a$interpolate(),
+        "Please spesify the following argument(s): alpha_v, lambda_v.",
+        fixed = TRUE
+    )
+
+    # error because alpha_u is a scalar during initialization
+    expect_error(
+        a$interpolate(alpha_u = 0.2323),
+        "Invalid index in SFPCA::interpolate: alpha_u"
+    )
+
+    # error because alpha_u should not be specified
+    expect_error(
+        a$interpolate(alpha_v = 0.23, lambda_v = 0.121, alpha_u = 1),
+        "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
+    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"
+    )
+
+    # 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"
+    )
+})