Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

adapt for matrix operations via lavaanC #402

Merged
merged 1 commit into from
Dec 18, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
119 changes: 111 additions & 8 deletions R/lav_matrix.R
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,9 @@ lav_matrix_vecr <- function(A) {
# faster way??
# nRow <- NROW(A); nCol <- NCOL(A)
# idx <- (seq_len(nCol) - 1L) * nRow + rep(seq_len(nRow), each = nCol)

if (lav_use_lavaanC() && is.numeric(A)) {
return(lavaanC::m_vecr(A))
}
lav_matrix_vec(t(A))
}

Expand All @@ -39,6 +41,9 @@ lav_matrix_vecr <- function(A) {
# M&N book: page 48-49
#
lav_matrix_vech <- function(S, diagonal = TRUE) {
if (lav_use_lavaanC() && is.numeric(S)) {
return(lavaanC::m_vech(S, diagonal))
}
ROW <- row(S)
COL <- col(S)
if (diagonal) S[ROW >= COL] else S[ROW > COL]
Expand All @@ -49,6 +54,9 @@ lav_matrix_vech <- function(S, diagonal = TRUE) {
# into a vector by stacking the *rows* of the matrix one after the
# other, but eliminating all supradiagonal elements
lav_matrix_vechr <- function(S, diagonal = TRUE) {
if (lav_use_lavaanC() && is.numeric(S)) {
return(lavaanC::m_vechr(S, diagonal))
}
S[lav_matrix_vechr_idx(n = NCOL(S), diagonal = diagonal)]
}

Expand All @@ -57,26 +65,31 @@ lav_matrix_vechr <- function(S, diagonal = TRUE) {
# into a vector by stacking the *columns* of the matrix one after the
# other, but eliminating all infradiagonal elements
lav_matrix_vechu <- function(S, diagonal = TRUE) {
if (lav_use_lavaanC() && is.numeric(S)) {
return(lavaanC::m_vechu(S, diagonal))
}
S[lav_matrix_vechu_idx(n = NCOL(S), diagonal = diagonal)]
}


# the vechru operator transforms a *symmetric* matrix
# into a vector by stacking the *rows* of the matrix one after the
# other, but eliminating all infradiagonal elements
#
# same as vech (but using upper-diagonal elements)
lav_matrix_vechru <- function(S, diagonal = TRUE) {
if (lav_use_lavaanC() && is.numeric(S)) {
return(lavaanC::m_vechru(S, diagonal))
}
S[lav_matrix_vechru_idx(n = NCOL(S), diagonal = diagonal)]
}




# return the *vector* indices of the lower triangular elements of a
# symmetric matrix of size 'n'
lav_matrix_vech_idx <- function(n = 1L, diagonal = TRUE) {
n <- as.integer(n)
if (lav_use_lavaanC() && n > 1L) {
return(lavaanC::m_vech_idx(n, diagonal))
}
if (n < 100L) {
ROW <- matrix(seq_len(n), n, n)
COL <- matrix(seq_len(n), n, n, byrow = TRUE)
Expand All @@ -101,6 +114,9 @@ lav_matrix_vech_idx <- function(n = 1L, diagonal = TRUE) {
# symmetric matrix of size 'n'
lav_matrix_vech_row_idx <- function(n = 1L, diagonal = TRUE) {
n <- as.integer(n)
if (lav_use_lavaanC() && n > 1L) {
return(lavaanC::m_vech_row_idx(n, diagonal))
}
if (diagonal) {
unlist(lapply(seq_len(n), seq.int, n))
} else {
Expand All @@ -112,6 +128,9 @@ lav_matrix_vech_row_idx <- function(n = 1L, diagonal = TRUE) {
# symmetric matrix of size 'n'
lav_matrix_vech_col_idx <- function(n = 1L, diagonal = TRUE) {
n <- as.integer(n)
if (lav_use_lavaanC() && n > 1L) {
return(lavaanC::m_vech_col_idx(n, diagonal))
}
if (!diagonal) {
n <- n - 1L
}
Expand All @@ -125,6 +144,9 @@ lav_matrix_vech_col_idx <- function(n = 1L, diagonal = TRUE) {
# symmetric matrix of size 'n' -- ROW-WISE
lav_matrix_vechr_idx <- function(n = 1L, diagonal = TRUE) {
n <- as.integer(n)
if (lav_use_lavaanC() && n > 1L) {
return(lavaanC::m_vechr_idx(n, diagonal))
}
if (n < 100L) {
ROW <- matrix(seq_len(n), n, n)
COL <- matrix(seq_len(n), n, n, byrow = TRUE)
Expand All @@ -149,6 +171,9 @@ lav_matrix_vechr_idx <- function(n = 1L, diagonal = TRUE) {
# symmetric matrix of size 'n' -- COLUMN-WISE
lav_matrix_vechu_idx <- function(n = 1L, diagonal = TRUE) {
n <- as.integer(n)
if (lav_use_lavaanC() && n > 1L) {
return(lavaanC::m_vechu_idx(n, diagonal))
}
if (n < 100L) {
ROW <- matrix(seq_len(n), n, n)
COL <- matrix(seq_len(n), n, n, byrow = TRUE)
Expand All @@ -166,6 +191,9 @@ lav_matrix_vechu_idx <- function(n = 1L, diagonal = TRUE) {
# symmetric matrix of size 'n' -- ROW-WISE
lav_matrix_vechru_idx <- function(n = 1L, diagonal = TRUE) {
n <- as.integer(n)
if (lav_use_lavaanC() && n > 1L) {
return(lavaanC::m_vechru_idx(n, diagonal))
}
if (n < 100L) {
# FIXME!! make this more efficient (without creating 3 n*n matrices!)
ROW <- matrix(seq_len(n), n, n)
Expand Down Expand Up @@ -196,6 +224,9 @@ lav_matrix_vechru_idx <- function(n = 1L, diagonal = TRUE) {
lav_matrix_vech_reverse <- lav_matrix_vechru_reverse <-
lav_matrix_upper2full <-
function(x, diagonal = TRUE) {
if (lav_use_lavaanC()) {
return(lavaanC::m_vech_reverse(x, diagonal))
}
# guess dimensions
if (diagonal) {
p <- (sqrt(1 + 8 * length(x)) - 1) / 2
Expand All @@ -218,6 +249,9 @@ lav_matrix_vech_reverse <- lav_matrix_vechru_reverse <-
# reconstruct S
lav_matrix_vechr_reverse <- lav_matrix_vechu_reverse <-
lav_matrix_lower2full <- function(x, diagonal = TRUE) {
if (lav_use_lavaanC()) {
return(lavaanC::m_vechr_reverse(x, diagonal))
}
# guess dimensions
if (diagonal) {
p <- (sqrt(1 + 8 * length(x)) - 1) / 2
Expand All @@ -239,19 +273,27 @@ lav_matrix_vechr_reverse <- lav_matrix_vechu_reverse <-
# matrix of size 'n'
lav_matrix_diag_idx <- function(n = 1L) {
# if(n < 1L) return(integer(0L))
n <- as.integer(n)
if (lav_use_lavaanC()) {
return(lavaanC::m_diag_idx(n))
}
1L + (seq_len(n) - 1L) * (n + 1L)
}


# return the *vector* indices of the diagonal elements of the LOWER part
# of a symmatrix matrix of size 'n'
lav_matrix_diagh_idx <- function(n = 1L) {
n <- as.integer(n)
if (n < 1L) {
return(integer(0L))
}
if (n == 1L) {
return(1L)
}
if (lav_use_lavaanC()) {
return(lavaanC::m_diagh_idx(n))
}
c(1L, cumsum(n:2L) + 1L)
}

Expand All @@ -262,6 +304,9 @@ lav_matrix_antidiag_idx <- function(n = 1L) {
if (n < 1L) {
return(integer(0L))
}
if (lav_use_lavaanC()) {
return(lavaanC::m_antidiag_idx(n))
}
1L + seq_len(n) * (n - 1L)
}

Expand Down Expand Up @@ -294,6 +339,10 @@ lav_matrix_vech_which_idx <- function(n = 1L, diagonal = TRUE,
return(integer(0L))
}
n <- as.integer(n)
idx <- as.integer(idx)
if (lav_use_lavaanC()) {
return(lavaanC::m_vech_which_idx(n, diagonal, idx, type, add.idx.at.start))
}
A <- matrix(FALSE, n, n)
if (type == "and") {
A[idx, idx] <- TRUE
Expand All @@ -319,6 +368,10 @@ lav_matrix_vech_match_idx <- function(n = 1L, diagonal = TRUE,
return(integer(0L))
}
n <- as.integer(n)
idx <- as.integer(idx)
if (lav_use_lavaanC()) {
return(lavaanC::m_vech_match_idx(n, diagonal, idx))
}
pstar <- n * (n + 1) / 2
A <- lav_matrix_vech_reverse(seq_len(pstar))
B <- A[idx, idx, drop = FALSE]
Expand All @@ -329,7 +382,9 @@ lav_matrix_vech_match_idx <- function(n = 1L, diagonal = TRUE,
lav_matrix_is_diagonal <- function(A = NULL) {
A <- as.matrix.default(A)
stopifnot(nrow(A) == ncol(A))

if (lav_use_lavaanC()) {
return(lavaanC::m_is_diagonal(A))
}
diag(A) <- 0
all(A == 0)
}
Expand Down Expand Up @@ -413,6 +468,10 @@ lav_matrix_is_diagonal <- function(A = NULL) {
# dup3: using col idx only
# D7 <- dup(7L); x<- apply(D7, 1, function(x) which(x > 0)); matrix(x,7,7)
.dup3 <- function(n = 1L) {
n <- as.integer(n)
if (lav_use_lavaanC()) {
return(lavaanC::m_duplication(n))
}
if ((n < 1L) || (round(n) != n)) {
lav_msg_stop(gettext("n must be a positive integer"))
}
Expand Down Expand Up @@ -468,6 +527,10 @@ lav_matrix_duplication <- .dup3
# - it returns a matrix of size p^2 * (p*(p-1))/2
# - the columns corresponding to the diagonal elements have been removed
lav_matrix_duplication_cor <- function(n = 1L) {
n <- as.integer(n)
if (lav_use_lavaanC()) {
return(lavaanC::m_duplication_cor(n))
}
out <- lav_matrix_duplication(n = n)
diag.idx <- lav_matrix_diagh_idx(n = n)
out[, -diag.idx, drop = FALSE]
Expand All @@ -477,6 +540,9 @@ lav_matrix_duplication_cor <- function(n = 1L) {
# sqrt(nrow(A)) is an integer
# A is not symmetric, and not even square, only n^2 ROWS
lav_matrix_duplication_pre <- function(A = matrix(0, 0, 0)) {
if (lav_use_lavaanC()) {
return(lavaanC::m_duplication_pre(A))
}
# number of rows
n2 <- NROW(A)

Expand Down Expand Up @@ -524,6 +590,9 @@ lav_matrix_duplication_dup_pre2 <- function(A = matrix(0, 0, 0)) {
# A is not symmetric, and not even square, only n^2 ROWS
# correlation version: ignoring diagonal elements
lav_matrix_duplication_cor_pre <- function(A = matrix(0, 0, 0)) {
if (lav_use_lavaanC()) {
return(lavaanC::m_duplication_cor_pre(A))
}
# number of rows
n2 <- NROW(A)

Expand All @@ -548,6 +617,9 @@ lav_matrix_duplication_cor_pre <- function(A = matrix(0, 0, 0)) {
# sqrt(ncol(A)) must be an integer
# A is not symmetric, and not even square, only n^2 COLUMNS
lav_matrix_duplication_post <- function(A = matrix(0, 0, 0)) {
if (lav_use_lavaanC()) {
return(lavaanC::m_duplication_post(A))
}
# number of columns
n2 <- NCOL(A)

Expand All @@ -573,6 +645,9 @@ lav_matrix_duplication_post <- function(A = matrix(0, 0, 0)) {
# A is not symmetric, and not even square, only n^2 COLUMNS
# correlation version: ignoring the diagonal elements
lav_matrix_duplication_cor_post <- function(A = matrix(0, 0, 0)) {
if (lav_use_lavaanC()) {
return(lavaanC::m_duplication_cor_post(A))
}
# number of columns
n2 <- NCOL(A)

Expand All @@ -597,6 +672,9 @@ lav_matrix_duplication_cor_post <- function(A = matrix(0, 0, 0)) {
# compute t(D) %*% A %*% D (without explicitly computing D)
# A must be a square matrix and sqrt(ncol) an integer
lav_matrix_duplication_pre_post <- function(A = matrix(0, 0, 0)) {
if (lav_use_lavaanC()) {
return(lavaanC::m_duplication_pre_post(A))
}
# number of columns
n2 <- NCOL(A)

Expand All @@ -623,6 +701,9 @@ lav_matrix_duplication_pre_post <- function(A = matrix(0, 0, 0)) {
# A must be a square matrix and sqrt(ncol) an integer
# correlation version: ignoring diagonal elements
lav_matrix_duplication_cor_pre_post <- function(A = matrix(0, 0, 0)) {
if (lav_use_lavaanC()) {
return(lavaanC::m_duplication_cor_pre_post(A))
}
# number of columns
n2 <- NCOL(A)

Expand Down Expand Up @@ -772,6 +853,10 @@ lav_matrix_elimination_pre_post <- function(A = matrix(0, 0, 0)) {

# create DUP.ginv without transpose
.dup_ginv2 <- function(n = 1L) {
if (lav_use_lavaanC()) {
n <- as.integer(n)
return(lavaanC::m_duplication_ginv(n))
}
if ((n < 1L) || (round(n) != n)) {
lav_msg_stop(gettext("n must be a positive integer"))
}
Expand Down Expand Up @@ -840,7 +925,9 @@ lav_matrix_duplication_ginv_post <- function(A = matrix(0, 0, 0)) {
# for square matrices only, with ncol = nrow = n^2
lav_matrix_duplication_ginv_pre_post <- function(A = matrix(0, 0, 0)) {
A <- as.matrix.default(A)

if (lav_use_lavaanC()) {
return(lavaanC::m_duplication_ginv_pre_post(A))
}
# number of columns
n2 <- NCOL(A)

Expand Down Expand Up @@ -902,6 +989,12 @@ lav_matrix_duplication_ginv_cor_pre_post <- function(A = matrix(0, 0, 0)) {

# first attempt
.com1 <- function(m = 1L, n = 1L) {
if (lav_use_lavaanC()) {
m <- as.integer(m)
n <- as.integer(n)
return(lavaanC::m_commutation(m, n))
}

if ((m < 1L) || (round(m) != m)) {
lav_msg_stop(gettext("n must be a positive integer"))
}
Expand Down Expand Up @@ -929,6 +1022,10 @@ lav_matrix_commutation <- .com1
# = permuting the rows of A
lav_matrix_commutation_pre <- function(A = matrix(0, 0, 0)) {
A <- as.matrix(A)

if (lav_use_lavaanC()) {
return(lavaanC::m_commutation_pre(A))
}

# number of rows of A
n2 <- nrow(A)
Expand All @@ -952,6 +1049,10 @@ lav_matrix_commutation_pre <- function(A = matrix(0, 0, 0)) {
# = permuting the columns of A
lav_matrix_commutation_post <- function(A = matrix(0, 0, 0)) {
A <- as.matrix(A)

if (lav_use_lavaanC()) {
return(lavaanC::m_commutation_post(A))
}

# number of columns of A
n2 <- ncol(A)
Expand All @@ -975,7 +1076,9 @@ lav_matrix_commutation_post <- function(A = matrix(0, 0, 0)) {
# = permuting both the rows AND columns of A
lav_matrix_commutation_pre_post <- function(A = matrix(0, 0, 0)) {
A <- as.matrix(A)

if (lav_use_lavaanC()) {
return(lavaanC::m_commutation_pre_post(A))
}
# number of columns of A
n2 <- NCOL(A)

Expand Down
11 changes: 9 additions & 2 deletions R/lav_model_information.R
Original file line number Diff line number Diff line change
Expand Up @@ -153,15 +153,22 @@ lav_model_information_expected <- function(lavmodel = NULL,
)
Info.group[[g]] <- fg * Info.g
} else {
# ldw_trace(paste(sum(Delta[[g]] == 0),"/",length(Delta[[g]])))
# compute information for this group
if (lavmodel@estimator %in% c("DWLS", "ULS")) {
# diagonal weight matrix
Delta2 <- sqrt(A1[[g]]) * Delta[[g]]
Info.group[[g]] <- fg * crossprod(Delta2)
} else {
# full weight matrix
Info.group[[g]] <-
fg * (crossprod(Delta[[g]], A1[[g]]) %*% Delta[[g]])
if (lav_use_lavaanC()) {
Info.group[[g]] <-
fg * lavaanC::m_prod(
lavaanC::m_crossprod(Delta[[g]], A1[[g]], "L"), Delta[[g]], "R")
} else {
Info.group[[g]] <-
fg * (crossprod(Delta[[g]], A1[[g]]) %*% Delta[[g]])
}
}
}
} # g
Expand Down
2 changes: 1 addition & 1 deletion R/lav_syntax_parser_cr.R
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ ldw_parse_model_string_cr <- function(model.syntax = "",
)

if (lav_use_lavaanC()) {
flat <- lavaanC::lav_parse_model_string_c(modelsrc)
flat <- lavaanC::parse_model_string(modelsrc)
} else {
flat <- lav_parse_model_string_r(modelsrc)
}
Expand Down
Loading