R/PCA.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

## File for code on Principal Component Analysis

#' Principal Component Analysis for Distributional Data
#'
#' @description
#' Performs Simplicial Functional Principal Component Analysis (SFPCA) on
#' distributional data where each observation is an entire probability distribution.
#' This function extends the classical `princomp` to work with BData objects,
#' enabling dimensionality reduction and exploratory analysis of datasets of
#' distributions.
#'
#' @param x A BData object containing the distributional data. Can be any subclass
#'     such as BDataSplines, BDataDensities, BDataPDFs, BDataHistograms, or
#'     BDataSamples. The data will be internally converted to spline coefficients
#'     for the analysis.
#' @param cor Logical. If `TRUE`, the analysis uses the correlation matrix instead
#'     of the covariance matrix. This standardizes each coefficient to unit variance
#'     before analysis, which can be useful when coefficients have very different
#'     scales. Default is `FALSE`.
#' @param scores Logical. If `TRUE` (default), the principal component scores are
#'     computed for each observation. Scores represent the coordinates of each
#'     distribution in the reduced PC space.
#' @param covmat An optional pre-computed covariance matrix. If provided, this
#'     matrix is used instead of computing the covariance from `x`. Useful for
#'     robust or regularized covariance estimates. **Note:** When provided with a
#'     non-orthonormal basis, this should be the covariance in the transformed
#'     (orthonormalized) space.
#' @param center An optional centering vector. If `NULL` (default), the column
#'     means of the coefficient matrix are used. Providing a custom center allows
#'     for centering around a specific reference distribution.
#' @param group An optional factor or vector specifying group membership for each
#'     observation. If provided, PCA is performed separately for each group level
#'     (unless `multivariate = TRUE`), and the result is an object of class
#'     `"princomp.BData.grouped"` containing a list of PCA results. This is useful
#'     for comparing principal components across different subgroups (e.g., males
#'     vs females, treatment vs control).
#' @param multivariate Logical. Only relevant when `group` is provided. If `FALSE`
#'     (default), PCA is performed separately for each group. If `TRUE`, performs
#'     **multivariate SFPCA** (2D-SFPCA per Hron et al. 2016) where clr-representations
#'     are stacked across groups and a single joint PCA is performed. This captures
#'     the covariation between groups and yields common eigenfunctions spanning all
#'     groups. Requires that observations are **paired** across groups (same number
#'     of observations per group, same observation index = same subject).
#' @param ... Additional arguments passed to methods.
#'
#' @return An object of class `c("princomp.BData", "princomp")` containing:
#' \describe{
#'   \item{sdev}{Numeric vector of standard deviations of the principal components,
#'       i.e., the square roots of the eigenvalues of the covariance matrix.
#'       These quantify the amount of variation captured by each component.}
#'   \item{loadings}{Matrix of eigenfunctions in coefficient form (class `"loadings"`).
#'       Each column is an eigenvector representing a principal component direction
#'       in the space of spline coefficients. Rows correspond to basis coefficients.
#'       **Note:** For non-orthonormal bases, these are back-transformed to the
#'       original coefficient space. Kept for compatibility with `stats::princomp`.}
#'   \item{eigenfunctions}{Alias for `loadings`. This is the proper SFPCA terminology
#'       (Hron et al. 2016).}
#'   \item{center}{Numeric vector containing the mean of each spline coefficient,
#'       used for centering the data before PCA.
#'       This represents the "average distribution" in coefficient space.}
#'   \item{scale}{Numeric vector of scalings applied to each coefficient. All ones
#'       unless `cor=TRUE`, in which case these are the standard deviations used
#'       for standardization.}
#'   \item{n.obs}{Integer giving the number of observations (distributions) in
#'       the dataset.}
#'   \item{scores}{Matrix of principal component scores (if `scores=TRUE`).
#'       Each row corresponds to an observation, each column to a principal
#'       component. The score indicates how far each distribution deviates from
#'       the mean along that PC direction.}
#'   \item{call}{The matched function call.}
#'   \item{Center}{The mean distribution as a BDataSplines object. This can be
#'       plotted or converted to other representations (PDF, CDF, etc.).}
#'   \item{Loadings}{The eigenfunctions as a BDataSplines object. Each row
#'       represents one principal component as a functional perturbation.
#'       Kept for compatibility.}
#'   \item{Eigenfunctions}{Alias for `Loadings`. This is the proper SFPCA terminology.}
#'   \item{basis}{The BayesSpaceBasis used for the spline representation.}
#'   \item{OrigClass}{Character string recording the original class of the input
#'       data, useful for back-conversion.}
#'   \item{orthonormal}{Logical indicating whether the basis was orthonormal
#'       (`TRUE`) or required Gram matrix correction (`FALSE`).}
#'   \item{GramMatrix}{The Gram matrix S from the basis. Identity matrix if
#'       the basis is orthonormal.}
#' }
#'
#' @details
#' ## Conceptual Background
#'
#' In distributional data analysis, each observation is not a single number but
#' an entire probability distribution. Examples include:
#' \itemize{
#'   \item Particle size distributions from different samples
#'   \item Income distributions across different regions
#'   \item Age distributions of populations over time
#'   \item Spectral distributions in chemistry or physics
#' }
#'
#' PCA for such data aims to find the main "modes of variation" - the principal
#' ways in which distributions in the dataset differ from each other.
#'
#' ## Method
#'
#' The function works by:
#' \enumerate{
#'   \item Converting all input distributions to a common spline basis representation
#'       using \code{\link{as.BDataSplines}}. This yields a matrix where each row
#'       is an observation and each column is a spline basis coefficient.
#'   \item **Handling non-orthonormal basis:** If the basis has a non-identity Gram
#'       matrix S, the coefficients are transformed via \eqn{C' = C \cdot S^{1/2}}{C' = C * S^(1/2)}
#'       to orthonormalize the representation before PCA.
#'   \item Computing the covariance (or correlation) matrix of these (transformed)
#'       coefficients.
#'   \item Performing eigendecomposition to find the principal directions.
#'   \item Projecting the centered data onto these directions to obtain scores.
#'   \item **Back-transforming loadings:** If transformation was applied, loadings
#'       are back-transformed via \eqn{L = S^{-1/2} \cdot L'}{L = S^(-1/2) * L'} so they
#'       remain in the original coefficient space.
#' }
#'
#' ## Non-Orthonormal Basis Handling
#'
#' **Important:** Many spline bases, including ZB-splines (zero-integral B-splines),
#' are **not orthonormal**. For such bases, standard PCA on raw coefficients gives
#' incorrect results because the Euclidean inner product in coefficient space does
#' not correspond to the L2 inner product in function space.
#'
#' This function automatically detects non-orthonormal bases by checking the Gram
#' matrix S stored in the \code{\link{BayesSpaceBasis}}. If \eqn{S \neq I}{S != I},
#' the function applies the correction described in Hron et al. (2016, Eq. 6):
#'
#' \deqn{N^{-1} C' C M b = \rho b}{N^(-1) C'CM b = rho * b}
#'

#' where M is the Gram matrix. This is implemented via Cholesky transformation:
#' \enumerate{
#'   \item Compute \eqn{S^{1/2}}{S^(1/2)} via Cholesky decomposition: \eqn{S = L L^T}{S = L * L^T}
#'   \item Transform coefficients: \eqn{\tilde{C} = C \cdot L}{C_tilde = C * L}
#'   \item Perform standard PCA on \eqn{\tilde{C}}{C_tilde}
#'   \item Back-transform loadings: \eqn{V = L^{-1} \cdot \tilde{V}}{V = L^(-1) * V_tilde}
#' }
#'
#' The scores remain unchanged since they are computed with the orthonormalized
#' representation. The loadings are transformed back to the original coefficient
#' space for interpretability.
#'
#' **For users:** If you want to use an orthonormal basis, you can either:
#' \itemize{
#'   \item Use an explicitly orthonormal spline basis
#'   \item Provide the correct Gram matrix S when constructing the
#'       \code{\link{BayesSpaceBasis}} (the function will handle the correction)
#' }
#'
#' The Gram matrix S defines the inner product in coefficient space:
#' \eqn{\langle c_1, c_2 \rangle_S = c_1^T S c_2}{<c1, c2>_S = c1^T * S * c2}.
#' If S = I (identity), the basis is orthonormal and no correction is needed.
#'
#' This approach is analogous to \code{\link[compositions]{princomp.acomp}} from

