Skip to content

Commit

Permalink
minor speedups all around
Browse files Browse the repository at this point in the history
  • Loading branch information
@jr-leary7
jr-leary7 committed Oct 10, 2024
1 parent 4f2e2d8 commit 75c7d93
Show file tree
Hide file tree
Showing 11 changed files with 47 additions and 36 deletions.
2 changes: 2 additions & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,9 @@ importFrom(MASS,theta.mm)
importFrom(Matrix,Matrix)
importFrom(Matrix,bdiag)
importFrom(Matrix,colSums)
importFrom(Matrix,diag)
importFrom(Matrix,rowMeans)
importFrom(Matrix,solve)
importFrom(Matrix,t)
importFrom(Rcpp,sourceCpp)
importFrom(bigstatsr,as_FBM)
Expand Down
4 changes: 2 additions & 2 deletions R/RcppExports.R
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

eigenMapMatrixInvert <- function(A, n_cores) {
eigenMapMatrixInvert <- function(A, n_cores = 1L) {
.Call(`_scLANE_eigenMapMatrixInvert`, A, n_cores)
}

eigenMapMatMult <- function(A, B, n_cores) {
eigenMapMatMult <- function(A, B, n_cores = 1L) {
.Call(`_scLANE_eigenMapMatMult`, A, B, n_cores)
}

4 changes: 2 additions & 2 deletions R/marge2.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ marge2 <- function(X_pred = NULL,
q <- ncol(X_pred) # Number of predictor variables
B <- as.matrix(rep(1, NN)) # Start with the intercept model.
# estimate "known" theta for GEE models
if (is.gee) {
if (is.gee || !is.glmm) {
theta_hat <- MASS::theta.mm(y = Y,
mu = mean(Y),
dfr = NN - 1)
Expand Down Expand Up @@ -840,7 +840,7 @@ marge2 <- function(X_pred = NULL,
data = model_df,
method = "glm.fit2",
link = log,
init.theta = 1,
init.theta = theta_hat,
y = FALSE,
model = FALSE)
}
Expand Down
18 changes: 9 additions & 9 deletions R/score_fun_gee.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#' @author David I. Warton
#' @author Jack Leary
#' @importFrom stats lm.fit
#' @importFrom Matrix diag t solve
#' @importFrom MASS ginv
#' @description Calculate the score statistic for a GEE model.
#' @param Y The response variable. Defaults to NULL.
Expand Down Expand Up @@ -51,22 +52,21 @@ score_fun_gee <- function(Y = NULL,
B.est <- matrix(0, nrow = p, ncol = 1)
Sigma22 <- matrix(0, nrow = p, ncol = p)
J21 <- Sigma21 <- matrix(0, nrow = p, ncol = p1)

for (i in seq(N)) {
k <- sum(n_vec[seq(i)])
VS.est_i <- VS.est_list[[i]]
AWA.est_i <- AWA.est_list[[i]]
J2_i <- J2_list[[i]]
Sigma2_i <- Sigma2_list[[i]]
D.est_i <- eigenMapMatMult(A = diag((mu.est[(sum(n_vec1[seq(i)]) + 1):k]), nrow = n_vec[i], ncol = n_vec[i]),
D.est_i <- eigenMapMatMult(A = Matrix::diag((mu.est[(sum(n_vec1[seq(i)]) + 1):k]), nrow = n_vec[i], ncol = n_vec[i]),
B = XA[(sum(n_vec1[seq(i)]) + 1):k, ],
n_cores = 1)
D_est_i_transpose <- t(D.est_i)
D_est_i_transpose <- Matrix::t(D.est_i)
J21 <- J21 + eigenMapMatMult(A = D_est_i_transpose,
B = t(J2_i),
B = Matrix::t(J2_i),
n_cores = 1)
Sigma21 <- Sigma21 + eigenMapMatMult(A = D_est_i_transpose,
B = t(Sigma2_i),
B = Matrix::t(Sigma2_i),
n_cores = 1)
B.est <- B.est + eigenMapMatMult(A = D_est_i_transpose,
B = VS.est_i,
Expand All @@ -82,12 +82,12 @@ score_fun_gee <- function(Y = NULL,
B = J11.inv,
n_cores = 1)
temp_prod_1 <- eigenMapMatMult(A = temp_prod_1,
B = t(Sigma21),
B = Matrix::t(Sigma21),
n_cores = 1)
temp_prod_2 <- eigenMapMatMult(A = Sigma21,
B = J11.inv,
n_cores = 1)
J21_transpose <- t(J21)
J21_transpose <- Matrix::t(J21)
temp_prod_2 <- eigenMapMatMult(A = temp_prod_2,
B = J21_transpose,
n_cores = 1)
Expand All @@ -98,11 +98,11 @@ score_fun_gee <- function(Y = NULL,
B = J21_transpose,
n_cores = 1)
Sigma <- Sigma22 - temp_prod_1 - temp_prod_2 + temp_prod_3
Sigma_inv <- try({ eigenMapMatrixInvert(Sigma, n_cores = 1) }, silent = TRUE)
Sigma_inv <- try({ Matrix::solve(Sigma) }, silent = TRUE)
if (inherits(Sigma_inv, "try-error")) {
Sigma_inv <- MASS::ginv(Sigma)
}
temp_prod <- eigenMapMatMult(A = t(B.est),
temp_prod <- eigenMapMatMult(A = Matrix::t(B.est),
B = Sigma_inv,
n_cores = 1)
score <- eigenMapMatMult(A = temp_prod,
Expand Down
11 changes: 6 additions & 5 deletions R/score_fun_glm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#' @author David I. Warton
#' @author Jack Leary
#' @importFrom stats lm.fit
#' @importFrom Matrix t
#' @importFrom MASS ginv
#' @description Calculate the score statistic for a GLM model.
#' @param Y The response variable. Defaults to NULL.
Expand Down Expand Up @@ -40,28 +41,28 @@ score_fun_glm <- function(Y = NULL,
B_list_i <- eigenMapMatMult(A = B1_list_i,
B = XA,
n_cores = 1)
D_list_i <- eigenMapMatMult(A = t(XA),
D_list_i <- eigenMapMatMult(A = Matrix::t(XA),
B = (XA * c(mu.est^2 / V.est)),
n_cores = 1)
temp_prod <- eigenMapMatMult(A = t(B_list_i),
temp_prod <- eigenMapMatMult(A = Matrix::t(B_list_i),
B = A_list_i,
n_cores = 1)
temp_prod <- eigenMapMatMult(A = temp_prod,
B = B_list_i,
n_cores = 1)
inv.XVX_22 <- D_list_i - temp_prod
B.est <- eigenMapMatMult(A = t(mu.est * VS.est_i),
B.est <- eigenMapMatMult(A = Matrix::t(mu.est * VS.est_i),
B = XA,
n_cores = 1)
XVX_22 <- try({ eigenMapMatrixInvert(inv.XVX_22, n_cores = 1) }, silent = TRUE)
XVX_22 <- try({ solve(inv.XVX_22) }, silent = TRUE)
if (inherits(XVX_22, "try-error")) {
XVX_22 <- MASS::ginv(inv.XVX_22)
}
temp_prod <- eigenMapMatMult(A = B.est,
B = XVX_22,
n_cores = 1)
score <- eigenMapMatMult(A = temp_prod,
B = t(B.est),
B = Matrix::t(B.est),
n_cores = 1)
}
res <- list(score = score)
Expand Down
25 changes: 13 additions & 12 deletions R/stat_out_score_gee_null.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#' @importFrom geeM geem
#' @importFrom MASS negative.binomial ginv
#' @importFrom stats fitted.values
#' @importFrom Matrix diag t solve
#' @description A function that calculates parts of the score statistic for GEEs only (it is used for the full path for forward selection).
#' @param Y The response variable Defaults to NULL.
#' @param B_null The design matrix matrix under the null model Defaults to NULL.
Expand Down Expand Up @@ -50,10 +51,10 @@ stat_out_score_gee_null <- function(Y = NULL,
k <- sum(n_vec[seq(i)])
# set up working correlation matrix structure
if (cor.structure == "independence") {
R_alpha <- diag(1, n_vec[i], n_vec[i])
R_alpha <- Matrix::diag(1, n_vec[i], n_vec[i])
} else if (cor.structure == "exchangeable") {
R_alpha <- matrix(alpha_est, nrow = n_vec[i], ncol = n_vec[i])
R_alpha <- R_alpha + diag(c(1 - alpha_est), nrow = n_vec[i], ncol = n_vec[i])
R_alpha <- R_alpha + Matrix::diag(c(1 - alpha_est), nrow = n_vec[i], ncol = n_vec[i])
} else if (cor.structure == "ar1") {
R_alpha <- alpha_est^outer(seq(n_vec[i]), seq(n_vec[i]), \(x, y) abs(x - y))
} else {
Expand All @@ -62,32 +63,32 @@ stat_out_score_gee_null <- function(Y = NULL,
# compute iteration values for each statistic
temp_seq_n <- (sum(n_vec1[seq(i)]) + 1):k
mu_est_sub <- mu_est[temp_seq_n]
diag_sqrt_V_est <- diag(sqrt(V_est[temp_seq_n]),
nrow = n_vec[i],
ncol = n_vec[i])
diag_sqrt_V_est <- Matrix::diag(sqrt(V_est[temp_seq_n]),
nrow = n_vec[i],
ncol = n_vec[i])
temp_prod <- eigenMapMatMult(A = diag_sqrt_V_est,
B = R_alpha,
n_cores = 1)
V_est_i <- eigenMapMatMult(A = temp_prod,
B = diag_sqrt_V_est,
n_cores = 1)
V_est_i_inv <- eigenMapMatrixInvert(A = V_est_i, n_cores = 1)
S_est_i <- t(Y)[temp_seq_n] - mu_est_sub
S_est_i <- Matrix::t(Y)[temp_seq_n] - mu_est_sub
temp_prod <- eigenMapMatMult(A = S_est_i,
B = t(S_est_i),
B = Matrix::t(S_est_i),
n_cores = 1)
temp_prod <- eigenMapMatMult(A = V_est_i_inv,
B = temp_prod,
n_cores = 1)
AWA_est_i <- eigenMapMatMult(A = temp_prod,
B = V_est_i_inv,
n_cores = 1)
D_est_i <- eigenMapMatMult(A = diag(mu_est_sub,
nrow = n_vec[i],
ncol = n_vec[i]),
D_est_i <- eigenMapMatMult(A = Matrix::diag(mu_est_sub,
nrow = n_vec[i],
ncol = n_vec[i]),
B = B_null[temp_seq_n, ],
n_cores = 1)
D_est_i_transpose <- t(D_est_i)
D_est_i_transpose <- Matrix::t(D_est_i)
temp_prod <- eigenMapMatMult(A = D_est_i_transpose,
B = V_est_i_inv,
n_cores = 1)
Expand Down Expand Up @@ -116,7 +117,7 @@ stat_out_score_gee_null <- function(Y = NULL,
if (nrow(J11) == 1 && ncol(J11) == 1) {
J11_inv <- 1 / J11
} else {
J11_inv <- try({ eigenMapMatrixInvert(J11, n_cores = 1) }, silent = TRUE)
J11_inv <- try({ Matrix::solve(J11) }, silent = TRUE)
if (inherits(J11_inv, "try-error")) {
J11_inv <- MASS::ginv(J11)
}
Expand Down
9 changes: 5 additions & 4 deletions R/stat_out_score_glm_null.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#' @author Jack Leary
#' @importFrom gamlss gamlss
#' @importFrom stats fitted.values
#' @importFrom Matrix diag t
#' @importFrom MASS ginv
#' @description A function that calculates parts of the score statistic for GLMs only (it is used for the full path for forward selection).
#' @param Y : The response variable Defaults to NULL.
Expand All @@ -28,8 +29,8 @@ stat_out_score_glm_null <- function(Y = NULL, B_null = NULL) {
mu_est <- as.matrix(stats::fitted.values(ests))
V_est <- mu_est * (1 + mu_est * sigma_est) # Type I NB variance = mu (1 + mu * sigma); sigma = 1 / theta
VS_est_list <- (Y - mu_est) / V_est
mu_V_diag <- diag(c(mu_est^2 / V_est))
temp_prod <- eigenMapMatMult(A = t(B_null),
mu_V_diag <- Matrix::diag(c(mu_est^2 / V_est))
temp_prod <- eigenMapMatMult(A = Matrix::t(B_null),
B = mu_V_diag,
n_cores = 1)
A_list_inv <- eigenMapMatMult(A = temp_prod,
Expand All @@ -38,12 +39,12 @@ stat_out_score_glm_null <- function(Y = NULL, B_null = NULL) {
if (ncol(A_list_inv) == 1 && nrow(A_list_inv) == 1) {
A_list <- 1 / A_list_inv
} else {
A_list <- try({ eigenMapMatrixInvert(A_list_inv, n_cores = 1) }, silent = TRUE)
A_list <- try({ solve(A_list_inv) }, silent = TRUE)
if (inherits(A_list, "try-error")) {
A_list <- MASS::ginv(A_list_inv)
}
}
B1_list <- eigenMapMatMult(A = t(B_null),
B1_list <- eigenMapMatMult(A = Matrix::t(B_null),
B = mu_V_diag,
n_cores = 1)
res <- list(VS.est_list = VS_est_list,
Expand Down
3 changes: 3 additions & 0 deletions src/Makevars
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,6 @@
## With R 3.1.0 or later, you can uncomment the following line to tell R to
## enable compilation with C++11 (or even C++14) where available
CXX_STD = CXX11

## enable CPU speedups
CXXFLAGS += -O3 -march=native
3 changes: 3 additions & 0 deletions src/Makevars.win
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,6 @@
## With R 3.1.0 or later, you can uncomment the following line to tell R to
## enable compilation with C++11 (or even C++14) where available
CXX_STD = CXX11

## enable CPU speedups
CXXFLAGS += -O3 -march=native
2 changes: 1 addition & 1 deletion src/eigenMapInvert.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::export]]

SEXP eigenMapMatrixInvert(const Eigen::Map<Eigen::MatrixXd> A, int n_cores){
SEXP eigenMapMatrixInvert(const Eigen::Map<Eigen::MatrixXd> A, int n_cores = 1){
Eigen::setNbThreads(n_cores);
Eigen::MatrixXd B = A.llt().solve(Eigen::MatrixXd::Identity(A.rows(), A.cols()));
return Rcpp::wrap(B);
Expand Down
2 changes: 1 addition & 1 deletion src/eigenMapMatMult.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@

SEXP eigenMapMatMult(const Eigen::Map<Eigen::MatrixXd> A,
const Eigen::Map<Eigen::MatrixXd> B,
int n_cores){
int n_cores = 1) {
Eigen::setNbThreads(n_cores);
Eigen::MatrixXd C;
C.noalias() = A * B;
Expand Down

0 comments on commit 75c7d93

Please sign in to comment.