R/RcppExports.R

In EFAtools: Fast and Flexible Implementations of Exploratory Factor Analysis Tools

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

.fit_ml_cpp <- function(R, n_fac, start, lower) {
    .Call(`_EFAtools_fit_ml_cpp`, R, n_fac, start, lower)
}

.error_ml <- function(psi, R, n_fac) {
    .Call(`_EFAtools_error_ml`, psi, R, n_fac)
}

.grad_ml <- function(psi, R, n_fac) {
    .Call(`_EFAtools_grad_ml`, psi, R, n_fac)
}

.fit_uls_cpp <- function(R, n_fac, start, lower) {
    .Call(`_EFAtools_fit_uls_cpp`, R, n_fac, start, lower)
}

.uls_residuals <- function(psi, R, n_fac) {
    .Call(`_EFAtools_uls_residuals`, psi, R, n_fac)
}

.grad_uls <- function(psi, R, n_fac) {
    .Call(`_EFAtools_grad_uls`, psi, R, n_fac)
}

.fit_dwls_cpp <- function(R, n_fac, W) {
    .Call(`_EFAtools_fit_dwls_cpp`, R, n_fac, W)
}

.dwls_residuals <- function(par, R, n_fac, W) {
    .Call(`_EFAtools_dwls_residuals`, par, R, n_fac, W)
}

.grad_dwls <- function(par, R, n_fac, W) {
    .Call(`_EFAtools_grad_dwls`, par, R, n_fac, W)
}

#' Compute number of non-matching indicator-to-factor correspondences
#'
#' @param x numeric matrix. A matrix of pattern coefficients.
#' @param y numeric matrix. A second matrix of coefficients.
#' @param thresh numeric. The threshold to classify a pattern coefficient as substantial.
#'
#' @keywords internal
.factor_corres <- function(x, y, thresh = 0.3) {
    .Call(`_EFAtools_factor_corres`, x, y, thresh)
}

#' Oblique Procrustes target rotation using a k x k inner objective
#'
#' Compute an oblique target rotation for a loading matrix using a
#' `targetQ`-compatible parameterization and a `k x k` objective.
#'
#' The rotated loading matrix is defined as
#' `L = A %*% solve(t(T))`, and the corresponding factor correlation matrix is
#' `Phi = t(T) %*% T`. The optimization is carried out over the transformation
#' matrix `T` under the oblique normalization constraint `diag(t(T) %*% T) = 1`.
#'
#' Non-invertible candidate transformations are rejected rather than evaluated
#' through a pseudo-inverse.
#'
#' Additional random starts may be requested. To reduce runtime, the solver uses
#' a two-stage strategy for extra starts: cheap objective screening, followed by
#' short triage optimization, followed by full optimization only for starts that
#' improve on the current incumbent by at least `triage_improve_tol`.
#'
#' @param A Numeric matrix. Loading matrix to be rotated.
#' @param B Numeric matrix. Target loading matrix with the same dimensions as
#'   `A`.
#' @param S_r Optional numeric `k x k` matrix containing `crossprod(A)`.
#'   Supplying this is useful when the same `A` is rotated repeatedly. Ignored
#'   when `normalize = TRUE` because the normalized cross-product is different.
#' @param T_init_r Optional numeric `k x k` starting transformation matrix.
#'   If `NULL`, the identity matrix is used for the primary start.
#' @param eps Numeric scalar. Convergence tolerance for the projected-gradient
#'   norm.
#' @param maxit Integer scalar. Maximum number of full projected-gradient
#'   updates.
#' @param max_line_search Integer scalar. Maximum number of step-halving
#'   attempts after the initial trial step in each line-search phase.
#' @param step0 Numeric scalar. Initial step size used in the projected-gradient
#'   update.
#' @param normalize Logical scalar. If `TRUE`, apply Kaiser normalization to the
#'   loadings (only) before rotation and reverse it afterwards; the target is left
#'   unnormalized, matching `GPArotation::targetQ(normalize = TRUE)`.
#' @param random_starts Integer scalar. Number of additional random starts.
#' @param screen_keep Integer scalar. Number of screened random starts retained
#'   for triage optimization.
#' @param triage_maxit Integer scalar. Number of short optimization iterations
#'   used in the triage stage.
#' @param triage_improve_tol Numeric scalar. Relative improvement required for a
#'   triaged start to be promoted to full optimization.
#'
#' @returns A named list containing the rotated loadings, transformation matrix,
#'   factor correlation matrix, target criterion value, convergence diagnostics,
#'   line-search diagnostics, and multi-start summaries.
#'
#' @details
#' The routine is intended for repeated oblique target rotations in workflows
#' such as bootstrap alignment or consensus alignment of exploratory factor
#' solutions across multiply imputed datasets. It follows the same oblique
#' transformation convention as `GPArotation::targetQ()`.
#'
#' @references
#' Bernaards, C. A., & Jennrich, R. I. (2005). Gradient projection algorithms
#' and software for arbitrary rotation criteria in factor analysis.
#' *Educational and Psychological Measurement*, 65, 676-696.
#'
#' Browne, M. W. (2001). An overview of analytic rotation in exploratory factor
#' analysis. *Multivariate Behavioral Research*, 36, 111-150.
#'
#' @keywords internal
#'
.oblique_procrustes <- function(A, B, S_r = NULL, T_init_r = NULL, eps = 1e-5, maxit = 1000L, max_line_search = 10L, step0 = 1.0, normalize = FALSE, random_starts = 0L, screen_keep = 2L, triage_maxit = 25L, triage_improve_tol = 0.0) {
    .Call(`_EFAtools_oblique_procrustes`, A, B, S_r, T_init_r, eps, maxit, max_line_search, step0, normalize, random_starts, screen_keep, triage_maxit, triage_improve_tol)
}

