Skip to content

Commit

Permalink
Merge pull request #406 from lucdw/master
Browse files Browse the repository at this point in the history
adapt function for use m_kronecker_* in lavaanC  if available
  • Loading branch information
yrosseel authored Dec 23, 2024
2 parents e65e642 + 4af12f5 commit 7fcecf6
Show file tree
Hide file tree
Showing 7 changed files with 150 additions and 41 deletions.
35 changes: 28 additions & 7 deletions R/lav_mvnorm.R
Original file line number Diff line number Diff line change
Expand Up @@ -743,11 +743,19 @@ lav_mvnorm_information_expected <- function(Y = NULL, # unused!
)
}

if(correlation) {
I22 <- 0.5 * lav_matrix_duplication_cor_pre_post(Sigma.inv %x%
Sigma.inv)
if (lav_use_lavaanC()) {
if (correlation) {
I22 <- lavaanC::m_kronecker_dup_cor_pre_post(Sigma.inv,
multiplicator = 0.5)
} else {
I22 <- lavaanC::m_kronecker_dup_pre_post(Sigma.inv, multiplicator = 0.5)
}
} else {
if (correlation) {
I22 <- 0.5 * lav_matrix_duplication_cor_pre_post(Sigma.inv %x% Sigma.inv)
} else {
I22 <- 0.5 * lav_matrix_duplication_pre_post(Sigma.inv %x% Sigma.inv)
}
}

# fixed.x?
Expand Down Expand Up @@ -835,7 +843,11 @@ lav_mvnorm_information_observed_samplestats <- function(
}

AAA <- Sigma.inv %*% (2 * W.tilde - Sigma) %*% Sigma.inv
I22 <- (1 / 2) * lav_matrix_duplication_pre_post(Sigma.inv %x% AAA)
if (lav_use_lavaanC()) {
I22 <- lavaanC::m_kronecker_dup_pre_post(Sigma.inv, AAA, 0.5)
} else {
I22 <- (1 / 2) * lav_matrix_duplication_pre_post(Sigma.inv %x% AAA)
}

if (meanstructure) {
out <- rbind(
Expand Down Expand Up @@ -937,8 +949,13 @@ lav_mvnorm_inverted_information_expected <- function(Y = NULL, # unused!
# take difference
R <- Sigma - YbarX.aug

SS <- 2 * lav_matrix_duplication_ginv_pre_post(Sigma %x% Sigma)
RR <- 2 * lav_matrix_duplication_ginv_pre_post(R %x% R)
if (lav_use_lavaanC()) {
SS <- lavaanC::m_kronecker_dup_ginv_pre_post(Sigma, multiplicator = 2.0)
RR <- lavaanC::m_kronecker_dup_ginv_pre_post(R, multiplicator = 2.0)
} else {
SS <- 2 * lav_matrix_duplication_ginv_pre_post(Sigma %x% Sigma)
RR <- 2 * lav_matrix_duplication_ginv_pre_post(R %x% R)
}
I22 <- SS - RR

if (meanstructure) {
Expand All @@ -948,7 +965,11 @@ lav_mvnorm_inverted_information_expected <- function(Y = NULL, # unused!
out <- I22
}
} else {
I22 <- 2 * lav_matrix_duplication_ginv_pre_post(Sigma %x% Sigma)
if (lav_use_lavaanC()) {
I22 <- lavaanC::m_kronecker_dup_ginv_pre_post(Sigma, multiplicator = 2.0)
} else {
I22 <- 2 * lav_matrix_duplication_ginv_pre_post(Sigma %x% Sigma)
}
if (meanstructure) {
I11 <- Sigma
out <- lav_matrix_bdiag(I11, I22)
Expand Down
24 changes: 20 additions & 4 deletions R/lav_mvnorm_cluster.R
Original file line number Diff line number Diff line change
Expand Up @@ -1012,7 +1012,11 @@ lav_mvnorm_cluster_information_expected <- function(Lp = NULL,
omega.j.inv <- solve(omega.j)

I11.j <- omega.j.inv
I22.j <- 0.5 * lav_matrix_duplication_pre_post(omega.j.inv %x% omega.j.inv)
if (lav_use_lavaanC()) {
I22.j <- lavaanC::m_kronecker_dup_pre_post(omega.j.inv, multiplicator = 0.5)
} else {
I22.j <- 0.5 * lav_matrix_duplication_pre_post(omega.j.inv %x% omega.j.inv)
}
I.j <- lav_matrix_bdiag(I11.j, I22.j)
info.j <- t(Delta.j) %*% I.j %*% Delta.j

Expand All @@ -1028,7 +1032,11 @@ lav_mvnorm_cluster_information_expected <- function(Lp = NULL,
Sigma.W.inv.tilde[ov.idx[[1]], ov.idx[[1]]] <- Sigma.W.inv

I11.w <- Sigma.W.inv.tilde
I22.w <- 0.5 * lav_matrix_duplication_pre_post(Sigma.W.inv.tilde %x% Sigma.W.inv.tilde)
if (lav_use_lavaanC()) {
I22.W <- lavaanC::m_kronecker_dup_pre_post(Sigma.W.inv.tilde, multiplicator = 0.5)
} else {
I22.w <- 0.5 * lav_matrix_duplication_pre_post(Sigma.W.inv.tilde %x% Sigma.W.inv.tilde)
}
I.w <- lav_matrix_bdiag(I11.w, I22.w)
information.w <- (nobs - nclusters) *
(t(Delta.W.tilde) %*% I.w %*% Delta.W.tilde)
Expand Down Expand Up @@ -1160,7 +1168,11 @@ lav_mvnorm_cluster_information_expected_delta <- function(Lp = NULL,
omega.j.inv <- solve(omega.j)

I11.j <- omega.j.inv
I22.j <- 0.5 * lav_matrix_duplication_pre_post(omega.j.inv %x% omega.j.inv)
if (lav_use_lavaanC()) {
I22.j <- lavaanC::m_kronecker_dup_pre_post(omega.j.inv, multiplicator = 0.5)
} else {
I22.j <- 0.5 * lav_matrix_duplication_pre_post(omega.j.inv %x% omega.j.inv)
}
I.j <- lav_matrix_bdiag(I11.j, I22.j)
info.j <- t(Delta.j) %*% I.j %*% Delta.j

Expand All @@ -1173,7 +1185,11 @@ lav_mvnorm_cluster_information_expected_delta <- function(Lp = NULL,
Sinv.method = Sinv.method
)
I11.w <- Sigma.W.inv
I22.w <- 0.5 * lav_matrix_duplication_pre_post(Sigma.W.inv %x% Sigma.W.inv)
if (lav_use_lavaanC()) {
I22.w <- lavaanC::m_kronecker_dup_pre_post(Sigma.W.inv, multiplicator = 0.5)
} else {
I22.w <- 0.5 * lav_matrix_duplication_pre_post(Sigma.W.inv %x% Sigma.W.inv)
}
I.w <- lav_matrix_bdiag(I11.w, I22.w)

# force zero for means both.idx in within part
Expand Down
31 changes: 25 additions & 6 deletions R/lav_mvnorm_h1.R
Original file line number Diff line number Diff line change
Expand Up @@ -262,11 +262,21 @@ lav_mvnorm_h1_information_expected <- function(

I11 <- sample.cov.inv
if(correlation) {
I22 <- 0.5 * lav_matrix_duplication_cor_pre_post(sample.cov.inv %x%
if (lav_use_lavaanC()) {
I22 <- lavaanC::m_kronecker_dup_cor_pre_post(sample.cov.inv,
multiplicator = 0.5)
} else {
I22 <- 0.5 * lav_matrix_duplication_cor_pre_post(sample.cov.inv %x%
sample.cov.inv)
}
} else {
I22 <- 0.5 * lav_matrix_duplication_pre_post(sample.cov.inv %x%
if (lav_use_lavaanC()) {
I22 <- lavaanC::m_kronecker_dup_pre_post(sample.cov.inv,
multiplicator = 0.5)
} else {
I22 <- 0.5 * lav_matrix_duplication_pre_post(sample.cov.inv %x%
sample.cov.inv)
}
}

# fixed.x?
Expand Down Expand Up @@ -333,9 +343,14 @@ lav_mvnorm_h1_information_observed_samplestats <- function(
I11[, x.idx] <- 0
}

I22 <- 0.5 * lav_matrix_duplication_pre_post(sample.cov.inv %x%
if (lav_use_lavaanC()) {
I22 <- lavaanC::m_kronecker_dup_pre_post(sample.cov.inv,
multiplicator = 0.5)
} else {
I22 <- 0.5 * lav_matrix_duplication_pre_post(sample.cov.inv %x%
sample.cov.inv)

}

# fixed.x?
if (length(x.idx) > 0L) {
pstar.x <- lav_matrix_vech_which_idx(
Expand Down Expand Up @@ -470,8 +485,12 @@ lav_mvnorm_h1_inverted_information_observed <- function(
)
} else {
I11 <- sample.cov
I22 <- 2 * lav_matrix_duplication_ginv_pre_post(sample.cov %x% sample.cov)

if (lav_use_lavaanC()) {
I22 <- lavaanC::m_kronecker_dup_ginv_pre_post(sample.cov,
multiplicator = 2.0)
} else {
I22 <- 2 * lav_matrix_duplication_ginv_pre_post(sample.cov %x% sample.cov)
}
Gamma.NT <- lav_matrix_bdiag(I11, I22)
}

Expand Down
31 changes: 23 additions & 8 deletions R/lav_mvnorm_missing.R
Original file line number Diff line number Diff line change
Expand Up @@ -1037,8 +1037,11 @@ lav_mvnorm_missing_logl_hessian_samplestats <-

i11 <- S.inv
i21 <- lav_matrix_duplication_pre(tmp21 %x% S.inv)
i22 <- (1 / 2) * lav_matrix_duplication_pre_post(S.inv %x% tmp22)

if (lav_use_lavaanC()) {
i22 <- lavaanC::m_kronecker_dup_pre_post(S.inv, tmp22, 0.5)
} else {
i22 <- (1 / 2) * lav_matrix_duplication_pre_post(S.inv %x% tmp22)
}
H11 <- H11 + pat.freq * i11
H21 <- H21 + pat.freq * i21
H22 <- H22 + pat.freq * i22
Expand Down Expand Up @@ -1128,8 +1131,12 @@ lav_mvnorm_missing_information_expected <- function(Y = NULL,
S.inv <- matrix(0, P, P)
S.inv[var.idx, var.idx] <- sigma.inv

S2.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)

if (lav_use_lavaanC()) {
S2.inv <- lavaanC::m_kronecker_dup_pre_post(S.inv, multiplicator = 0.5)
} else {
S2.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
}

if (!is.null(wt)) {
FREQ <- sum(wt[Mp$case.idx[[p]]])
} else {
Expand Down Expand Up @@ -1360,8 +1367,12 @@ lav_mvnorm_missing_information_both <- function(Y = NULL,
}

if (information == "expected") {
S2.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)

if (lav_use_lavaanC()) {
S2.inv <- lavaanC::m_kronecker_dup_pre_post(S.inv, multiplicator = 0.5)
} else {
S2.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
}

I11 <- I11 + FREQ * S.inv
I22 <- I22 + FREQ * S2.inv
} else {
Expand All @@ -1379,8 +1390,12 @@ lav_mvnorm_missing_information_both <- function(Y = NULL,

i11 <- S.inv
i21 <- lav_matrix_duplication_pre(tmp21 %x% S.inv)
i22 <- (1 / 2) * lav_matrix_duplication_pre_post(S.inv %x% tmp22)

if (lav_use_lavaanC()) {
i22 <- lavaanC::m_kronecker_dup_pre_post(S.inv, tmp22, 0.5)
} else {
i22 <- (1 / 2) * lav_matrix_duplication_pre_post(S.inv %x% tmp22)
}

I11 <- I11 + pat.freq * i11
I21 <- I21 + pat.freq * i21
I22 <- I22 + pat.freq * i22
Expand Down
16 changes: 12 additions & 4 deletions R/lav_mvreg.R
Original file line number Diff line number Diff line change
Expand Up @@ -517,8 +517,12 @@ lav_mvreg_information_expected <- function(Y = NULL, # not used

# expected information
I11 <- res.cov.inv %x% sample.xx
I22 <- 0.5 * lav_matrix_duplication_pre_post(res.cov.inv %x% res.cov.inv)

if (lav_use_lavaanC()) {
I22 <- lavaanC::m_kronecker_dup_pre_post(res.cov.inv, multiplicator = 0.5)
} else {
I22 <- 0.5 * lav_matrix_duplication_pre_post(res.cov.inv %x% res.cov.inv)
}

lav_matrix_bdiag(I11, I22)
}

Expand Down Expand Up @@ -605,8 +609,12 @@ lav_mvreg_information_observed_samplestats <-
H12 <- t(H21)

AAA <- res.cov.inv %*% (2 * W.tilde - res.cov) %*% res.cov.inv
H22 <- (1 / 2) * lav_matrix_duplication_pre_post(res.cov.inv %x% AAA)

if (lav_use_lavaanC()) {
H22 <- lavaanC::m_kronecker_dup_pre_post(res.cov.inv, AAA, 0.5)
} else {
H22 <- (1 / 2) * lav_matrix_duplication_pre_post(res.cov.inv %x% AAA)
}

out <- rbind(
cbind(H11, H12),
cbind(H21, H22)
Expand Down
32 changes: 25 additions & 7 deletions R/lav_samplestats_gamma.R
Original file line number Diff line number Diff line change
Expand Up @@ -323,7 +323,11 @@ lav_samplestats_Gamma_NT <- function(Y = NULL, # should include
if (!conditional.x) {
# unconditional - stochastic x
if (!fixed.x) {
Gamma <- 2 * lav_matrix_duplication_ginv_pre_post(S %x% S)
if (lav_use_lavaanC()) {
Gamma <- lavaanC::m_kronecker_dup_ginv_pre_post(S, multiplicator = 2.0)
} else {
Gamma <- 2 * lav_matrix_duplication_ginv_pre_post(S %x% S)
}
if (meanstructure) {
Gamma <- lav_matrix_bdiag(S, Gamma)
}
Expand All @@ -346,8 +350,13 @@ lav_samplestats_Gamma_NT <- function(Y = NULL, # should include
# take difference
R <- S - YbarX.aug

Gamma.S <- 2 * lav_matrix_duplication_ginv_pre_post(S %x% S)
Gamma.R <- 2 * lav_matrix_duplication_ginv_pre_post(R %x% R)
if (lav_use_lavaanC()) {
Gamma.S <- lavaanC::m_kronecker_dup_ginv_pre_post(S, multiplicator = 2.0)
Gamma.R <- lavaanC::m_kronecker_dup_ginv_pre_post(R, multiplicator = 2.0)
} else {
Gamma.S <- 2 * lav_matrix_duplication_ginv_pre_post(S %x% S)
Gamma.R <- 2 * lav_matrix_duplication_ginv_pre_post(R %x% R)
}
Gamma <- Gamma.S - Gamma.R

if (meanstructure) {
Expand All @@ -368,8 +377,12 @@ lav_samplestats_Gamma_NT <- function(Y = NULL, # should include
B <- S[-x.idx, x.idx, drop = FALSE]
C <- S[x.idx, x.idx, drop = FALSE]
Cov.YbarX <- A - B %*% solve(C) %*% t(B)
Gamma <- 2 * lav_matrix_duplication_ginv_pre_post(Cov.YbarX %x% Cov.YbarX)

if (lav_use_lavaanC()) {
Gamma <- lavaanC::m_kronecker_dup_ginv_pre_post(Cov.YbarX, multiplicator = 2.0)
} else {
Gamma <- 2 * lav_matrix_duplication_ginv_pre_post(Cov.YbarX %x% Cov.YbarX)
}

if (meanstructure || slopestructure) {
MY <- M[-x.idx]
MX <- M[x.idx]
Expand Down Expand Up @@ -716,8 +729,13 @@ lav_samplestats_Gamma <- function(Y, # Y+X if cond!
# unbiased?
if (unbiased) {
# normal-theory Gamma (cov only)
GammaNT.cov <- 2 * lav_matrix_duplication_ginv_pre_post(COV %x% COV)

if (lav_use_lavaanC()) {
GammaNT.cov <- lavaanC::m_kronecker_dup_ginv_pre_post(COV,
multiplicator = 2.0)
} else {
GammaNT.cov <- 2 * lav_matrix_duplication_ginv_pre_post(COV %x% COV)
}

if (meanstructure) {
Gamma.cov <- Gamma[-(1:p), -(1:p), drop = FALSE]
Gamma.mean.cov <- Gamma[1:p, -(1:p), drop = FALSE]
Expand Down
22 changes: 17 additions & 5 deletions R/lav_samplestats_igamma.R
Original file line number Diff line number Diff line change
Expand Up @@ -75,16 +75,24 @@ lav_samplestats_Gamma_inverse_NT <- function(Y = NULL,
if (!conditional.x) {
# unconditional - stochastic x
if (!fixed.x) {
Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
if (lav_use_lavaanC()) {
Gamma.inv <- lavaanC::m_kronecker_dup_pre_post(S.inv, multiplicator = 0.5)
} else {
Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
}
if (meanstructure) {
Gamma.inv <- lav_matrix_bdiag(S.inv, Gamma.inv)
}

# unconditional - fixed x
} else {
# handle fixed.x = TRUE
Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)

if (lav_use_lavaanC()) {
Gamma.inv <- lavaanC::m_kronecker_dup_pre_post(S.inv, multiplicator = 0.5)
} else {
Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S.inv %x% S.inv)
}

# zero rows/cols corresponding with x/x combinations
nvar <- NROW(ICOV)
pstar <- nvar * (nvar + 1) / 2
Expand Down Expand Up @@ -112,8 +120,12 @@ lav_samplestats_Gamma_inverse_NT <- function(Y = NULL,

S11 <- S.inv[-x.idx, -x.idx, drop = FALSE]

Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S11 %x% S11)

if (lav_use_lavaanC()) {
Gamma.inv <- lavaanC::m_kronecker_dup_pre_post(S11, multiplicator = 0.5)
} else {
Gamma.inv <- 0.5 * lav_matrix_duplication_pre_post(S11 %x% S11)
}

if (meanstructure || slopestructure) {
C <- S[x.idx, x.idx, drop = FALSE]
MY <- M[-x.idx]
Expand Down

0 comments on commit 7fcecf6

Please sign in to comment.