#' the compositions package, which performs PCA on compositional data using
#' log-ratio coordinates.
#'
#' ## Interpretation
#'
#' \itemize{
#'   \item \strong{Eigenvalues/sdev}: The variance explained by each component.
#'       Large values indicate important modes of variation.
#'   \item \strong{Eigenfunctions}: Directions in coefficient space representing
#'       functional perturbations. When added to/subtracted from the mean, they
#'       show how distributions change along each PC.
#'   \item \strong{Scores}: Coordinates of each observation in the PC space.
#'       Similar scores indicate similar distribution shapes.
#' }
#'
#' Use \code{\link{plot.princomp.BData}} with `type="eigenfunctions"` to visualize
#' how the distribution shape changes along each principal component direction.
#'
#' ## Important Considerations
#'
#' \itemize{
#'   \item The number of meaningful components is limited by both the number of
#'       observations and the basis dimension.
#'   \item If your data shows that PC1 explains nearly 100\% of variance, the
#'       distributions may vary primarily in one characteristic (e.g., only in
#'       scale or only in location). This is not a bug but a property of your data.
#'   \item For meaningful PCA, ensure your dataset has genuine multivariate
#'       variation. Use \code{\link{B2GaussianDistribution}} with an appropriate
#'       covariance structure for simulation studies.
#' }
#'
#' @seealso
#' \code{\link{plot.princomp.BData}} for visualization options including
#'     screeplots, score plots, biplots, and functional loadings plots.
#'
#' \code{\link{biplot.princomp.BData}} for creating biplots showing both
#'     observations and loadings.
#'
#' \code{\link{predict.princomp.BData}} for projecting new distributions onto
#'     the principal components.
#'
#' \code{\link{summary.princomp.BData}} for a summary of variance explained.
#'
#' \code{\link{as.BDataSplines}} for the conversion to spline representation.
#'
#' \code{\link[compositions]{princomp.acomp}} for the analogous method for
#'     compositional data in the compositions package.
#'
#' @references
#' Hron, K., Menafoglio, A., Templ, M., Hruzova, K., & Filzmoser, P. (2016).
#' Simplicial principal component analysis for density functions in Bayes spaces.
#' \emph{Computational Statistics & Data Analysis}, 94, 330-350.
#'
#' @author Matthias Templ
#' @family PCA
#' @export
#' @importFrom stats princomp
#' @importFrom Matrix bdiag
#'
#' @examples
#' ## Example 1: Basic PCA on simulated distributional data
#' ## -------------------------------------------------------
#' ## We use B2GaussianDistribution to generate distributions with
#' ## genuine multivariate variation (not just scaling differences).
#'
#' # Define a Gaussian process on Bayes space with exponential covariance
#' meanDist <- BDataPDFs(list(dunif), basis = B2Unif01Basis7)
#' covFunction <- function(x, y) exp(-4 * abs(x - y)^1.7)
#' gaussProc <- B2GaussianDistribution(meanDist, covFunction)
#'
#' # Simulate 30 distributions from this process
#' simData <- simulate(gaussProc, nsim = 30, seed = 42)
#' splineData <- as.BDataSplines(simData)
#'
#' # Perform PCA
#' pcaResult <- princomp(splineData)
#'
#' # View results
#' pcaResult
#' summary(pcaResult)
#'
#' ## Example 2: Visualizing and interpreting the PCA results
#' ## ---------------------------------------------------------
#'
#' # Screeplot: How much variance does each PC explain?
#' plot(pcaResult, type = "screeplot")
#' # Interpretation: The first two components (PC1 and PC2) capture most of the
#' # variance (around 70%), with a clear "elbow" after PC3. This suggests that
#' # the main patterns in the distributional data can be described by 2-3 modes
#' # of variation.
#'
#' # Score plot: How do observations cluster in PC space?
#' plot(pcaResult, type = "scores")
#' # Interpretation: Each point represents one distribution. Distributions that
#' # are close together in this plot have similar shapes. Look for clusters
#' # (groups of similar distributions), outliers (isolated points), or gradients
#' # (continuous variation along a direction).
#'
#' # To identify individual distributions, use labels:
#' plot(pcaResult, type = "scores", labels = TRUE)
#' # Now each distribution is labeled with its observation number, making it
#' # easy to trace outliers or interesting cases back to the original data.
#'
#' # Loadings plot: How do distributions change along each PC?
#' # This is THE KEY PLOT for interpretation!
#' plot(pcaResult, type = "loadings", npcs = 2)
#' # Interpretation: The black line shows the mean distribution. For each PC,
#' # the solid colored line shows what a distribution looks like when it has a
#' # HIGH score on that PC (positive direction), while the dashed line shows a
#' # LOW score (negative direction). For example, if PC1+ is shifted right and
#' # PC1- is shifted left compared to the mean, then PC1 captures location
#' # variation. If PC1+ is wider and PC1- is narrower, then PC1 captures spread.
#'
#' # Biplot: Combine scores and loadings in one plot
#' biplot(pcaResult)
#' # Interpretation: Numbers show observations, red arrows show how spline
#' # coefficients contribute to the PCs. Observations in similar positions have
#' # similar distributions. The arrows are harder to interpret directly since
#' # they represent abstract basis coefficients - use the loadings plot above
#' # for functional interpretation.
#'
#' ## Example 3: Predicting scores for new data
#' ## -------------------------------------------
#' # Generate new distributions from the same process
#' newData <- as.BDataSplines(simulate(gaussProc, nsim = 5, seed = 123))
#'
#' # Project onto the existing PC space
#' newScores <- predict(pcaResult, newData)
#' print(newScores)
#'
#' ## Example 4: PCA on real age distribution data
#' ## ----------------------------------------------
#' # Load age distribution data from Upper Austrian communes
#' data("age")
#'
#' # Create a basis for age distributions (support 0-100 years)
#' rmAge <- ReferenceMeasureUniform(0, 100)
#' B2Age <- BayesSpace(rmAge)
#' basisAge <- BayesSpaceBasis(B2Age, nbasis = 7)
#'
#' # Convert proportions to BDataDensityHistograms
#' # The data are proportions in 5-year age bins (last bin is 10 years: 90-100)
#' breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 100)
#' ageHist <- BDataDensityHistograms(age$distributions, breaks = breaks, basis = basisAge)
#'
#' # Convert: histograms -> densities -> splines
#' ageDensities <- as.BDataDensities(ageHist)
#' ageSplines <- as.BDataSplines(ageDensities)
#'
#' # Perform PCA
#' pcaAge <- princomp(ageSplines)
#' summary(pcaAge)
#'
#' # Screeplot: Most variance is in PC1 and PC2
#' plot(pcaAge, type = "screeplot")
#'
#' # Score plot colored by sex
#' sexColors <- ifelse(age$data$sex == "male", "blue", "red")
#' plot(pcaAge, type = "scores", col = sexColors)
#' legend("topright", c("Male", "Female"), col = c("blue", "red"), pch = 19)
#' # Interpretation: Males and females separate along PC1, indicating that the
#' # primary mode of variation captures sex differences in age structure
#' # (e.g., females tend to live longer, creating different elderly proportions).
#'
#' # Loadings plot: What do the PCs represent?
#' plot(pcaAge, type = "loadings", npcs = 2)
#' # Interpretation: PC1 captures the balance between younger and older ages.
#' # PC2 may capture differences in the working-age population bulge.
#'
#' ## Example 5: Grouped PCA - compare males vs females
#' ## ---------------------------------------------------
#' # Perform PCA separately for each sex group
#' pcaGrouped <- princomp(ageSplines, group = age$data$sex)
#' pcaGrouped
#'
#' # Compare screeplots: Do males and females have similar dimensionality?
#' plot(pcaGrouped, type = "screeplot")
#'
#' # Compare loadings: Are the modes of variation the same?
#' plot(pcaGrouped, type = "loadings", npcs = 2)
#' # Interpretation: If the loadings look similar across groups, the same
#' # types of variation exist in both. If they differ, each group has
#' # distinct patterns of age distribution variation.
#'
#' # Compare score distributions with boxplots
#' plot(pcaGrouped, type = "boxplot", npcs = 3)
#'
#' ## Example 6: Multivariate SFPCA (2D-SFPCA) - joint analysis across groups
#' ## -------------------------------------------------------------------------
#' # Unlike grouped PCA which analyzes groups separately, multivariate SFPCA
#' # stacks the representations and finds COMMON principal directions across groups.
#' # This requires paired data: same observation index = same subject across groups.
#'
#' pcaMulti <- princomp(ageSplines, group = age$data$sex, multivariate = TRUE)
#' pcaMulti
#'
#' # Screeplot: Total variance across both male and female distributions
#' plot(pcaMulti, type = "screeplot")
#'
#' # Score plot: Each point is a commune, scored on joint PCs
#' plot(pcaMulti, type = "scores")
#'
#' # Eigenfunctions: Shows how each group component of the joint eigenfunction
#' # perturbs the mean. Left column = male, right column = female.
#' plot(pcaMulti, type = "eigenfunctions", npcs = 2)
#' # Interpretation: If PC1 shows similar patterns for both groups, they vary
#' # together. If patterns differ, the PC captures group-specific variation
#' # or interactions between groups.
#'
princomp.BData <- function(x, cor = FALSE, scores = TRUE, covmat = NULL,
                           center = NULL, group = NULL, multivariate = FALSE, ...) {
    cl <- match.call()
    cl[[1L]] <- as.name("princomp")

    # Handle grouped PCA
    if (!is.null(group)) {
        return(.princompGrouped(x, group = group, cor = cor, scores = scores,
                                covmat = covmat, center = center,
                                multivariate = multivariate, ...))
    }

    # Store original class for later conversion back
    origClass <- class(x)[1] # TODO: exploit common interface gsi.orig(x)

    # Get basis
    basis <- getBayesSpaceBasis(x)

    # Convert to spline representation
    xSplines <- as.BDataSplines(x, basis = basis)
    xMat <- unclass(xSplines)

    n <- nrow(xMat)
    p <- ncol(xMat)

    # Get the Gram matrix S from the basis
    # S defines the inner product: <c1, c2>_S = c1^T S c2
    # If S = I (identity), the basis is orthonormal
    S <- basis$S
    if (is.null(S)) {
        S <- diag(p)
    }

    # Check if basis is orthonormal (S is identity matrix)
    isOrthonormal <- isTRUE(all.equal(S, diag(nrow(S)),
                                       check.attributes = FALSE,
                                       tolerance = .Machine$double.eps^0.5))

    # For non-orthonormal basis, we need to transform coefficients
    # to account for the Gram matrix (Hron et al. 2016, Eq. 6)
    # Transform: C_tilde = C %*% S^(1/2)
    # Then standard PCA on C_tilde, then back-transform loadings
    if (!isOrthonormal) {
        # Compute S^(1/2) via Cholesky decomposition: S = L %*% t(L)
        # So S^(1/2) = L (lower triangular Cholesky factor)
        cholS <- tryCatch(
            chol(S),
            error = function(e) {
                # If Cholesky fails, S might not be positive definite
                # Fall back to eigendecomposition
                eigS <- eigen(S, symmetric = TRUE)
                if (any(eigS$values <= 0)) {
                    stop("Gram matrix S is not positive definite. ",
                         "Cannot perform PCA with this basis.")
                }
                eigS$vectors %*% diag(sqrt(eigS$values)) %*% t(eigS$vectors)
            }
        )
        # cholS is upper triangular from R's chol(), so S = t(cholS) %*% cholS
        # We use L = t(cholS) as the "square root"
        sqrtS <- t(cholS)
        sqrtS_inv <- solve(sqrtS)

        # Transform coefficient matrix: each row c_i -> c_i %*% sqrtS
        xMatTransformed <- xMat %*% sqrtS

        message("Note: Basis is non-orthonormal (Gram matrix S != I). ",
                "Applying transformation for correct SFPCA (see Hron et al. 2016).")
    } else {
        xMatTransformed <- xMat
        sqrtS <- diag(p)
        sqrtS_inv <- diag(p)
    }

    # Compute center if not provided (in original space)
    if (is.null(center)) {
        center <- colMeans(xMat)
    }

    # Center the transformed data
    if (!isOrthonormal) {
        centerTransformed <- as.vector(center %*% sqrtS)
        xCenteredTransformed <- sweep(xMatTransformed, 2, centerTransformed, "-")
    } else {
        xCenteredTransformed <- sweep(xMatTransformed, 2, center, "-")
    }

    # Compute covariance matrix if not provided (on transformed data)
    if (is.null(covmat)) {
        if (cor) {
            covmat <- stats::cor(xMatTransformed)
        } else {
            covmat <- stats::cov(xMatTransformed)
        }
    }

    # Eigendecomposition
    eig <- eigen(covmat, symmetric = TRUE)

    # Remove components with near-zero eigenvalues
    tol <- max(dim(covmat)) * max(abs(eig$values)) * .Machine$double.eps
    keep <- eig$values > tol
    if (!any(keep)) {
        stop("All eigenvalues are effectively zero")
    }

    # Standard deviations (sqrt of eigenvalues)
    sdev <- sqrt(eig$values[keep])
    names(sdev) <- paste0("Comp.", seq_along(sdev))

    # Loadings matrix (in transformed space)
    loadingsTransformed <- eig$vectors[, keep, drop = FALSE]

    # Back-transform loadings to original coefficient space with the
    # INVERSE TRANSPOSE of the whitening factor: with S = M %*% t(M) and
    # whitened coefficients C %*% M, a whitened direction v_tilde
    # corresponds to the original-coordinate loading v = t(M)^{-1} v_tilde,
    # which satisfies <c, v>_S = c_tilde %*% v_tilde and
    # t(v) %*% S %*% v = 1 (S-orthonormal eigenfunctions).
    # M^{-1} v_tilde is wrong for the non-symmetric Cholesky factor.
    if (!isOrthonormal) {
        loadings <- t(sqrtS_inv) %*% loadingsTransformed
    } else {
        loadings <- loadingsTransformed
    }
    dimnames(loadings) <- list(paste0("Coef.", seq_len(p)),
                               paste0("Comp.", seq_along(sdev)))
    class(loadings) <- "loadings"

    # Compute scores if requested
    # Scores are computed in transformed space: scores = C_centered_transformed %*% V_transformed
    # This is equivalent to scores in original space with Gram matrix correction
    if (scores) {
        sc <- xCenteredTransformed %*% loadingsTransformed
        dimnames(sc) <- list(rownames(xMat), colnames(loadings))
    } else {
        sc <- NULL
    }

    # Create the mean as a BData object (Center) - in original space
    Center <- BDataSplines(matrix(center, nrow = 1), basis = basis)

    # Create loadings as BDataSplines (each column is a loading)
    # Loadings are in original coefficient space
    LoadingsBData <- BDataSplines(t(loadings), basis = basis)

    # Scale if correlation was used
    if (cor) {
        scale <- apply(xMatTransformed, 2, stats::sd)
    } else {
        scale <- rep(1, p)
    }
    names(scale) <- paste0("Coef.", seq_len(p))

    # Construct result
    # Note: "loadings" kept for compatibility with stats::princomp
    # "eigenfunctions" is the proper SFPCA terminology (Hron et al. 2016)
    result <- list(
        sdev = sdev,
        loadings = loadings,
        eigenfunctions = loadings,
        center = center,
        scale = scale,
        n.obs = n,
        scores = sc,
        call = cl,
        Center = Center,
        Loadings = LoadingsBData,
        Eigenfunctions = LoadingsBData,
        basis = basis,
        OrigClass = origClass,
        orthonormal = isOrthonormal,
        GramMatrix = S
    )
    attr(result, "orig") = x
    class(result) <- c("princomp.BData", "princomp")
    return(result)
}