#' Batched oblique Procrustes target rotation over a cube of loading matrices
#'
#' Align each slice of a loading-matrix cube to a single shared target using the
#' same oblique target rotation as `.oblique_procrustes()`, in one call. This
#' removes the per-replicate marshalling overhead of looping `efa_procrustes()` in R
#' over bootstrap or multiple-imputation arrays.
#'
#' Each slice `A[, , i]` is aligned to `B`. For a single-factor cube the alignment
#' reduces to the closed-form sign match `T = sign(crossprod(A_i, B))` with factor
#' correlation `1`, matching the one-factor short-circuit in `efa_procrustes()`. For
#' two or more factors the slice is warm-started from the closed-form orthogonal
#' Procrustes solution (mirroring `efa_procrustes()`) and optimized with the same
#' multi-start oblique solver as `.oblique_procrustes()`. Random starts are drawn
#' serially with `R::rnorm` in the calling process.
#'
#' Slices are aligned independently. A slice that cannot be aligned (a non-finite
#' loading matrix, a failed warm-start decomposition, an invalid fit, or any
#' linear-algebra exception) is reported with `valid = FALSE` and `NA` for the
#' loadings, factor correlations, and all other per-slice diagnostics, rather than
#' aborting the whole call, so one degenerate replicate does not discard the rest.
#'
#' @param A Numeric array of dimension `n x m x b`: the `b` loading matrices to
#'   align.
#' @param B Numeric `n x m` target loading matrix shared across all slices.
#' @param eps Numeric scalar. Convergence tolerance for the projected-gradient
#'   norm.
#' @param maxit Integer scalar. Maximum number of full projected-gradient updates.
#' @param max_line_search Integer scalar. Maximum number of step-halving attempts
#'   after the initial trial step in each line-search phase.
#' @param step0 Numeric scalar. Initial step size used in the projected-gradient
#'   update.
#' @param normalize Logical scalar. If `TRUE`, apply Kaiser normalization to the
#'   loadings (only) before rotation and reverse it afterwards, leaving the target
#'   unnormalized (ignored for single-factor slices).
#' @param random_starts Integer scalar. Number of additional random starts per
#'   slice.
#' @param screen_keep Integer scalar. Number of screened random starts retained
#'   for triage optimization.
#' @param triage_maxit Integer scalar. Number of short optimization iterations
#'   used in the triage stage.
#' @param triage_improve_tol Numeric scalar. Relative improvement required for a
#'   triaged start to be promoted to full optimization.
#'
#' @returns A named list with the aligned-loadings array `loadings` (`n x m x b`),
#'   the factor-correlation array `Phi` (`m x m x b`), and the per-slice
#'   diagnostics `valid`, `convergence`, `value`, `iterations`, and
#'   `line_search_failed`.
#'
#' @keywords internal
#'
.oblique_procrustes_batch <- function(A, B, eps = 1e-5, maxit = 1000L, max_line_search = 10L, step0 = 1.0, normalize = FALSE, random_starts = 0L, screen_keep = 2L, triage_maxit = 25L, triage_improve_tol = 0.0) {
    .Call(`_EFAtools_oblique_procrustes_batch`, A, B, eps, maxit, max_line_search, step0, normalize, random_starts, screen_keep, triage_maxit, triage_improve_tol)
}

#' Perform the iterative PAF procedure
#'
#' Function called from within PAF so usually no call to this is needed by the user.
#' Provides a C++ implementation of the PAF procedure
#'
#' @param h2 numeric. The initial communality estimates.
#' @param criterion double. The convergence criterion to use.
#' @param R matrix. The correlation matrix with the initial communality estimates in the diagonal.
#' @param n_fac numeric. The number of factors to extract.
#' @param abs_eig logical. Whether absolute eigenvalues should be used to compute the loadings.
#' @param crit_type numeric. Whether maximum absolute differences (crit_type = 1), or sum of differences (crit_type = 2) should be used
#' @param max_iter numeric. The number of iterations after which to end the procedure if no convergence has been reached by then.
#'
#' @keywords internal
.paf_iter <- function(h2, criterion, R, n_fac, abs_eig, crit_type, max_iter) {
    .Call(`_EFAtools_paf_iter`, h2, criterion, R, n_fac, abs_eig, crit_type, max_iter)
}

#' Parallel analysis on simulated data.
#'
#' Function called from within efa_parallel() so usually no call to this is needed by the user.
#' Provides a C++ implementation of the efa_parallel() simulation procedure
#'
#' @param n_datasets numeric. Number of datasets with dimensions (N, n_vars) to simulate.
#' @param n_vars numeric. Number of variables / indicators in dataset.
#' @param N numeric. Number of cases / observations in dataset.
#' @param eigen_type numeric. Whether PCA (eigen_type = 1; i.e., leaving diagonal of correlation matrix at 1), SMC (eigen_type = 2; i.e., setting diagonal of correlation matrix to SMCs), or both from the same simulated datasets (eigen_type = 3), in which case the returned matrix holds the PCA eigenvalues in the first n_vars columns and the SMC eigenvalues in the next n_vars.
#' @param maxit numeric. Maximum iterations to perform after which to abort.
#'
#' @keywords internal
.parallel_sim <- function(n_datasets, n_vars, N, eigen_type, maxit = 10000L) {
    .Call(`_EFAtools_parallel_sim`, n_datasets, n_vars, N, eigen_type, maxit)
}

.polychoric_cpp <- function(x, acov, nearest_pd) {
    .Call(`_EFAtools_polychoric_cpp`, x, acov, nearest_pd)
}

.bvn_rect_cpp <- function(a0, a1, b0, b1, rho) {
    .Call(`_EFAtools_bvn_rect_cpp`, a0, a1, b0, b1, rho)
}

#' Orthogonal Crawford-Ferguson factor rotation
#'
#' Rotate a loading matrix orthogonally under the Crawford-Ferguson criterion using a
#' gradient-projection optimizer along the orthogonal (Stiefel) manifold.
#'
#' The criterion value `f` and its gradient `dQ/dL` at the rotated loadings
#' `L = A %*% T` define the search; the engine maps the gradient to the orthogonal
#' transformation `T`, projects it onto the tangent space, performs a
#' non-monotone line search, and retracts back onto the orthogonal group via a
#' polar (singular value) projection. `kappa = 0` is the quartimax criterion and
#' `kappa = ncol(A) / (2 * nrow(A))` is the equamax criterion.
#'
#' Additional random orthogonal starts may be requested. To bound runtime the solver
#' screens each random start by its objective, runs a short triage optimization on the
#' best-screened starts, and fully optimizes only those that improve on the current
#' incumbent by at least `triage_improve_tol`.
#'
#' @param L Numeric matrix. The unrotated loading matrix (variables by factors).
#' @param kappa Numeric scalar in `[0, 1]`. The Crawford-Ferguson parameter.
#' @param eps Numeric scalar. Convergence tolerance for the projected-gradient norm.
#' @param normalize Logical scalar. If `TRUE`, apply Kaiser normalization before
#'   rotation and reverse it afterwards.
#' @param random_starts Integer scalar. Number of additional random orthogonal starts.
#' @param maxit Integer scalar. Maximum number of projected-gradient updates.
#' @param max_line_search Integer scalar. Maximum number of step-halving attempts after
#'   the initial trial step in each line-search phase.
#' @param step0 Numeric scalar. Initial step size used in the projected-gradient update.
#' @param screen_keep Integer scalar. Number of screened random starts retained for
#'   triage optimization.
#' @param triage_maxit Integer scalar. Number of short optimization iterations used in
#'   the triage stage.
#' @param triage_improve_tol Numeric scalar. Relative improvement required for a triaged
#'   start to be promoted to full optimization.
#'
#' @returns A named list with the rotated loadings, the orthogonal rotation matrix `Th`
#'   (with `L %*% Th` reproducing the rotated loadings), the attained criterion value, and
#'   the convergence and validity flags. The list additionally reports the
#'   criterion value reached at each optimized start in `all_values`, with a per-start convergence flag in `all_converged`.
#'
#' @references
#' Bernaards, C. A., & Jennrich, R. I. (2005). Gradient projection algorithms and
#' software for arbitrary rotation criteria in factor analysis. *Educational and
#' Psychological Measurement*, 65, 676-696.
#'
#' Crawford, C. B., & Ferguson, G. A. (1970). A general rotation criterion and its use
#' in orthogonal rotation. *Psychometrika*, 35, 321-332.
#'
#' @keywords internal
.rotate_cf_orth <- function(L, kappa, eps = 1e-5, normalize = TRUE, random_starts = 0L, maxit = 1000L, max_line_search = 10L, step0 = 1.0, screen_keep = 5L, triage_maxit = 25L, triage_improve_tol = 0.0) {
    .Call(`_EFAtools_rotate_cf_orth`, L, kappa, eps, normalize, random_starts, maxit, max_line_search, step0, screen_keep, triage_maxit, triage_improve_tol)
}