# Internal function: Grouped PCA
# Routes to either separate PCA per group or multivariate (joint) SFPCA
.princompGrouped <- function(x, group, cor = FALSE, scores = TRUE,
                             covmat = NULL, center = NULL,
                             multivariate = FALSE, ...) {
    if (multivariate) {
        return(.princompMultivariate(x, group = group, cor = cor,
                                     scores = scores, covmat = covmat,
                                     center = center, ...))
    }

    # Get coefficient matrix and basis
    xMat <- unclass(x)
    basis <- getBayesSpaceBasis(x)

    if (length(group) != nrow(xMat)) {
        stop("Length of 'group' must match number of observations in 'x'")
    }

    group <- as.factor(group)
    levels <- levels(group)

    # Perform PCA for each group separately
    results <- lapply(levels, function(lvl) {
        idx <- which(group == lvl)
        # Subset the matrix and recreate BDataSplines
        xSubMat <- xMat[idx, , drop = FALSE]
        xSub <- BDataSplines(xSubMat, basis = basis)
        princomp.BData(xSub, cor = cor, scores = scores,
                       covmat = covmat, center = center, group = NULL, ...)
    })
    names(results) <- levels

    # Store group information
    structure(
        results,
        class = c("princomp.BData.grouped", "list"),
        group = group,
        levels = levels,
        call = match.call()
    )
}


# Internal function: Multivariate SFPCA (2D-SFPCA)
# Performs joint PCA on stacked clr-representations across groups
# Following the approach in Hron et al. (2016) and Menafoglio's implementation
#
# Key difference from separate grouped PCA:
# - Stacks coefficients: for observation i, [coefs_group1, coefs_group2, ...]
# - Single covariance on stacked (n_groups * n_coefs) dimensional space
# - Eigenfunctions span all groups jointly
# - Captures covariation BETWEEN groups
#
# Requires paired observations (same indices = same subject across groups)
.princompMultivariate <- function(x, group, cor = FALSE, scores = TRUE,
                                   covmat = NULL, center = NULL, ...) {
    cl <- match.call()
    cl[[1L]] <- as.name("princomp")

    # Get coefficient matrix and basis
    xMat <- unclass(x)
    basis <- getBayesSpaceBasis(x)
    origClass <- class(x)[1]  # TODO: exploit common interface gsi.orig(x)
    p <- ncol(xMat)  # Number of coefficients per distribution

    if (length(group) != nrow(xMat)) {
        stop("Length of 'group' must match number of observations in 'x'")
    }

    group <- as.factor(group)
    levels <- levels(group)
    nGroups <- length(levels)

    # Get indices for each group
    groupIndices <- lapply(levels, function(lvl) which(group == lvl))
    groupSizes <- sapply(groupIndices, length)

    # Check that all groups have the same number of observations (paired data)
    if (length(unique(groupSizes)) != 1) {
        stop("Multivariate SFPCA requires paired observations: ",
             "all groups must have the same number of observations.\n",
             "Group sizes: ", paste(levels, "=", groupSizes, collapse = ", "))
    }
    nObs <- groupSizes[1]

    # Get the Gram matrix S from the basis
    S <- basis$S
    if (is.null(S)) {
        S <- diag(p)
    }

    # Check if basis is orthonormal
    isOrthonormal <- isTRUE(all.equal(S, diag(nrow(S)),
                                       check.attributes = FALSE,
                                       tolerance = .Machine$double.eps^0.5))

    # For non-orthonormal basis, compute transformation matrices
    if (!isOrthonormal) {
        cholS <- tryCatch(
            chol(S),
            error = function(e) {
                eigS <- eigen(S, symmetric = TRUE)
                if (any(eigS$values <= 0)) {
                    stop("Gram matrix S is not positive definite.")
                }
                eigS$vectors %*% diag(sqrt(eigS$values)) %*% t(eigS$vectors)
            }
        )
        sqrtS <- t(cholS)
        sqrtS_inv <- solve(sqrtS)

        # For stacked coefficients, the Gram matrix is block-diagonal
        # S_stacked = diag(S, S, ..., S) with nGroups blocks
        sqrtS_stacked <- as.matrix(Matrix::bdiag(rep(list(sqrtS), nGroups)))
        sqrtS_inv_stacked <- as.matrix(Matrix::bdiag(rep(list(sqrtS_inv), nGroups)))

        message("Note: Basis is non-orthonormal (Gram matrix S != I). ",
                "Applying transformation for correct multivariate SFPCA.")
    } else {
        sqrtS_stacked <- diag(p * nGroups)
        sqrtS_inv_stacked <- diag(p * nGroups)
    }

    # Stack coefficients: for each observation i, concatenate coefficients
    # from all groups [coefs_group1_i, coefs_group2_i, ...]
    # Result: nObs x (nGroups * p) matrix
    xStacked <- matrix(NA, nrow = nObs, ncol = nGroups * p)
    colnames(xStacked) <- paste0(rep(levels, each = p), ".Coef.", rep(1:p, nGroups))

    for (g in seq_along(levels)) {
        idx <- groupIndices[[g]]
        cols <- ((g - 1) * p + 1):(g * p)
        xStacked[, cols] <- xMat[idx, , drop = FALSE]
    }

    # Apply Gram matrix transformation if needed
    if (!isOrthonormal) {
        xStackedTransformed <- xStacked %*% sqrtS_stacked
    } else {
        xStackedTransformed <- xStacked
    }

    # Compute center if not provided
    if (is.null(center)) {
        centerStacked <- colMeans(xStacked)
    } else {
        # center should be provided as stacked vector
        centerStacked <- center
    }

    # Center the transformed data
    if (!isOrthonormal) {
        centerTransformed <- as.vector(centerStacked %*% sqrtS_stacked)
        xCenteredTransformed <- sweep(xStackedTransformed, 2, centerTransformed, "-")
    } else {
        xCenteredTransformed <- sweep(xStackedTransformed, 2, centerStacked, "-")
    }

    # Compute covariance matrix if not provided
    if (is.null(covmat)) {
        if (cor) {
            covmat <- stats::cor(xStackedTransformed)
        } else {
            covmat <- stats::cov(xStackedTransformed)
        }
    }

    # Eigendecomposition
    eig <- eigen(covmat, symmetric = TRUE)

    # Remove components with near-zero eigenvalues
    tol <- max(dim(covmat)) * max(abs(eig$values)) * .Machine$double.eps
    keep <- eig$values > tol
    if (!any(keep)) {
        stop("All eigenvalues are effectively zero")
    }

    # Standard deviations
    sdev <- sqrt(eig$values[keep])
    names(sdev) <- paste0("Comp.", seq_along(sdev))

    # Loadings (in transformed space)
    loadingsTransformed <- eig$vectors[, keep, drop = FALSE]

    # Normalize eigenfunctions so that L2 norm = 1
    # In transformed (orthonormal) space, this is just Euclidean norm
    # But we normalize per-group to match Alessandra's approach
    for (i in seq_len(ncol(loadingsTransformed))) {
        # L2 norm using trapezoidal approximation would require evaluation
        # For coefficient space, we use the Gram-weighted norm
        # In transformed space, it's just Euclidean
        norm_i <- sqrt(sum(loadingsTransformed[, i]^2))
        loadingsTransformed[, i] <- loadingsTransformed[, i] / norm_i
    }

    # Back-transform loadings to original coefficient space with the
    # inverse TRANSPOSE (see the comment in princomp.BData)
    if (!isOrthonormal) {
        loadings <- t(sqrtS_inv_stacked) %*% loadingsTransformed
    } else {
        loadings <- loadingsTransformed
    }
    dimnames(loadings) <- list(colnames(xStacked),
                               paste0("Comp.", seq_along(sdev)))
    class(loadings) <- "loadings"

    # Compute scores: project centered data onto eigenvectors
    if (scores) {
        sc <- xCenteredTransformed %*% loadingsTransformed
        dimnames(sc) <- list(paste0("Obs.", seq_len(nObs)), colnames(loadings))
    } else {
        sc <- NULL
    }

    # Create the mean as a list of BDataSplines objects (one per group)
    CenterList <- lapply(seq_along(levels), function(g) {
        cols <- ((g - 1) * p + 1):(g * p)
        BDataSplines(matrix(centerStacked[cols], nrow = 1), basis = basis)
    })
    names(CenterList) <- levels

    # Create loadings as list of BDataSplines objects (one per group)
    # Each eigenfunction has components for each group
    LoadingsList <- lapply(seq_along(levels), function(g) {
        cols <- ((g - 1) * p + 1):(g * p)
        BDataSplines(t(loadings[cols, , drop = FALSE]), basis = basis)
    })
    names(LoadingsList) <- levels

    # Scale if correlation was used
    if (cor) {
        scale <- apply(xStackedTransformed, 2, stats::sd)
    } else {
        scale <- rep(1, p * nGroups)
    }
    names(scale) <- colnames(xStacked)

    # Construct result
    result <- list(
        sdev = sdev,
        loadings = loadings,
        eigenfunctions = loadings,
        center = centerStacked,
        scale = scale,
        n.obs = nObs,
        scores = sc,
        call = cl,
        Center = CenterList,
        Loadings = LoadingsList,
        Eigenfunctions = LoadingsList,
        basis = basis,
        OrigClass = origClass,
        orthonormal = isOrthonormal,
        GramMatrix = S,
        # Multivariate-specific fields
        levels = levels,
        nGroups = nGroups,
        nCoefs = p
    )

    class(result) <- c("princomp.BData.multivariate", "princomp.BData", "princomp")
    return(result)
}


#' Print a grouped distributional PCA
#'
#' @param x a `princomp.BData.grouped` object (from [princomp.BData] with a `group`).
#' @param ... further arguments, currently unused.
#' @returns `x`, invisibly. Prints, per group, the number of observations and the
#'     share of variance explained by the first principal component.
#' @seealso [princomp.BData]
#' @export
#' @method print princomp.BData.grouped
#' @examples
#' \donttest{
#' data("age")
#' basisAge <- BayesSpaceBasis(BayesSpace(ReferenceMeasureUniform(0, 100)), nbasis = 7)
#' breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 100)
#' ageHist <- BDataDensityHistograms(age$distributions, breaks = breaks, basis = basisAge)
#' ageSplines <- as.BDataSplines(as.BDataDensities(ageHist))
#' pcaGrouped <- princomp(ageSplines, group = age$data$sex)
#' print(pcaGrouped)
#' }
print.princomp.BData.grouped <- function(x, ...) {
    levels <- attr(x, "levels")
    cat("Grouped PCA on distributional data\n")
    cat("Groups:", paste(levels, collapse = ", "), "\n\n")

    for (lvl in levels) {
        cat("--- Group:", lvl, "---\n")
        cat("  Observations:", x[[lvl]]$n.obs, "\n")
        vars <- x[[lvl]]$sdev^2
        pct1 <- round(100 * vars[1] / sum(vars), 1)
        cat("  PC1 explains:", pct1, "% of variance\n")
    }
    invisible(x)
}