#' Oblique oblimin factor rotation
#'
#' Rotate a loading matrix obliquely under the oblimin criterion using a
#' gradient-projection optimizer along the oblique (column-normalized) manifold.
#'
#' The criterion value `f` and its gradient `dQ/dL` at the rotated loadings
#' `L = A %*% solve(t(T))` define the search; the engine maps the gradient to the
#' transformation `T` on the manifold `diag(t(T) %*% T) = 1`, projects it onto the tangent
#' space, performs a non-monotone line search, and retracts back onto the manifold
#' by column normalization. `gam = 0` is the quartimin criterion.
#'
#' Additional random starts may be requested. To bound runtime the solver screens each
#' random start by its objective, runs a short triage optimization on the best-screened
#' starts, and fully optimizes only those that improve on the current incumbent by at
#' least `triage_improve_tol`.
#'
#' @inheritParams .rotate_cf_orth
#' @param gam Numeric scalar. The oblimin parameter; `gam = 0` is the quartimin criterion
#'   and the recommended setting. Larger values increasingly reward correlated factors and
#'   can drive the solution toward factor collapse, where the rotation may also fail to
#'   converge; inspect `Phi` before interpreting a fit with `gam > 0`.
#'
#' @returns A named list with the rotated loadings, the transformation matrix `Th`
#'   (with `L %*% t(solve(Th))` reproducing the rotated loadings), the factor correlation
#'   matrix `Phi` (`t(Th) %*% Th`), the attained criterion value, and the convergence and
#'   validity flags. The list additionally reports the criterion value reached
#'   at each optimized start in `all_values`, with a per-start convergence flag in `all_converged`.
#'
#' @references
#' Bernaards, C. A., & Jennrich, R. I. (2005). Gradient projection algorithms and
#' software for arbitrary rotation criteria in factor analysis. *Educational and
#' Psychological Measurement*, 65, 676-696.
#'
#' Jennrich, R. I., & Sampson, P. F. (1966). Rotation for simple loadings.
#' *Psychometrika*, 31, 313-323.
#'
#' @keywords internal
.rotate_oblimin <- function(L, gam = 0.0, eps = 1e-5, normalize = TRUE, random_starts = 0L, maxit = 1000L, max_line_search = 10L, step0 = 1.0, screen_keep = 5L, triage_maxit = 25L, triage_improve_tol = 0.0) {
    .Call(`_EFAtools_rotate_oblimin`, L, gam, eps, normalize, random_starts, maxit, max_line_search, step0, screen_keep, triage_maxit, triage_improve_tol)
}

#' Orthogonal geomin factor rotation
#'
#' Rotate a loading matrix orthogonally under the geomin criterion using a
#' gradient-projection optimizer along the orthogonal (Stiefel) manifold.
#'
#' The criterion value `f` and its gradient `dQ/dL` at the rotated loadings
#' `L = A %*% T` define the search; the engine maps the gradient to the orthogonal
#' transformation `T`, projects it onto the tangent space, performs a non-monotone
#' line search, and retracts back onto the orthogonal group via a polar (singular value)
#' projection. The geomin criterion sums the per-variable geometric mean of the squared
#' loadings offset by `delta`; it is prone to local minima, so additional random starts are
#' recommended.
#'
#' Additional random orthogonal starts may be requested. To bound runtime the solver screens
#' each random start by its objective, runs a short triage optimization on the best-screened
#' starts, and fully optimizes only those that improve on the current incumbent by at least
#' `triage_improve_tol`.
#'
#' @inheritParams .rotate_cf_orth
#' @param delta Numeric scalar. The geomin offset added to the squared loadings; must be a
#'   positive finite scalar. `delta = 0.01` is the usual default.
#'
#' @returns A named list with the rotated loadings, the orthogonal rotation matrix `Th`
#'   (with `L %*% Th` reproducing the rotated loadings), the attained criterion value, and the
#'   convergence and validity flags. The list additionally reports the criterion
#'   value reached at each optimized start in `all_values`, with a per-start convergence flag in `all_converged`.
#'
#' @references
#' Bernaards, C. A., & Jennrich, R. I. (2005). Gradient projection algorithms and
#' software for arbitrary rotation criteria in factor analysis. *Educational and
#' Psychological Measurement*, 65, 676-696.
#'
#' Browne, M. W. (2001). An overview of analytic rotation in exploratory factor analysis.
#' *Multivariate Behavioral Research*, 36, 111-150.
#'
#' @keywords internal
.rotate_geomin_orth <- function(L, delta = 0.01, eps = 1e-5, normalize = TRUE, random_starts = 0L, maxit = 1000L, max_line_search = 10L, step0 = 1.0, screen_keep = 5L, triage_maxit = 25L, triage_improve_tol = 0.0) {
    .Call(`_EFAtools_rotate_geomin_orth`, L, delta, eps, normalize, random_starts, maxit, max_line_search, step0, screen_keep, triage_maxit, triage_improve_tol)
}

#' Oblique geomin factor rotation
#'
#' Rotate a loading matrix obliquely under the geomin criterion using a gradient-projection
#' optimizer along the oblique (column-normalized) manifold.
#'
#' The criterion value `f` and its gradient `dQ/dL` at the rotated loadings
#' `L = A %*% solve(t(T))` define the search; the engine maps the gradient to the
#' transformation `T` on the manifold `diag(t(T) %*% T) = 1`, projects it onto the tangent
#' space, performs a non-monotone line search, and retracts back onto the manifold by
#' column normalization. The geomin criterion sums the per-variable geometric mean of the
#' squared loadings offset by `delta`; it is prone to local minima, so additional random
#' starts are recommended.
#'
#' Additional random starts may be requested. To bound runtime the solver screens each random
#' start by its objective, runs a short triage optimization on the best-screened starts, and
#' fully optimizes only those that improve on the current incumbent by at least
#' `triage_improve_tol`.
#'
#' @inheritParams .rotate_cf_orth
#' @param delta Numeric scalar. The geomin offset added to the squared loadings; must be a
#'   positive finite scalar. `delta = 0.01` is the usual default.
#'
#' @returns A named list with the rotated loadings, the transformation matrix `Th`
#'   (with `L %*% t(solve(Th))` reproducing the rotated loadings), the factor correlation
#'   matrix `Phi` (`t(Th) %*% Th`), the attained criterion value, and the convergence and
#'   validity flags. The list additionally reports the criterion value reached
#'   at each optimized start in `all_values`, with a per-start convergence flag in `all_converged`.
#'
#' @references
#' Bernaards, C. A., & Jennrich, R. I. (2005). Gradient projection algorithms and
#' software for arbitrary rotation criteria in factor analysis. *Educational and
#' Psychological Measurement*, 65, 676-696.
#'
#' Browne, M. W. (2001). An overview of analytic rotation in exploratory factor analysis.
#' *Multivariate Behavioral Research*, 36, 111-150.
#'
#' @keywords internal
.rotate_geomin_oblq <- function(L, delta = 0.01, eps = 1e-5, normalize = TRUE, random_starts = 0L, maxit = 1000L, max_line_search = 10L, step0 = 1.0, screen_keep = 5L, triage_maxit = 25L, triage_improve_tol = 0.0) {
    .Call(`_EFAtools_rotate_geomin_oblq`, L, delta, eps, normalize, random_starts, maxit, max_line_search, step0, screen_keep, triage_maxit, triage_improve_tol)
}

#' Orthogonal Bentler factor rotation
#'
#' Rotate a loading matrix orthogonally under Bentler's invariant pattern simplicity criterion
#' using a gradient-projection optimizer along the orthogonal (Stiefel) manifold.
#'
#' The criterion value `f` and its gradient `dQ/dL` at the rotated loadings
#' `L = A %*% T` define the search; the engine maps the gradient to the orthogonal
#' transformation `T`, projects it onto the tangent space, performs a non-monotone line
#' search, and retracts back onto the orthogonal group via a polar (singular value) projection.
#' The Bentler criterion measures the departure of the cross-products of squared loadings from a
#' diagonal pattern; it is prone to local minima, so additional random starts are recommended.
#'
#' Additional random orthogonal starts may be requested. To bound runtime the solver screens
#' each random start by its objective, runs a short triage optimization on the best-screened
#' starts, and fully optimizes only those that improve on the current incumbent by at least
#' `triage_improve_tol`.
#'
#' @inheritParams .rotate_cf_orth
#'
#' @returns A named list with the rotated loadings, the orthogonal rotation matrix `Th`
#'   (with `L %*% Th` reproducing the rotated loadings), the attained criterion value, and the
#'   convergence and validity flags. The list additionally reports the criterion
#'   value reached at each optimized start in `all_values`, with a per-start convergence flag in `all_converged`.
#'
#' @references
#' Bentler, P. M. (1977). Factor simplicity index and transformations. *Psychometrika*, 42,
#' 277-295.
#'
#' Bernaards, C. A., & Jennrich, R. I. (2005). Gradient projection algorithms and
#' software for arbitrary rotation criteria in factor analysis. *Educational and
#' Psychological Measurement*, 65, 676-696.
#'
#' @keywords internal
.rotate_bentler_orth <- function(L, eps = 1e-5, normalize = TRUE, random_starts = 0L, maxit = 1000L, max_line_search = 10L, step0 = 1.0, screen_keep = 5L, triage_maxit = 25L, triage_improve_tol = 0.0) {
    .Call(`_EFAtools_rotate_bentler_orth`, L, eps, normalize, random_starts, maxit, max_line_search, step0, screen_keep, triage_maxit, triage_improve_tol)
}