#' Print a multivariate (2D-)SFPCA on distributional data
#'
#' @param x a `princomp.BData.multivariate` object (from [princomp.BData] with
#'     `multivariate = TRUE`).
#' @param ... further arguments, currently unused.
#' @returns `x`, invisibly. Prints the groups, observations and dimensions of the
#'     joint analysis together with the variance explained.
#' @seealso [princomp.BData]
#' @export
#' @method print princomp.BData.multivariate
#' @examples
#' \donttest{
#' data("age")
#' basisAge <- BayesSpaceBasis(BayesSpace(ReferenceMeasureUniform(0, 100)), nbasis = 7)
#' breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 100)
#' ageHist <- BDataDensityHistograms(age$distributions, breaks = breaks, basis = basisAge)
#' ageSplines <- as.BDataSplines(as.BDataDensities(ageHist))
#' pcaMulti <- princomp(ageSplines, group = age$data$sex, multivariate = TRUE)
#' print(pcaMulti)
#' }
print.princomp.BData.multivariate <- function(x, ...) {
    cat("Multivariate SFPCA (2D-SFPCA) on distributional data\n")
    cat("Groups:", paste(x$levels, collapse = ", "), "\n")
    cat("Observations per group:", x$n.obs, "\n")
    cat("Coefficients per group:", x$nCoefs, "\n")
    cat("Total dimensions:", x$nCoefs * x$nGroups, "\n")
    if (!x$orthonormal) {
        cat("Note: Non-orthonormal basis, Gram matrix correction applied\n")
    }
    cat("\n")

    # Variance explained
    vars <- x$sdev^2
    propVar <- vars / sum(vars)
    cumProp <- cumsum(propVar)

    nshow <- min(5, length(vars))
    cat("Importance of components:\n")
    cat(sprintf("  %-12s", ""), paste(sprintf("%-10s", paste0("PC", 1:nshow)), collapse = ""), "\n")
    cat(sprintf("  %-12s", "Std. Dev."),
        paste(sprintf("%-10.4f", x$sdev[1:nshow]), collapse = ""), "\n")
    cat(sprintf("  %-12s", "Prop. Var."),
        paste(sprintf("%-10.4f", propVar[1:nshow]), collapse = ""), "\n")
    cat(sprintf("  %-12s", "Cum. Prop."),
        paste(sprintf("%-10.4f", cumProp[1:nshow]), collapse = ""), "\n")

    invisible(x)
}


#' Plot Method for Multivariate SFPCA of Distributional Data
#'
#' @description
#' Creates visualization plots for multivariate SFPCA (2D-SFPCA) results.
#' Eigenfunctions span multiple groups and are shown with one panel per group.
#'
#' @param x An object of class `princomp.BData.multivariate` as returned by
#'     \code{\link{princomp.BData}} with `group` and `multivariate = TRUE`.
#' @param type Character string specifying the type of plot:
#'     \describe{
#'       \item{"screeplot"}{Variance explained by each component}
#'       \item{"scores"}{Scatter plot of PC scores}
#'       \item{"eigenfunctions"}{Eigenfunctions shown per group (default)}
#'       \item{"loadings"}{Alias for eigenfunctions}
#'     }
#' @param npcs Number of principal components to display. Default is 2.
#' @param pcs Which PCs to use for scores plot. Default is c(1, 2).
#' @param as Representation for eigenfunction plots: "PDF", "clr", "density", etc.
#' @param ... Additional graphical arguments.
#'
#' @return Invisibly returns the multivariate PCA object.
#'
#' @details
#' For multivariate SFPCA, eigenfunctions are joint across all groups. The
#' eigenfunctions plot shows how each group's component of the eigenfunction
#' affects the distribution shape.
#'
#' @seealso
#' \code{\link{princomp.BData}} for performing multivariate SFPCA.
#'
#' @author Matthias Templ
#' @family PCA
#' @export
#' @method plot princomp.BData.multivariate
#' @importFrom graphics par mtext layout
#'
#' @examples
#' \donttest{
#' data("age")
#' basisAge <- BayesSpaceBasis(BayesSpace(ReferenceMeasureUniform(0, 100)), nbasis = 7)
#' breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 100)
#' ageHist <- BDataDensityHistograms(age$distributions, breaks = breaks, basis = basisAge)
#' ageSplines <- as.BDataSplines(as.BDataDensities(ageHist))
#' pcaMulti <- princomp(ageSplines, group = age$data$sex, multivariate = TRUE)
#' plot(pcaMulti, type = "eigenfunctions", npcs = 2)
#' }
plot.princomp.BData.multivariate <- function(x,
                                              type = c("screeplot", "variance",
                                                       "scores", "eigenfunctions",
                                                       "loadings"),
                                              npcs = 2, pcs = c(1, 2),
                                              as = "PDF", ...) {
    type <- match.arg(type)

    # Alias loadings -> eigenfunctions
    if (type == "loadings") type <- "eigenfunctions"
    if (type == "variance") type <- "screeplot"

    levels <- x$levels
    nGroups <- x$nGroups

    if (type == "screeplot") {
        # Standard screeplot
        vars <- x$sdev^2
        propVar <- 100 * vars / sum(vars)
        nshow <- min(10, length(vars))

        graphics::barplot(propVar[1:nshow],
                          names.arg = paste0("PC", 1:nshow),
                          ylab = "% Variance Explained",
                          main = "Multivariate SFPCA Screeplot",
                          col = "steelblue", ...)
    } else if (type == "scores") {
        # Score scatter plot
        if (is.null(x$scores)) {
            stop("Scores not available. Rerun princomp with scores = TRUE.")
        }

        sc <- x$scores
        graphics::plot(sc[, pcs[1]], sc[, pcs[2]],
                       xlab = paste0("PC", pcs[1]),
                       ylab = paste0("PC", pcs[2]),
                       main = "Multivariate SFPCA Scores",
                       pch = 19, ...)
        graphics::abline(h = 0, v = 0, col = "gray", lty = 2)
    } else if (type == "eigenfunctions") {
        # Plot eigenfunctions for each group
        # Following the approach in .plotLoadingsBData
        oldpar <- graphics::par(no.readonly = TRUE)
        on.exit(graphics::par(oldpar))

        # Layout: npcs rows, nGroups columns
        graphics::par(mfrow = c(npcs, nGroups), mar = c(4, 4, 3, 1))

        basis <- x$basis

        for (pc in 1:npcs) {
            for (g in seq_along(levels)) {
                lvl <- levels[g]

                # Get mean coefficients for this group
                meanCoefs <- unclass(x$Center[[lvl]])[1, ]

                # Get eigenfunction coefficients for this group
                eigCoefs <- unclass(x$Eigenfunctions[[lvl]])[pc, ]

                # Scale by standard deviation
                sdPC <- x$sdev[pc]
                scaledEig <- sdPC * eigCoefs

                # Create BDataSplines for mean, plus, and minus
                meanBData <- BDataSplines(matrix(meanCoefs, nrow = 1), basis = basis)
                posCoefs <- meanCoefs + scaledEig
                negCoefs <- meanCoefs - scaledEig
                posBData <- BDataSplines(matrix(posCoefs, nrow = 1), basis = basis)
                negBData <- BDataSplines(matrix(negCoefs, nrow = 1), basis = basis)

                # Variance explained
                pct <- round(100 * x$sdev[pc]^2 / sum(x$sdev^2), 1)
                mainTitle <- paste0(lvl, ": PC", pc, " (", pct, "%)")

                # Plot using BDataSplines plot method; ... lets callers set
                # e.g. xlab/ylab on the panels
                if (as == "clr") {
                    plot(clr(meanBData), col = "black", lwd = 2, lty = 1,
                         main = mainTitle, ...)
                    plot(clr(posBData), add = TRUE, col = "blue", lwd = 1.5, lty = 1)
                    plot(clr(negBData), add = TRUE, col = "red", lwd = 1.5, lty = 2)
                } else {
                    plot(meanBData, as = as, close = TRUE, col = "black", lwd = 2,
                         lty = 1, main = mainTitle, ...)
                    plot(posBData, as = as, close = TRUE, add = TRUE, col = "blue",
                         lwd = 1.5, lty = 1)
                    plot(negBData, as = as, close = TRUE, add = TRUE, col = "red",
                         lwd = 1.5, lty = 2)
                }
            }
        }
    }

    invisible(x)
}


#' Plot Method for Grouped PCA of Distributional Data
#'
#' @description
#' Creates multi-panel comparison plots for grouped PCA results. Each group
#' is shown in a separate panel, allowing visual comparison of PCA results
#' across different subgroups of the data.
#'
#' @param x An object of class `princomp.BData.grouped` as returned by
#'     \code{\link{princomp.BData}} with a `group` argument.
#' @param type Character string specifying the type of plot. Same options as
#'     \code{\link{plot.princomp.BData}}: `"screeplot"`, `"variance"`,
#'     `"scores"`, `"loadings"`, `"covariance"`, or `"reconstruction"`.
#' @param mfrow Optional vector of c(nrow, ncol) for panel layout. If `NULL`,
#'     panels are arranged automatically in a single row.
#' @param main Character string for overall title. If `NULL`, group names
#'     are used as individual panel titles.
#' @param ... Additional arguments passed to \code{\link{plot.princomp.BData}}.
#'
#' @return Invisibly returns the grouped PCA object.
#'
#' @details
#' This function sets up a multi-panel plotting layout and calls
#' \code{\link{plot.princomp.BData}} for each group. The layout is automatically
#' determined based on the number of groups, but can be overridden using `mfrow`.
#'
#' Common use cases:
#' \itemize{
#'   \item Compare loadings across groups to see if different subgroups have
#'       different modes of variation.
#'   \item Compare screeplots to assess if dimensionality differs by group.
#'   \item Compare score plots to see internal structure within each group.
#' }
#'
#' @seealso
#' \code{\link{princomp.BData}} for performing grouped PCA.
#' \code{\link{plot.princomp.BData}} for single-group plot options.
#'
#' @author Matthias Templ
#' @family PCA
#' @export
#' @method plot princomp.BData.grouped
#' @importFrom graphics par mtext
#'
#' @examples
#' # Create grouped PCA using age distribution data
#' data("age")
#'
#' # Set up basis
#' rmAge <- ReferenceMeasureUniform(0, 100)
#' B2Age <- BayesSpace(rmAge)
#' basisAge <- BayesSpaceBasis(B2Age, nbasis = 7)
#'
#' # Create histograms and convert to splines
#' breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 100)
#' ageHist <- BDataDensityHistograms(age$distributions, breaks = breaks, basis = basisAge)
#' ageSplines <- as.BDataSplines(as.BDataDensities(ageHist))
#'
#' # Perform grouped PCA by sex
#' pcaGrouped <- princomp(ageSplines, group = age$data$sex)
#' pcaGrouped
#'
#' # Compare screeplots by group
#' plot(pcaGrouped, type = "screeplot")
#'
#' # Compare loadings by group
#' plot(pcaGrouped, type = "loadings", npcs = 2)
#'
#' # Custom layout (2 rows)
#' plot(pcaGrouped, type = "scores", mfrow = c(2, 1))
#'
plot.princomp.BData.grouped <- function(x, type = c("screeplot", "variance",
                                                     "scores", "eigenfunctions",
                                                     "loadings", "covariance",
                                                     "reconstruction", "boxplot"),
                                         mfrow = NULL, main = NULL, ...) {
    type <- match.arg(type)

    # "loadings" is kept as alias for "eigenfunctions" for backwards compatibility
    if (type == "loadings") type <- "eigenfunctions"

    levels <- attr(x, "levels")

    # Validate input
    if (is.null(levels) || length(levels) == 0) {
        stop("Invalid grouped PCA object: no 'levels' attribute found. ",
             "Use princomp() with 'group' argument to create grouped PCA.")
    }

    nGroups <- length(levels)

    # Set up panel layout
    oldpar <- par(no.readonly = TRUE)
    on.exit(par(oldpar))

    if (is.null(mfrow)) {
        # Automatic layout: single row
        mfrow <- c(1, nGroups)
    }
    par(mfrow = mfrow)

    # Plot each group
    for (i in seq_along(levels)) {
        lvl <- levels[i]
        groupMain <- if (is.null(main)) lvl else paste(main, "-", lvl)

        plot(x[[lvl]], type = type, main = groupMain, ...)
    }

    invisible(x)
}