#' Oblique Bentler factor rotation
#'
#' Rotate a loading matrix obliquely under Bentler's invariant pattern simplicity criterion
#' using a gradient-projection optimizer along the oblique (column-normalized) manifold.
#'
#' The criterion value `f` and its gradient `dQ/dL` at the rotated loadings
#' `L = A %*% solve(t(T))` define the search; the engine maps the gradient to the
#' transformation `T` on the manifold `diag(t(T) %*% T) = 1`, projects it onto the tangent
#' space, performs a non-monotone line search, and retracts back onto the manifold by
#' column normalization. The Bentler criterion measures the departure of the cross-products of
#' squared loadings from a diagonal pattern; it is prone to local minima, so additional random
#' starts are recommended.
#'
#' Additional random starts may be requested. To bound runtime the solver screens each random
#' start by its objective, runs a short triage optimization on the best-screened starts, and
#' fully optimizes only those that improve on the current incumbent by at least
#' `triage_improve_tol`.
#'
#' @inheritParams .rotate_cf_orth
#'
#' @returns A named list with the rotated loadings, the transformation matrix `Th`
#'   (with `L %*% t(solve(Th))` reproducing the rotated loadings), the factor correlation
#'   matrix `Phi` (`t(Th) %*% Th`), the attained criterion value, and the convergence and
#'   validity flags. The list additionally reports the criterion value reached
#'   at each optimized start in `all_values`, with a per-start convergence flag in `all_converged`.
#'
#' @references
#' Bentler, P. M. (1977). Factor simplicity index and transformations. *Psychometrika*, 42,
#' 277-295.
#'
#' Bernaards, C. A., & Jennrich, R. I. (2005). Gradient projection algorithms and
#' software for arbitrary rotation criteria in factor analysis. *Educational and
#' Psychological Measurement*, 65, 676-696.
#'
#' @keywords internal
.rotate_bentler_oblq <- function(L, eps = 1e-5, normalize = TRUE, random_starts = 0L, maxit = 1000L, max_line_search = 10L, step0 = 1.0, screen_keep = 5L, triage_maxit = 25L, triage_improve_tol = 0.0) {
    .Call(`_EFAtools_rotate_bentler_oblq`, L, eps, normalize, random_starts, maxit, max_line_search, step0, screen_keep, triage_maxit, triage_improve_tol)
}

#' Orthogonal bifactor factor rotation
#'
#' Rotate a loading matrix orthogonally under the Jennrich-Bentler bifactor criterion using a
#' gradient-projection optimizer along the orthogonal (Stiefel) manifold.
#'
#' The criterion value `f` and its gradient `dQ/dL` at the rotated loadings
#' `L = A %*% T` define the search; the engine maps the gradient to the orthogonal
#' transformation `T`, projects it onto the tangent space, performs a non-monotone line
#' search, and retracts back onto the orthogonal group via a polar (singular value) projection.
#' The first factor is treated as a general factor and is exempt from the penalty; the criterion
#' measures the between-group-factor cross-products of the squared loadings, so it is minimized
#' when each variable loads on the general factor plus at most one group factor. The criterion is
#' prone to local minima, so additional random starts are recommended.
#'
#' Additional random orthogonal starts may be requested. To bound runtime the solver screens
#' each random start by its objective, runs a short triage optimization on the best-screened
#' starts, and fully optimizes only those that improve on the current incumbent by at least
#' `triage_improve_tol`.
#'
#' @inheritParams .rotate_cf_orth
#'
#' @returns A named list with the rotated loadings, the orthogonal rotation matrix `Th`
#'   (with `L %*% Th` reproducing the rotated loadings), the attained criterion value, and the
#'   convergence and validity flags. The list additionally reports the criterion
#'   value reached at each optimized start in `all_values`, with a per-start convergence flag in `all_converged`.
#'
#' @references
#' Bernaards, C. A., & Jennrich, R. I. (2005). Gradient projection algorithms and
#' software for arbitrary rotation criteria in factor analysis. *Educational and
#' Psychological Measurement*, 65, 676-696.
#'
#' Jennrich, R. I., & Bentler, P. M. (2011). Exploratory bi-factor analysis. *Psychometrika*,
#' 76, 537-549.
#'
#' @keywords internal
.rotate_bifactor_orth <- function(L, eps = 1e-5, normalize = TRUE, random_starts = 0L, maxit = 1000L, max_line_search = 10L, step0 = 1.0, screen_keep = 5L, triage_maxit = 25L, triage_improve_tol = 0.0) {
    .Call(`_EFAtools_rotate_bifactor_orth`, L, eps, normalize, random_starts, maxit, max_line_search, step0, screen_keep, triage_maxit, triage_improve_tol)
}

#' Oblique bifactor factor rotation
#'
#' Rotate a loading matrix obliquely under the Jennrich-Bentler bifactor criterion using a
#' gradient-projection optimizer along the oblique (column-normalized) manifold.
#'
#' The criterion value `f` and its gradient `dQ/dL` at the rotated loadings
#' `L = A %*% solve(t(T))` define the search; the engine maps the gradient to the
#' transformation `T` on the manifold `diag(t(T) %*% T) = 1`, projects it onto the tangent
#' space, performs a non-monotone line search, and retracts back onto the manifold by
#' column normalization. The first factor is treated as a general factor and is exempt from the
#' penalty; the criterion measures the between-group-factor cross-products of the squared
#' loadings, so it is minimized when each variable loads on the general factor plus at most one
#' group factor. The criterion is prone to local minima, so additional random starts are
#' recommended.
#'
#' Additional random starts may be requested. To bound runtime the solver screens each random
#' start by its objective, runs a short triage optimization on the best-screened starts, and
#' fully optimizes only those that improve on the current incumbent by at least
#' `triage_improve_tol`.
#'
#' @inheritParams .rotate_cf_orth
#'
#' @returns A named list with the rotated loadings, the transformation matrix `Th`
#'   (with `L %*% t(solve(Th))` reproducing the rotated loadings), the factor correlation
#'   matrix `Phi` (`t(Th) %*% Th`), the attained criterion value, and the convergence and
#'   validity flags. The list additionally reports the criterion value reached
#'   at each optimized start in `all_values`, with a per-start convergence flag in `all_converged`.
#'
#' @references
#' Bernaards, C. A., & Jennrich, R. I. (2005). Gradient projection algorithms and
#' software for arbitrary rotation criteria in factor analysis. *Educational and
#' Psychological Measurement*, 65, 676-696.
#'
#' Jennrich, R. I., & Bentler, P. M. (2011). Exploratory bi-factor analysis. *Psychometrika*,
#' 76, 537-549.
#'
#' @keywords internal
.rotate_bifactor_oblq <- function(L, eps = 1e-5, normalize = TRUE, random_starts = 0L, maxit = 1000L, max_line_search = 10L, step0 = 1.0, screen_keep = 5L, triage_maxit = 25L, triage_improve_tol = 0.0) {
    .Call(`_EFAtools_rotate_bifactor_oblq`, L, eps, normalize, random_starts, maxit, max_line_search, step0, screen_keep, triage_maxit, triage_improve_tol)
}

#' Oblique simplimax factor rotation
#'
#' Rotate a loading matrix obliquely under the simplimax criterion using a gradient-projection
#' optimizer along the oblique (column-normalized) manifold.
#'
#' The criterion value `f` and its gradient `dQ/dL` at the rotated loadings
#' `L = A %*% solve(t(T))` define the search; the engine maps the gradient to the
#' transformation `T` on the manifold `diag(t(T) %*% T) = 1`, projects it onto the tangent
#' space, performs a non-monotone line search, and retracts back onto the manifold by column
#' normalization. The simplimax criterion sums the `k` smallest squared loadings, so it is
#' minimized when the `k` "close-to-zero" loadings are driven toward zero; the count `k` is a
#' tuning parameter. Because the set of `k` smallest loadings is reselected at every evaluation,
#' the criterion is only piecewise smooth: its gradient jumps as loadings cross the kth-smallest threshold, so the
#' line search accepts a step whenever it decreases the largest objective over a short window of
#' recent iterations (a non-monotone test; Grippo, Lampariello, & Lucidi, 1986), letting the
#' optimizer step across the kinks where a strictly monotone descent would stall.
#'
#' The criterion is strongly prone to local minima, so the solver fully optimizes the identity
#' start together with `random_starts` random orthogonal starts and keeps the solution with the
#' lowest criterion value. Fully optimizing every start -- rather than the screen-and-triage
#' strategy used for the smooth criteria, which assumes the rational start lies in the global
#' basin -- is the standard remedy for the local minima of complexity-based rotation criteria
#' (Kiers, 1994; Browne, 2001).
#'
#' @inheritParams .rotate_cf_orth
#' @param k Integer scalar. The number of "close-to-zero" loadings the criterion targets; must
#'   be in `[1, nrow(L) * ncol(L)]`. `k = nrow(L)` is the usual default.
#' @param eps Numeric scalar. Convergence tolerance for the projected-gradient norm. Because the
#'   simplimax criterion is only piecewise smooth, the projected gradient need not reach this
#'   tolerance at the optimum; convergence is then reported when the criterion value stalls (the
#'   non-monotone search described above), so `eps` mainly governs the smooth phases of the search.
#' @param random_starts Integer scalar. Number of random orthogonal starts fully optimized in
#'   addition to the identity start.
#' @param maxit Integer scalar. Maximum number of projected-gradient updates per start.
#'
#' @returns A named list with the rotated loadings, the transformation matrix `Th`
#'   (with `L %*% t(solve(Th))` reproducing the rotated loadings), the factor correlation
#'   matrix `Phi` (`t(Th) %*% Th`), the attained criterion value, and the convergence and
#'   validity flags. The list additionally reports the criterion value reached
#'   at each optimized start in `all_values`, with a per-start convergence flag in `all_converged`.
#'
#' @references
#' Bernaards, C. A., & Jennrich, R. I. (2005). Gradient projection algorithms and
#' software for arbitrary rotation criteria in factor analysis. *Educational and
#' Psychological Measurement*, 65, 676-696.
#'
#' Browne, M. W. (2001). An overview of analytic rotation in exploratory factor analysis.
#' *Multivariate Behavioral Research*, 36, 111-150.
#'
#' Grippo, L., Lampariello, F., & Lucidi, S. (1986). A nonmonotone line search technique for
#' Newton's method. *SIAM Journal on Numerical Analysis*, 23, 707-716.
#'
#' Kiers, H. A. L. (1994). Simplimax: Oblique rotation to an optimal target with simple
#' structure. *Psychometrika*, 59, 567-579.
#'
#' @keywords internal
.rotate_simplimax_oblq <- function(L, k, eps = 1e-5, normalize = TRUE, random_starts = 0L, maxit = 1000L, max_line_search = 10L, step0 = 1.0) {
    .Call(`_EFAtools_rotate_simplimax_oblq`, L, k, eps, normalize, random_starts, maxit, max_line_search, step0)
}