#' @describeIn princomp.BData Print method displaying key PCA results
#' @param x A princomp.BData object to print.
#' @param ... Additional arguments (currently unused).
#' @export
#' @method print princomp.BData
print.princomp.BData <- function(x, ...) {
    cat("Call:\n")
    print(x$call)
    cat("\nStandard deviations:\n")
    print(x$sdev)
    vars <- x$sdev^2
    pct1 <- round(100 * vars[1] / sum(vars), 1)
    pct2 <- if (length(vars) > 1) round(100 * sum(vars[1:min(2, length(vars))]) / sum(vars), 1) else pct1
    cat("\nPC1 explains", pct1, "% of the variance")
    if (length(vars) > 1) {
        cat(", PC1-2 explain", pct2, "%")
    }
    cat("\nBasis dimension:", ncol(x$loadings), "\n")
    cat("Number of observations:", x$n.obs, "\n")

    # Report on orthonormality
    if (!is.null(x$orthonormal)) {
        if (x$orthonormal) {
            cat("Basis: orthonormal (no Gram matrix correction needed)\n")
        } else {
            cat("Basis: non-orthonormal (Gram matrix correction applied)\n")
        }
    }

    invisible(x)
}


#' Plot Method for PCA of Distributional Data
#'
#' @description
#' A versatile plotting function for visualizing principal component analysis
#' results on distributional data. Provides multiple plot types for comprehensive
#' exploration of the PCA results, including the unique ability to visualize
#' loadings as functional perturbations of the mean distribution.
#'
#' @param x An object of class `princomp.BData` as returned by
#'     \code{\link{princomp.BData}}.
#' @param type Character string specifying the type of plot to produce:
#'   \describe{
#'     \item{`"screeplot"`}{Bar chart showing the percentage of variance explained
#'         by each principal component. Useful for deciding how many components
#'         to retain.}
#'     \item{`"variance"`}{Line plot showing cumulative variance explained as
#'         more components are included. Horizontal reference lines at 80\%, 90\%,
#'         and 95\% help identify cutoff points.}
#'     \item{`"scores"`}{Scatter plot of observations in the space of two
#'         principal components. Useful for identifying clusters, outliers, or
#'         patterns in the data.}
#'     \item{`"eigenfunctions"`}{\strong{Key plot for interpretation.} Shows how the
#'         distribution shape changes along each principal component direction.
#'         The mean distribution is shown along with perturbations in the positive
#'         and negative PC directions. This is the standard SFPCA terminology
#'         (Hron et al. 2016).}
#'     \item{`"loadings"`}{Alias for `"eigenfunctions"` (kept for backwards
#'         compatibility).}
#'     \item{`"biplot"`}{Combines scores and eigenfunctions in a single plot. Calls
#'         \code{\link{biplot.princomp.BData}}.}
#'     \item{`"covariance"`}{Heatmap showing the functional covariance surface.
#'         This visualizes which parts of the domain co-vary: high values indicate
#'         that when one part of the distribution is above average, the other part
#'         tends to be as well. Useful for understanding the correlation structure
#'         across the support of the distributions.}
#'     \item{`"reconstruction"`}{Shows the reconstructed distributions using only
#'         the first `npcs` principal components. Useful for assessing how well
#'         a reduced number of components approximates the original data.}
#'     \item{`"boxplot"`}{Boxplots showing the distribution of scores for each
#'         principal component. Useful for visualizing the spread of observations
#'         along each PC direction and identifying potential outliers.}
#'   }
#' @param choices Integer vector of length 2 specifying which principal components
#'     to plot (for `type="scores"` and `type="biplot"`). Default is `c(1, 2)`.
#' @param scale Scaling factor for loadings in biplot (passed to
#'     \code{\link{biplot.princomp.BData}}). Default is 1.
#' @param pc.biplot Logical. If `TRUE`, use Gabriel's biplot scaling. Default
#'     is `FALSE`.
#' @param main Character string for plot title. If `NULL`, an appropriate default
#'     is generated based on the plot type.
#' @param col Colors for plotting. For loadings plot: first color is for the mean,
#'     subsequent colors for each PC. For scores plot: color for points/labels.
#' @param labels For `type="scores"`: Controls labeling of observations:
#'   \describe{
#'     \item{`FALSE`}{(Default) Show only points, no labels.}
#'     \item{`TRUE`}{Show observation numbers (1, 2, 3, ...) instead of points.}
#'     \item{character vector}{Use custom labels for each observation.}
#'   }
#'   This helps identify which distribution each point represents.
#' @param lwd Line width for functional plots (loadings). Default is 1.
#' @param lty Line type (currently not directly used; solid for positive direction,
#'     dashed for negative).
#' @param legend Logical. If `TRUE` (default), add a legend to the loadings plot.
#' @param ... Additional graphical parameters passed to underlying plot functions.
#' @param npcs Number of principal components to display in screeplot, variance
#'     plot, or loadings plot. Default is `min(10, length(x$sdev))`.
#' @param as For loadings plot: how to display the distributions. One of:
#'   \describe{
#'     \item{`"PDF"`}{Probability density function (default). Shows distributions
#'         as proper densities that integrate to 1.}
#'     \item{`"density"`}{Density with respect to the reference measure.}
#'     \item{`"logPDF"`}{Log of the probability density function.}
#'     \item{`"logDensity"`}{Log density with respect to reference measure.}
#'     \item{`"clr"`}{Centered log-ratio transformation. Shows the additive
#'         structure in Bayes space. The clr is centered around zero, making
#'         it easy to see positive and negative deviations from the geometric
#'         mean. This is the natural representation for Bayes space PCA.}
#'     \item{`"coef"`}{Raw spline coefficients (parallel coordinates plot).
#'         Useful for seeing which coefficients drive each PC.}
#'   }
#' @param scale.sdev Numeric scaling factor for the loadings plot. Controls how
#'     far from the mean the perturbed distributions are shown:
#'   \describe{
#'     \item{`1`}{Show distributions at ±1 standard deviation from mean (default).}
#'     \item{`2`}{Show distributions at ±2 standard deviations (more extreme).}
#'     \item{`NA`}{Show unit loadings without scaling by standard deviation.}
#'   }
#'
#' @return Invisibly returns the data used for plotting (varies by plot type).
#'
#' @details
#' ## Understanding the Eigenfunctions Plot
#'
#' The `type="eigenfunctions"` plot is the most important visualization for
#' interpreting SFPCA results on distributional data. It answers the question:
#' "How does the distribution shape change as we move along each principal component?"
#'
#' The plot shows:
#' \itemize{
#'   \item \strong{Black thick line}: The mean distribution (center of the data).
#'   \item \strong{Colored solid lines}: Distributions at +`scale.sdev` standard
#'       deviations from the mean along each PC direction (positive perturbation).
#'   \item \strong{Colored dashed lines}: Distributions at -`scale.sdev` standard
#'       deviations from the mean along each PC direction (negative perturbation).
#' }
#'
#' For example, if PC1 captures variation in the "location" of distributions, the
#' PC1+ line will be shifted right and PC1- will be shifted left compared to the
#' mean. If PC1 captures variation in "spread", PC1+ will be wider and PC1- will
#' be narrower.
#'
#' ## Choosing the Number of Components
#'
#' Use the screeplot and variance plot to decide how many components are meaningful:
#' \itemize{
#'   \item \strong{Screeplot}: Look for an "elbow" where the variance drops off.
#'       Components before the elbow are typically retained.
#'   \item \strong{Variance plot}: Choose enough components to explain a target
#'       percentage (e.g., 80\% or 90\%) of total variance.
#' }
#'
#' @seealso
#' \code{\link{princomp.BData}} for performing the PCA.
#'
#' \code{\link{biplot.princomp.BData}} for detailed control over biplots.
#'
#' \code{\link{screeplot.princomp.BData}} for the standalone screeplot function.
#'
#' @author Matthias Templ
#' @family PCA
#' @export
#' @method plot princomp.BData
#' @importFrom graphics par barplot lines legend axis abline arrows text points
#'
#' @examples
#' ## Setup: Create PCA result for examples
#' meanDist <- BDataPDFs(list(dunif), basis = B2Unif01Basis7)
#' gaussProc <- B2GaussianDistribution(meanDist, function(x,y) exp(-4*abs(x-y)^1.7))
#' pcaResult <- princomp(as.BDataSplines(simulate(gaussProc, 30, seed = 42)))
#'
#' ## Screeplot: Decide how many components to keep
#' plot(pcaResult, type = "screeplot")
#' # Interpretation: The bars show the percentage of total variance explained by
#' # each PC. Here, PC1 and PC2 dominate, suggesting that the main structure in
#' # the data can be captured by just two components. The "elbow" after PC2 or
#' # PC3 indicates where additional components add little explanatory power.
#'
#' ## Cumulative variance: What percentage is explained?
#' plot(pcaResult, type = "variance")
#' # Interpretation: This shows how total explained variance accumulates as more
#' # components are included. The dashed reference lines help identify how many
#' # components are needed to reach common thresholds (80%, 90%, 95%). Here,
#' # approximately 3 components explain over 90% of the variance.
#'
#' ## Score plot: Explore structure in the data
#' plot(pcaResult, type = "scores")
#' # Interpretation: Each point is one distribution projected onto PC1 and PC2.
#' # The spread along PC1 (x-axis) shows variation in the first mode; spread
#' # along PC2 (y-axis) shows the second mode. Clusters indicate groups of
#' # similar distributions; outliers are distributions with unusual shapes.
#'
#' ## Score plot with labels: Identify individual distributions
#' plot(pcaResult, type = "scores", labels = TRUE)
#' # Using labels=TRUE replaces points with observation numbers (1, 2, 3, ...),
#' # allowing you to identify which specific distribution each position
#' # represents. This is useful for tracking outliers or interesting cases
#' # back to the original data.
#'
#' ## Eigenfunctions plot: THE KEY INTERPRETIVE PLOT
#' plot(pcaResult, type = "eigenfunctions", npcs = 2)
#' # Interpretation: The black curve is the mean (average) distribution. For
#' # each PC, solid lines show how distributions look when they score HIGH on
#' # that component, dashed lines show LOW scores. Compare the colored curves
#' # to the black mean to see what "high PC1" or "low PC2" means in terms of
#' # actual distribution shape (e.g., shifted location, changed spread, altered
#' # skewness, or different modality).
#'
#' ## More extreme perturbations (2 standard deviations)
#' plot(pcaResult, type = "eigenfunctions", npcs = 2, scale.sdev = 2)
#' # Using scale.sdev=2 shows more extreme deviations from the mean, making
#' # the differences between PC directions more visible.
#'
#' ## Show loadings as log-PDFs (useful for seeing tail behavior)
#' plot(pcaResult, type = "loadings", npcs = 2, as = "logPDF")
#' # The log-PDF representation emphasizes differences in the tails of the
#' # distributions, which can be important for detecting outlier behavior.
#'
#' ## Show loadings in CLR (centered log-ratio) representation
#' plot(pcaResult, type = "loadings", npcs = 2, as = "clr")
#' # The CLR representation shows the additive structure in Bayes space.
#' # Values are centered around zero, making it easy to see which parts of the
#' # distribution are above or below the geometric mean. This is the natural
#' # representation for PCA in Bayes space.
#'
#' ## Covariance surface: See which parts of the domain co-vary
#' plot(pcaResult, type = "covariance")
#' # Interpretation: This heatmap shows functional covariance Cov(f(s), f(t))
#' # across the domain. High values (red) indicate that when one part of the
#' # distribution is above average, the other tends to be as well. The diagonal
#' # shows variance at each point. Off-diagonal patterns reveal correlation
#' # structure (e.g., if left and right tails move together or oppositely).
#'
#' ## Reconstruction: How well do k PCs approximate the data?
#' plot(pcaResult, type = "reconstruction", npcs = 2)
#' # Shows all distributions reconstructed using only the first 2 PCs.
#' # Compare to the original data to assess information loss.
#' # If reconstructions look similar to originals, 2 PCs are sufficient.
#'
#' plot(pcaResult, type = "reconstruction", npcs = 1)
#' # Using only 1 PC shows more loss of detail - useful for understanding
#' # what information each additional PC captures.
#'
#' ## Boxplot: Distribution of scores across PCs
#' plot(pcaResult, type = "boxplot", npcs = 4)
#' # Shows the spread of observations along each PC direction.
#' # Useful for identifying outliers (points beyond whiskers) and comparing
#' # the amount of variation captured by each component. The percentage
#' # labels show variance explained by each PC.
#'
plot.princomp.BData <- function(x, type = c("screeplot", "variance", "scores",
                                            "eigenfunctions", "loadings", "biplot",
                                            "covariance", "reconstruction", "boxplot"),
                                choices = 1:2, scale = 1, pc.biplot = FALSE,
                                main = NULL, col = NULL, lwd = 1, lty = 1,
                                legend = TRUE, labels = FALSE, ...,
                                npcs = min(10, length(x$sdev)),
                                as = c("PDF", "density", "logPDF", "logDensity", "clr", "coef"),
                                scale.sdev = 1) {
    type <- match.arg(type)
    as <- match.arg(as)

    # "loadings" is kept as alias for "eigenfunctions" for backwards compatibility
    if (type == "loadings") type <- "eigenfunctions"

    switch(type,
           screeplot = .plotScreeBData(x, npcs = npcs, main = main, ...),
           variance = .plotVarianceBData(x, npcs = npcs, main = main, ...),
           scores = .plotScoresBData(x, choices = choices, main = main,
                                     col = col, labels = labels, ...),
           eigenfunctions = .plotLoadingsBData(x, npcs = npcs, main = main,
                                         col = col, lwd = lwd, lty = lty,
                                         legend = legend, as = as,
                                         scale.sdev = scale.sdev, ...),
           biplot = biplot.princomp.BData(x, choices = choices, scale = scale,
                                          pc.biplot = pc.biplot, main = main,
                                          col = col, ...),
           covariance = .plotCovarianceBData(x, main = main, col = col, ...),
           reconstruction = .plotReconstructionBData(x, npcs = npcs, main = main,
                                                     col = col, as = as, ...),
           boxplot = .plotBoxplotBData(x, npcs = npcs, main = main, col = col, ...)
    )
}