#' Rotation Jacobians for analytic rotation standard errors
#'
#' Forward-difference the warm-started re-rotation map `A -> (rotated loadings, Phi)` over the
#' unrotated loadings `A` to obtain the rotation Jacobians used by the analytic standard errors for
#' rotated loadings (`se = "information"` in [efa_fit()]). The full `nrow(A) * ncol(A)` finite-
#' difference loop runs in compiled code, re-solving the rotation from the converged transformation
#' `T_init` at each perturbation; the caller forms `J V J'` in R.
#'
#' @param A Numeric matrix. The unrotated loading matrix at the solution.
#' @param T_init Numeric matrix. The converged transformation that warm-starts each re-rotation.
#' @param method Character scalar. The criterion family: one of `"cf"`, `"oblimin"`, `"geomin"`,
#'   `"bentler"`, `"bifactor"`.
#' @param param Numeric scalar. The criterion's tuning argument (`kappa` for `"cf"`, `gam` for
#'   `"oblimin"`, `delta` for `"geomin"`); ignored for `"bentler"` and `"bifactor"`.
#' @param normalize Logical scalar. Apply Kaiser normalization before rotation and reverse it after.
#' @param oblique Logical scalar. Use the oblique (column-normalized) manifold; otherwise orthogonal.
#' @param eps Numeric scalar. The forward-difference step on the loadings.
#' @param general_col Integer scalar. For `method = "bifactor"`, the zero-based column holding the
#'   general factor in the (factor-reordered) reported solution; ignored by the other criteria.
#'
#' @returns A named list with the Jacobian `J_L` (`pk x pk`), the re-rotated `base_loadings`, a
#'   validity flag, and -- when oblique -- the Jacobian `J_Phi` (`k^2 x pk`) and `base_Phi`.
#'
#' @references
#' Jennrich, R. I. (1973). Standard errors for obliquely rotated factor loadings.
#' *Psychometrika*, 38, 593-604.
#'
#' @keywords internal
.rotation_se_jacobian <- function(A, T_init, method, param, normalize, oblique, eps, general_col = 0L) {
    .Call(`_EFAtools_rotation_se_jacobian`, A, T_init, method, param, normalize, oblique, eps, general_col)
}

#' Draw multivariate-normal data from a population correlation matrix.
#'
#' Internal helper called from efa_simulate(). Draws `N` cases from a
#' `p`-variate normal with correlation (or covariance) `R` by post-multiplying a
#' matrix of standard normal deviates by a matrix square root `M` of `R` (with
#' M' M = R, so the rows of Z * M are N(0, R)). This is the same Z * M rule used by
#' the NEST reference simulation (.simulate_cfm_eigen): there `M` is the transposed
#' factor-score matrix, here it is a Cholesky or eigen square root.
#' A positive-definite `R` is factored by Cholesky; a positive-semidefinite but
#' singular `R` (which makes the Cholesky fail although it is still a valid
#' covariance, e.g. a no-factor block or a smoothed factor intercorrelation
#' matrix) falls back to a symmetric eigen square root.
#'
#' @param R numeric matrix. Population correlation/covariance matrix.
#' @param N integer. Number of cases to draw.
#' @param tol numeric. Eigenvalues below `-tol` mark `R` as indefinite.
#'
#' @keywords internal
.simulate_cfm_mvn <- function(R, N, tol = 1e-8) {
    .Call(`_EFAtools_simulate_cfm_mvn`, R, N, tol)
}

#' Reference eigenvalues for the efa_nest() simulation via the shared kernel.
#'
#' Internal helper called from efa_nest(). Simulates `nreps` datasets from an
#' `(nf - 1)`-factor reference model, given that model's loadings `Lambda` and
#' uniquenesses `Psi`, and returns the `nf`-th largest eigenvalue of each simulated
#' correlation matrix. The data are drawn with the shared Z * M rule (see
#' .simulate_cfm_mvn()) using the factor-score square root
#' `M = t([Lambda | diag(sqrt(Psi))])`, so a row `randn(1, nf - 1 + p) * M` is
#' N(0, Lambda Lambda' + diag(Psi)). Drawing `nf - 1 + p` standard normals and
#' post-multiplying by the factor-score matrix is faster than forming the model-
#' implied correlation matrix and drawing from it, and matches the position at which
#' efa_nest() reads the reference eigenvalue.
#'
#' @param nf integer. Position of the empirical eigenvalue being tested (1-based);
#'   the `nf`-th largest simulated eigenvalue is returned per replicate.
#' @param N integer. Number of cases / observations per simulated dataset.
#' @param Lambda numeric matrix. Loadings of the `(nf - 1)`-factor reference model
#'   (`p x (nf - 1)`); pass a `p x 0` matrix for the `nf == 1` null (identity) model.
#' @param Psi numeric vector. Uniquenesses (`1 - h2`) of the reference model.
#' @param nreps integer. Number of datasets to simulate.
#'
#' @keywords internal
.simulate_cfm_eigen <- function(nf, N, Lambda, Psi, nreps = 1000L) {
    .Call(`_EFAtools_simulate_cfm_eigen`, nf, N, Lambda, Psi, nreps)
}