# Internal function: Screeplot
.plotScreeBData <- function(x, npcs, main = NULL, ...) {
    if (is.null(main)) main <- "Screeplot"
    vars <- x$sdev^2
    vars <- vars[1:min(npcs, length(vars))]
    pct <- 100 * vars / sum(x$sdev^2)

    bp <- barplot(pct, names.arg = paste0("PC", seq_along(pct)),
                  main = main, ylab = "Percentage of Variance",
                  xlab = "Principal Component", ...)
    lines(bp, pct, type = "b", pch = 19, col = "blue")
    invisible(list(variance = vars, percent = pct))
}


# Internal function: Cumulative variance plot
.plotVarianceBData <- function(x, npcs, main = NULL, ...) {
    if (is.null(main)) main <- "Cumulative Variance Explained"
    vars <- x$sdev^2
    vars <- vars[1:min(npcs, length(vars))]
    cumvar <- cumsum(vars) / sum(x$sdev^2) * 100

    plot(seq_along(cumvar), cumvar, type = "b", pch = 19,
         main = main, xlab = "Number of Components",
         ylab = "Cumulative Variance (%)",
         xaxt = "n", ylim = c(0, 100), ...)
    axis(1, at = seq_along(cumvar))
    abline(h = c(80, 90, 95), lty = 2, col = "gray60")
    invisible(cumvar)
}


# Internal function: Score plot
# labels: FALSE = points only, TRUE = row numbers, character vector = custom labels
.plotScoresBData <- function(x, choices, main = NULL, col = NULL,
                             labels = FALSE, cex = NULL, pch = 19, ...) {
    if (is.null(x$scores)) {
        stop("No scores available. Re-run princomp with scores=TRUE.")
    }
    if (is.null(main)) {
        main <- paste("Scores on PC", choices[1], "vs PC", choices[2])
    }
    if (is.null(col)) col <- 1

    scores <- x$scores[, choices, drop = FALSE]
    varexp <- round(100 * x$sdev[choices]^2 / sum(x$sdev^2), 1)

    # Determine labeling mode
    n <- nrow(scores)
    useLabels <- !isFALSE(labels)

    # Set up plot
    plot(scores[, 1], scores[, 2], type = "n",
         main = main,
         xlab = paste0("PC", choices[1], " (", varexp[1], "%)"),
         ylab = paste0("PC", choices[2], " (", varexp[2], "%)"))
    abline(h = 0, v = 0, lty = 2, col = "gray60")

    if (useLabels) {
        # Determine what labels to use
        if (isTRUE(labels)) {
            # Use row numbers or rownames if available
            lbls <- rownames(scores)
            if (is.null(lbls)) lbls <- seq_len(n)
        } else {
            # Custom labels provided
            lbls <- labels
            if (length(lbls) != n) {
                warning("Number of labels does not match number of observations. ",
                        "Recycling labels.")
                lbls <- rep_len(lbls, n)
            }
        }

        # Add text labels
        if (is.null(cex)) cex <- 0.8
        text(scores[, 1], scores[, 2], labels = lbls, col = col, cex = cex)
    } else {
        # Original behavior: points only (pch is recycled across observations,
        # so a vector pch adds a redundant, color-vision-safe group encoding)
        points(scores[, 1], scores[, 2], pch = pch, col = col, cex = cex)
    }

    invisible(scores)
}


# Internal function: Loadings as functional plot
# Shows how the distribution changes along each principal component direction.
# The mean distribution is shown in black (thick line), and for each PC:
# - solid colored line: mean + scale.sdev * sdev * loading (positive direction)
# - dashed colored line: mean - scale.sdev * sdev * loading (negative direction)
# This represents the distribution of an observation that is scale.sdev standard
# deviations away from the mean along that PC direction.
.plotLoadingsBData <- function(x, npcs, main = NULL, col = NULL,
                               lwd = 1, lty = 1, legend = TRUE,
                               as = "PDF", scale.sdev = 1, ...) {
    if (is.null(main)) main <- paste("Eigenfunctions: Mean \u00b1 PC perturbations as", as)

    npcs <- min(npcs, length(x$sdev))
    basis <- x$basis

    # Default colors: black for mean, then standard colors for PCs
    if (is.null(col)) {
        col <- c(1, 2:(npcs + 1))
    }
    # Recycle colors if needed
    if (length(col) < npcs + 1) {
        col <- rep(col, length.out = npcs + 1)
    }

    # Get the mean coefficients
    meanCoefs <- x$center

    # Get the loadings (eigenvectors, unit length in coefficient space)
    loadingsMat <- x$loadings[, 1:npcs, drop = FALSE]

    # Scale loadings by standard deviation if requested
    # scale.sdev = 1 means show +/- 1 standard deviation from mean
    # scale.sdev = 2 means show +/- 2 standard deviations from mean
    if (!is.na(scale.sdev) && scale.sdev != 0) {
        loadingsMat <- sweep(loadingsMat, 2, x$sdev[1:npcs] * scale.sdev, "*")
    }

    # Create the mean as BDataSplines
    meanBData <- BDataSplines(matrix(meanCoefs, nrow = 1), basis = basis)

    # Collect all distributions to determine plot range
    allBData <- list(meanBData)
    for (i in seq_len(npcs)) {
        posCoefs <- meanCoefs + loadingsMat[, i]
        negCoefs <- meanCoefs - loadingsMat[, i]
        allBData[[length(allBData) + 1]] <- BDataSplines(matrix(posCoefs, nrow = 1), basis = basis)
        allBData[[length(allBData) + 1]] <- BDataSplines(matrix(negCoefs, nrow = 1), basis = basis)
    }

    # Plot mean first (thick black line)
    # Handle "clr" separately since it's not a standard plot.BDataSplines option
    if (as == "clr") {
        plotArgs <- list(col = col[1], lwd = lwd * 2, lty = 1, main = main, ...)
        do.call(plot, c(list(x = clr(meanBData)), plotArgs))
    } else {
        plotArgs <- list(close = TRUE, col = col[1], lwd = lwd * 2, lty = 1,
                         main = main, ...)
        do.call(plot, c(list(x = meanBData, as = as), plotArgs))
    }

    # Build legend
    legendLabels <- "Mean"
    legendCols <- col[1]
    legendLty <- 1
    legendLwd <- lwd * 2

    # Add perturbations for each PC
    for (i in seq_len(npcs)) {
        # Positive direction: mean + loading
        posCoefs <- meanCoefs + loadingsMat[, i]
        posBData <- BDataSplines(matrix(posCoefs, nrow = 1), basis = basis)

        # Negative direction: mean - loading
        negCoefs <- meanCoefs - loadingsMat[, i]
        negBData <- BDataSplines(matrix(negCoefs, nrow = 1), basis = basis)

        # Plot positive direction (solid line)
        if (as == "clr") {
            plot(clr(posBData), add = TRUE, col = col[i + 1], lwd = lwd, lty = 1)
        } else {
            plot(posBData, as = as, close = TRUE, add = TRUE,
                 col = col[i + 1], lwd = lwd, lty = 1)
        }

        # Plot negative direction (dashed line)
        if (as == "clr") {
            plot(clr(negBData), add = TRUE, col = col[i + 1], lwd = lwd, lty = 2)
        } else {
            plot(negBData, as = as, close = TRUE, add = TRUE,
                 col = col[i + 1], lwd = lwd, lty = 2)
        }

        # Update legend
        pct <- round(100 * x$sdev[i]^2 / sum(x$sdev^2), 1)
        legendLabels <- c(legendLabels,
                          paste0("PC", i, "+ (", pct, "%)"),
                          paste0("PC", i, "-"))
        legendCols <- c(legendCols, col[i + 1], col[i + 1])
        legendLty <- c(legendLty, 1, 2)
        legendLwd <- c(legendLwd, lwd, lwd)
    }

    if (legend) {
        legend("topright", legend = legendLabels, col = legendCols,
               lty = legendLty, lwd = legendLwd, bty = "n", cex = 0.7)
    }

    invisible(list(mean = meanBData, loadings = loadingsMat))
}


# Internal function: Covariance surface plot
# Shows the covariance structure of the functional data across the domain.
# This helps identify which parts of the support co-vary.
.plotCovarianceBData <- function(x, main = NULL, col = NULL, nGrid = 101, ...) {
    if (is.null(main)) main <- "Functional Covariance Surface"

    basis <- x$basis
    support <- getSupport(basis)
    lims=getSupport1DPlottingLimits(support)

    # Create evaluation grid (Support uses $a and $b for continuous range)
    grid <- seq(lims[1], lims[2], length.out = nGrid)

    # Evaluate orthonormalized basis functions at grid points
    # basis$basis returns raw B-splines (nGrid x nRawBasis)
    # We need to apply the transformation matrix A to get orthonormal basis
    # B_ortho = B_raw %*% t(A), giving nGrid x nCoef
    B_raw <- basis$basis(grid)
    B <- B_raw %*% t(basis$A)

    # Reconstruct coefficient covariance from PCA
    # Cov_coef = Loadings %*% diag(sdev^2) %*% t(Loadings)
    loadings <- x$loadings
    vars <- x$sdev^2
    CovCoef <- loadings %*% diag(vars) %*% t(loadings)

    # Functional covariance: Cov_func = B %*% Cov_coef %*% t(B)
    CovFunc <- B %*% CovCoef %*% t(B)

    # Default color palette (blue-white-red)
    if (is.null(col)) {
        col <- grDevices::hcl.colors(100, "Blue-Red 3")
    }

    # Plot as image
    graphics::image(grid, grid, CovFunc,
                    col = col,
                    main = main,
                    xlab = "Domain",
                    ylab = "Domain",
                    ...)

    # Add contour lines
    graphics::contour(grid, grid, CovFunc, add = TRUE, col = "gray30", lwd = 0.5)

    invisible(list(grid = grid, covariance = CovFunc))
}


# Internal function: Reconstruction plot
# Shows the reconstructed distributions using a reduced number of PCs.
# Reconstruction: center + sum(loadings[k] * scores[i,k]) for k = 1 to npcs
.plotReconstructionBData <- function(x, npcs, main = NULL, col = NULL,
                                     as = "PDF", ...) {
    if (is.null(x$scores)) {
        stop("No scores available. Re-run princomp with scores=TRUE.")
    }

    npcs <- min(npcs, length(x$sdev))
    if (is.null(main)) {
        main <- paste("Reconstruction using", npcs, "PCs")
    }

    basis <- x$basis
    n <- nrow(x$scores)

    # Default colors
    if (is.null(col)) {
        col <- grDevices::hcl.colors(n, "Set2")
    }
    if (length(col) < n) {
        col <- rep(col, length.out = n)
    }

    # Reconstruct coefficients: center + scores %*% t(loadings)
    # Using only first npcs components
    reconstructedCoefs <- matrix(x$center, nrow = n, ncol = length(x$center),
                                  byrow = TRUE)
    reconstructedCoefs <- reconstructedCoefs +
        x$scores[, 1:npcs, drop = FALSE] %*% t(x$loadings[, 1:npcs, drop = FALSE])

    # Create BDataSplines from reconstructed coefficients
    reconstructedBData <- BDataSplines(reconstructedCoefs, basis = basis)

    # Plot the reconstructed distributions
    if (as == "clr") {
        plot(clr(reconstructedBData), col = col, main = main, ...)
    } else {
        plot(reconstructedBData, as = as, close = TRUE, col = col, main = main, ...)
    }

    invisible(reconstructedBData)
}


# Internal function: Boxplot of score distributions
# Shows the distribution of scores for each principal component as boxplots.
# This helps visualize the spread and identify outliers in PC space.
.plotBoxplotBData <- function(x, npcs, main = NULL, col = NULL, ...) {
    if (is.null(x$scores)) {
        stop("No scores available. Re-run princomp with scores=TRUE.")
    }

    npcs <- min(npcs, length(x$sdev))
    if (is.null(main)) {
        main <- "Distribution of PC Scores"
    }

    # Get scores for the first npcs components
    scores <- x$scores[, 1:npcs, drop = FALSE]

    # Compute variance explained for labels
    varexp <- round(100 * x$sdev[1:npcs]^2 / sum(x$sdev^2), 1)
    labels <- paste0("PC", 1:npcs, "\n(", varexp, "%)")

    # Default colors
    if (is.null(col)) {
        col <- grDevices::hcl.colors(npcs, "Set2")
    }
    if (length(col) < npcs) {
        col <- rep(col, length.out = npcs)
    }

    # Create boxplot
    graphics::boxplot(as.data.frame(scores), names = labels, main = main,
                      col = col, ylab = "Score", ...)

    # Add horizontal line at zero
    graphics::abline(h = 0, lty = 2, col = "gray60")

    invisible(scores)
}


#' Biplot for PCA of Distributional Data
#'
#' @description
#' Creates a biplot that simultaneously displays the principal component scores
#' (observations) and loadings (variables/coefficients) in the same plot. This
#' visualization helps understand relationships between observations and the
#' underlying coefficient structure.
#'
#' @param x An object of class `princomp.BData` as returned by
#'     \code{\link{princomp.BData}}.
#' @param choices Integer vector of length 2 specifying which principal components
#'     to display. Default is `c(1, 2)` for the first two components.
#' @param scale Numeric scaling factor controlling the relative scaling of scores
#'     and loadings:
#'   \describe{
#'     \item{`0`}{(Default) Distance biplot: Loadings unscaled, scores scaled.
#'         Distances between observation labels approximate Euclidean distances
#'         in the full space.
#'         Best for seeing relationships between observations.}
#'     \item{`1`}{Correlation biplot: Loadings scaled by standard deviations.
#'         Inner products between loading arrows approximate covariances.
#'         Best for seeing relationships between variables.}
#'   }
#' @param pc.biplot Logical. If `TRUE`, use Gabriel's (1971) biplot scaling,
#'     which sets `scale=1`. Default is `FALSE`.
#' @param main Character string for plot title. Default is `"Biplot"`.
#' @param col Vector of two colors: first for observation labels, second for
#'     loading arrows. Default is `c("black", "red")`.
#' @param cex Vector of two character expansion factors: first for observation
#'     labels, second for loading labels. Default is `c(0.8, 0.8)`.
#' @param xlabs Character vector of labels for observations. Default uses
#'     observation numbers (1, 2, 3, ...).
#' @param ylabs Character vector of labels for loading arrows (coefficients).
#'     Default uses `"C1"`, `"C2"`, etc.
#' @param expand Numeric expansion factor for loading arrows. Increase this
#'     if the arrows are too short to see clearly. Default is 1.
#' @param arrow.len Length of arrow heads. Default is 0.1.
#' @param ... Additional graphical parameters passed to the underlying plot.
#'
#' @return Invisibly returns a list containing the scaled `scores` and `loadings`
#'     matrices as plotted.
#'
#' @details
#' ## Biplot Interpretation
#'
#' In a biplot:
#' \itemize{
#'   \item \strong{Observations} (numbers/text) represent individual distributions.
#'       Observations close together have similar distribution shapes.
#'   \item \strong{Loading arrows} (red arrows labeled C1, C2, etc.) represent
#'       spline basis coefficients. The direction and length indicate how each
#'       coefficient contributes to the principal components.
#'   \item \strong{Projecting} an observation onto a loading arrow approximates
#'       that observation's value on that coefficient (relative to the mean).
#'   \item \strong{Arrow angles}: Arrows pointing in similar directions indicate
#'       positively correlated coefficients; opposite directions indicate negative
#'       correlation; perpendicular arrows indicate uncorrelated coefficients.
#' }
#'
#' ## Note on Coefficient Interpretation
#'
#' Unlike PCA on raw data where variables have direct meaning, the "variables"
#' here are spline basis coefficients. These are abstract and their individual
#' interpretation may not be intuitive. For understanding how distributions
#' vary, the \code{\link{plot.princomp.BData}} function with `type="eigenfunctions"`
#' provides a more interpretable visualization showing how distribution shapes
#' change along each PC direction.
#'
#' ## Why Biplot Works for ZB-Splines
#'
#' The biplot is meaningful for the default ZB-spline basis because ZB-splines
#' have **local support**: each basis function is non-zero only in a specific
#' region of the domain. This means coefficient i primarily reflects the behavior
#' of the distribution near knot region i. When a loading arrow for coefficient
#' Ci points strongly in a PC direction, it indicates that variation in that
#' region of the domain is associated with that principal component.
#'
#' **Important:** If orthonormal splines were used instead, the biplot would
#' lose this interpretability because orthonormal basis functions have global
#' support (they are non-zero across the entire domain), so each coefficient
#' would mix information from everywhere rather than representing a local region.
#'
#' @seealso
#' \code{\link{princomp.BData}} for performing the PCA.
#'
#' \code{\link{plot.princomp.BData}} for other visualization options, especially
#'     `type="loadings"` for functional interpretation.
#'
#' \code{\link[stats]{biplot}} for the generic biplot function.
#'
#' @references
#' Gabriel, K. R. (1971). The biplot graphic display of matrices with application
#' to principal component analysis. \emph{Biometrika}, 58(3), 453-467.
#'
#' @author Matthias Templ
#' @family PCA
#' @export
#' @method biplot princomp.BData
#' @importFrom graphics arrows text points
#'
#' @examples
#' ## Create example PCA result
#' meanDist <- BDataPDFs(list(dunif), basis = B2Unif01Basis7)
#' gaussProc <- B2GaussianDistribution(meanDist, function(x,y) exp(-4*abs(x-y)^1.7))
#' pcaResult <- princomp(as.BDataSplines(simulate(gaussProc, 25, seed = 42)))
#'
#' ## Basic biplot
#' biplot(pcaResult)
#'
#' ## Biplot of PC2 vs PC3
#' biplot(pcaResult, choices = c(2, 3))
#'
#' ## Correlation biplot (scale=1)
#' biplot(pcaResult, scale = 1)
#'
#' ## Customized biplot
#' biplot(pcaResult,
#'        col = c("blue", "darkgreen"),
#'        cex = c(1, 0.7),
#'        expand = 1.5,
#'        main = "PCA Biplot of Distributional Data")
#'
biplot.princomp.BData <- function(x, choices = 1:2, scale = 0, pc.biplot = FALSE,
                                  main = NULL, col = c("black", "red"),
                                  cex = c(0.8, 0.8),
                                  xlabs = NULL, ylabs = NULL,
                                  expand = 1, arrow.len = 0.1, ...) {
    if (length(choices) != 2) {
        stop("choices must be a vector of length 2")
    }
    if (is.null(x$scores)) {
        stop("No scores available. Re-run princomp with scores=TRUE.")
    }

    if (is.null(main)) {
        main <- "Biplot"
    }

    # Get scores and loadings for chosen components
    scores <- x$scores[, choices, drop = FALSE]
    loadings <- x$loadings[, choices, drop = FALSE]

    # Apply scaling (following standard biplot convention)
    # scale=0: loadings unscaled, scores scaled by 1/sdev (distance biplot)
    # scale=1: loadings scaled by sdev, scores unscaled (correlation biplot)
    lam <- x$sdev[choices]
    if (pc.biplot) scale <- 1

    n <- x$n.obs
    lam_scale <- lam^scale

    # Scale scores and loadings
    scores <- t(t(scores) / lam^(1 - scale)) / sqrt(n)
    loadings <- t(t(loadings) * lam_scale)

    # Auto-scale loadings to be comparable to scores
    # This makes the biplot more readable
    score_range <- max(abs(scores))
    loading_range <- max(abs(loadings))
    if (loading_range > 0) {
        auto_scale <- score_range / loading_range * 0.8
    } else {
        auto_scale <- 1
    }
    loadings <- loadings * auto_scale * expand

    # Labels
    if (is.null(xlabs)) {
        xlabs <- seq_len(nrow(scores))
    }
    if (is.null(ylabs)) {
        ylabs <- rownames(loadings)
        if (is.null(ylabs)) ylabs <- paste0("C", seq_len(nrow(loadings)))
    }

    # Compute axis limits
    rangx <- range(scores[, 1], loadings[, 1])
    rangy <- range(scores[, 2], loadings[, 2])

    # Variance explained
    varexp <- round(100 * x$sdev[choices]^2 / sum(x$sdev^2), 1)

    # Plot
    plot(rangx, rangy, type = "n",
         main = main,
         xlab = paste0("PC", choices[1], " (", varexp[1], "%)"),
         ylab = paste0("PC", choices[2], " (", varexp[2], "%)"),
         asp = 1, ...)

    # Add grid lines
    abline(h = 0, v = 0, lty = 2, col = "gray70")

    # Plot observation scores
    if (length(col) < 2) col <- rep(col, 2)
    if (length(cex) < 2) cex <- rep(cex, 2)

    text(scores[, 1], scores[, 2], labels = xlabs, col = col[1], cex = cex[1])

    # Plot loadings as arrows
    for (i in seq_len(nrow(loadings))) {
        arrows(0, 0, loadings[i, 1], loadings[i, 2],
               col = col[2], length = arrow.len)
        text(loadings[i, 1] * 1.1, loadings[i, 2] * 1.1,
             labels = ylabs[i], col = col[2], cex = cex[2])
    }

    invisible(list(scores = scores, loadings = loadings))
}


#' Predict Scores for New Distributional Data
#'
#' @description
#' Projects new distributional data onto the principal component space defined
#' by a fitted PCA model. This allows new distributions to be positioned relative
#' to the original dataset without re-running the entire PCA.
#'
#' @param object A `princomp.BData` object containing the fitted PCA model.
#' @param newdata A BData object containing the new distributions to project.
#'     Can be any BData subclass; it will be converted to splines using the
#'     same basis as the original PCA.
#' @param ... Additional arguments (currently unused).
#'
#' @return A matrix of principal component scores with observations in rows and
#'     components in columns. Each row contains the coordinates of one new
#'     distribution in the PC space.
#'
#' @details
#' The prediction process:
#' \enumerate{
#'   \item Converts `newdata` to spline coefficients using the same basis as the
#'       original PCA.
#'   \item **For non-orthonormal bases:** Applies the same Gram matrix transformation
#'       that was used during fitting (see \code{\link{princomp.BData}} details).
#'   \item Centers the (transformed) coefficients using the center from the original PCA.
#'   \item Projects the centered coefficients onto the loading vectors to obtain
#'       scores.
#' }
#'
#' The resulting scores can be compared directly to the scores from the original
#' analysis. New observations with scores similar to original observations have
#' similar distribution shapes.
#'
#' **Note:** For non-orthonormal bases, the function automatically applies the
#' same transformation that was used during fitting. This ensures consistency
#' between training and prediction.
#'
#' @seealso
#' \code{\link{princomp.BData}} for fitting the PCA model.
#'
#' \code{\link{plot.princomp.BData}} with `type="scores"` for visualizing scores.
#'
#' @author Matthias Templ
#' @family PCA
#' @export
#' @method predict princomp.BData
#'
#' @examples
#' ## Fit PCA on training data
#' meanDist <- BDataPDFs(list(dunif), basis = B2Unif01Basis7)
#' gaussProc <- B2GaussianDistribution(meanDist, function(x,y) exp(-4*abs(x-y)^1.7))
#' trainData <- as.BDataSplines(simulate(gaussProc, 30, seed = 42))
#' pcaModel <- princomp(trainData)
#'
#' ## Generate new data and predict scores
#' newData <- as.BDataSplines(simulate(gaussProc, 5, seed = 999))
#' newScores <- predict(pcaModel, newData)
#' print(newScores)
#'
#' ## Compare: new observations vs training data
#' trainScores <- pcaModel$scores
#' plot(trainScores[,1], trainScores[,2], pch = 19, col = "gray50",
#'      xlab = "PC1", ylab = "PC2", main = "Training (gray) vs New (red)")
#' points(newScores[,1], newScores[,2], pch = 17, col = "red", cex = 1.5)
#' legend("topright", c("Training", "New"), pch = c(19, 17),
#'        col = c("gray50", "red"))
#'
predict.princomp.BData <- function(object, newdata, ...) {
    if (missing(newdata)) {
        return(object$scores)
    }else{
      stopifnot(hasSameBayesSpace(attr(object,"orig"),newdata)) # TODO_ exploit gsi.orig
    }

    # Convert new data to spline representation
    newSplines <- as.BDataSplines(newdata, basis = getBayesSpaceBasis(attr(object,"orig")))
    newMat <- unclass(newSplines)

    # For non-orthonormal basis, we need to transform the new data
    # the same way as the training data
    if (!isTRUE(object$orthonormal)) {
        S <- object$GramMatrix
        cholS <- chol(S)
        sqrtS <- t(cholS)

        # Transform new data
        newMatTransformed <- newMat %*% sqrtS

        # Center using the transformed center
        centerTransformed <- as.vector(object$center %*% sqrtS)
        newCenteredTransformed <- sweep(newMatTransformed, 2, centerTransformed, "-")

        # For scores, we need the loadings in transformed space
        # loadings (original) = t(sqrtS_inv) %*% loadingsTransformed
        # So: loadingsTransformed = t(sqrtS) %*% loadings
        loadingsTransformed <- t(sqrtS) %*% object$loadings

        # Project onto transformed loadings
        scores <- newCenteredTransformed %*% loadingsTransformed
    } else {
        # Orthonormal case: simple projection
        newCentered <- sweep(newMat, 2, object$center, "-")
        scores <- newCentered %*% object$loadings
    }

    return(scores)
}


#' Summary of PCA for Distributional Data
#'
#' @description
#' Provides a summary of the principal component analysis results, showing the
#' importance of each component in terms of standard deviation, proportion of
#' variance explained, and cumulative proportion.
#'
#' @param object A `princomp.BData` object.
#' @param ... Additional arguments (currently unused).
#'
#' @return An object of class `summary.princomp.BData` containing:
#' \describe{
#'   \item{sdev}{Standard deviations of the principal components.}
#'   \item{importance}{Matrix with three rows: standard deviation, proportion
#'       of variance, and cumulative proportion for each component.}
#' }
#'
#' @seealso \code{\link{princomp.BData}}
#' @author Matthias Templ
#' @family PCA
#' @export
#' @method summary princomp.BData
#'
#' @examples
#' meanDist <- BDataPDFs(list(dunif), basis = B2Unif01Basis7)
#' gaussProc <- B2GaussianDistribution(meanDist, function(x,y) exp(-4*abs(x-y)^1.7))
#' pcaResult <- princomp(as.BDataSplines(simulate(gaussProc, 25, seed = 42)))
#' summary(pcaResult)
#'
summary.princomp.BData <- function(object, ...) {
    vars <- object$sdev^2
    result <- list(
        sdev = object$sdev,
        importance = rbind(
            "Standard deviation" = object$sdev,
            "Proportion of Variance" = vars / sum(vars),
            "Cumulative Proportion" = cumsum(vars) / sum(vars)
        )
    )
    class(result) <- "summary.princomp.BData"
    result
}


#' @export
print.summary.princomp.BData <- function(x, digits = 3, ...) {
    cat("Importance of components:\n")
    print(round(x$importance, digits))
    invisible(x)
}


#' Screeplot for PCA of Distributional Data
#'
#' @description
#' Creates a screeplot showing the variance (or eigenvalue) associated with each
#' principal component. This is a standard tool for deciding how many principal
#' components to retain based on the "elbow" criterion.
#'
#' @param x A `princomp.BData` object.
#' @param npcs Number of principal components to include in the plot. Default
#'     is `min(10, length(x$sdev))`.
#' @param type Character string specifying the plot type:
#'   \describe{
#'     \item{`"barplot"`}{Bar chart of variances (default).}
#'     \item{`"lines"`}{Line plot connecting variance values.}
#'   }
#' @param main Plot title. Default is `"Screeplot"`.
#' @param ... Additional graphical parameters.
#'
#' @return Invisibly returns the variance values for each component.
#'
#' @details
#' The screeplot displays the eigenvalues (variances) of the principal components
#' in decreasing order. The name comes from the plot's resemblance to a scree
#' slope at the base of a mountain.
#'
#' ## How to Use
#'
#' Look for an "elbow" in the plot where the variance drops off sharply.
#' Components before the elbow are typically considered important and retained,
#' while those after the elbow contribute little to the overall structure.
#'
#' Alternative criteria for choosing components include:
#' \itemize{
#'   \item Retain components explaining a target percentage of variance (see
#'       `plot(x, type="variance")`).
#'   \item Kaiser criterion: retain components with variance > 1 (when using
#'       correlation matrix, `cor=TRUE`).
#' }
#'
#' @seealso
#' \code{\link{princomp.BData}} for the PCA analysis.
#'
#' \code{\link{plot.princomp.BData}} for additional plot types including
#'     cumulative variance plots.
#'
#' @author Matthias Templ
#' @family PCA
#' @export
#' @method screeplot princomp.BData
#' @importFrom stats screeplot
#'
#' @examples
#' meanDist <- BDataPDFs(list(dunif), basis = B2Unif01Basis7)
#' gaussProc <- B2GaussianDistribution(meanDist, function(x,y) exp(-4*abs(x-y)^1.7))
#' pcaResult <- princomp(as.BDataSplines(simulate(gaussProc, 25, seed = 42)))
#'
#' # Standard barplot screeplot
#' screeplot(pcaResult)
#'
#' # Line plot version
#' screeplot(pcaResult, type = "lines")
#'
screeplot.princomp.BData <- function(x, npcs = min(10, length(x$sdev)),
                                     type = c("barplot", "lines"),
                                     main = "Screeplot", ...) {
    type <- match.arg(type)

    vars <- x$sdev^2
    vars <- vars[seq_len(min(npcs, length(vars)))]

    if (type == "barplot") {
        barplot(vars, names.arg = paste0("PC", seq_along(vars)),
                main = main, ylab = "Variances", ...)
    } else {
        plot(seq_along(vars), vars, type = "b", pch = 19,
             main = main, xlab = "Component", ylab = "Variances",
             xaxt = "n", ...)
        axis(1, at = seq_along(vars))
    }

    invisible(vars)
}