From 4906e0da9b07e342cfe1c9cdc6b57b4c9dce2cf2 Mon Sep 17 00:00:00 2001 From: xuewei cao <36172337+xueweic@users.noreply.github.com> Date: Tue, 9 Sep 2025 00:18:37 -0400 Subject: [PATCH] important update To release the memory usage in LD-free mode, we pre-set LD = 1 instead of a diagonal matrix. However, using sum(LD) == 1 to check LD is pretty computational entensive. We changed to use length(LD) == 1 instead, which is saving a lot of computational time. --- R/colocboost_check_update_jk.R | 11 +++++++---- R/colocboost_inference.R | 8 ++++---- R/colocboost_init.R | 2 +- R/colocboost_update.R | 8 ++++---- 4 files changed, 16 insertions(+), 13 deletions(-) diff --git a/R/colocboost_check_update_jk.R b/R/colocboost_check_update_jk.R index a9917b0..623faa8 100644 --- a/R/colocboost_check_update_jk.R +++ b/R/colocboost_check_update_jk.R @@ -76,6 +76,8 @@ boost_check_update_jk_nofocal <- function(cb_model, cb_model_para, cb_data) { return(ifelse(length(jk_temp) == 1, jk_temp, sample(jk_temp, 1))) }) update_jk[c(1, pos.update + 1)] <- c(jk, jk_each) + + judge_each <- check_jk_jkeach(jk, jk_each, pos.update, model_update = model_update, @@ -493,12 +495,12 @@ get_LD_jk1_jk2 <- function(jk1, jk2, remain_jk = NULL) { if (!is.null(X)) { LD_temp <- suppressWarnings({ - cor(X[, c(jk1, jk2)]) + get_cormat(X[, c(jk1, jk2)]) }) LD_temp[which(is.na(LD_temp))] <- 0 LD_temp <- LD_temp[1, 2] } else if (!is.null(XtX)) { - if (sum(XtX) == 1){ + if (length(XtX) == 1){ LD_temp <- 0 } else { jk1.remain <- which(remain_jk == jk1) @@ -545,8 +547,9 @@ check_pair_jkeach <- function(jk_each, jk_equiv_loglik = 0.001) { data_update <- cb_data$data[pos.update] # -- check if jk_i ~ jk_j - change_each_pair <- matrix(NA, nrow = length(jk_each), ncol = length(jk_each)) + change_each_pair <- matrix(FALSE, nrow = length(jk_each), ncol = length(jk_each)) for (i in 1:length(jk_each)) { + jk_i <- jk_each[i] change_log_jk_i <- model_update[[i]]$change_loglike[jk_i] for (j in 1:length(jk_each)) { @@ -597,7 +600,7 @@ estimate_change_profile_res <- function(jk, } else { xty <- XtY / scaling_factor } - if (sum(xtx) == 1){ + if (length(xtx) == 1){ rtr <- yty - 2 * sum(beta_k * xty) + sum(beta_k^2) } else { rtr <- yty - 2 * sum(beta_k * xty) + sum((xtx %*% as.matrix(beta_k)) * beta_k) diff --git a/R/colocboost_inference.R b/R/colocboost_inference.R index 5226db0..dc53908 100644 --- a/R/colocboost_inference.R +++ b/R/colocboost_inference.R @@ -241,7 +241,7 @@ check_null_post <- function(cb_obj, } else { xty <- XtY / scaling_factor } - if (sum(xtx) == 1){ + if (length(xtx) == 1){ (yty - 2 * sum(cs_beta * xty) + sum(cs_beta^2)) * adj_dep } else { (yty - 2 * sum(cs_beta * xty) + sum((xtx %*% as.matrix(cs_beta)) * cs_beta)) * adj_dep @@ -293,14 +293,14 @@ check_null_post <- function(cb_obj, if (length(miss_idx) != 0) { xty <- XtY[-miss_idx] / scaling.factor res.tmp <- rep(0, length(XtY)) - if (sum(xtx) == 1){ + if (length(xtx) == 1){ res.tmp[-miss_idx] <- xty - cs_beta[-miss_idx] / beta_scaling } else { res.tmp[-miss_idx] <- xty - xtx %*% (cs_beta[-miss_idx] / beta_scaling) } } else { xty <- XtY / scaling.factor - if (sum(xtx) == 1){ + if (length(xtx) == 1){ res.tmp <- xty - (cs_beta / beta_scaling) } else { res.tmp <- xty - xtx %*% (cs_beta / beta_scaling) @@ -508,7 +508,7 @@ get_cos_evidence <- function(cb_obj, coloc_out, data_info) { } else { xty <- XtY / scaling_factor } - if (sum(xtx) == 1){ + if (length(xtx) == 1){ cos_profile <- (yty - 2 * sum(cs_beta * xty) + sum(cs_beta^2)) * adj_dep } else { cos_profile <- (yty - 2 * sum(cs_beta * xty) + sum((xtx %*% as.matrix(cs_beta)) * cs_beta)) * adj_dep diff --git a/R/colocboost_init.R b/R/colocboost_init.R index bd278ab..8685b01 100644 --- a/R/colocboost_init.R +++ b/R/colocboost_init.R @@ -396,7 +396,7 @@ get_correlation <- function(X = NULL, res = NULL, XtY = NULL, N = NULL, Xtr <- res / scaling_factor XtY <- XtY / scaling_factor } - if (sum(XtX) == 1){ + if (length(XtX) == 1){ var_r <- YtY - 2 * sum(beta_k * XtY) + sum(beta_k^2) } else { var_r <- YtY - 2 * sum(beta_k * XtY) + sum((XtX %*% as.matrix(beta_k)) * beta_k) diff --git a/R/colocboost_update.R b/R/colocboost_update.R index 5858a3e..6100828 100644 --- a/R/colocboost_update.R +++ b/R/colocboost_update.R @@ -111,7 +111,7 @@ colocboost_update <- function(cb_model, cb_model_para, cb_data) { if (length(cb_data$data[[i]]$variable_miss) != 0) { beta <- cb_model[[i]]$beta[-cb_data$data[[i]]$variable_miss] / beta_scaling xty <- cb_data$data[[i]]$XtY[-cb_data$data[[i]]$variable_miss] - if (sum(xtx) == 1){ + if (length(xtx) == 1){ cb_model[[i]]$res[-cb_data$data[[i]]$variable_miss] <- xty - scaling_factor * beta } else { cb_model[[i]]$res[-cb_data$data[[i]]$variable_miss] <- xty - scaling_factor * xtx %*% beta @@ -120,7 +120,7 @@ colocboost_update <- function(cb_model, cb_model_para, cb_data) { } else { beta <- cb_model[[i]]$beta / beta_scaling xty <- cb_data$data[[i]]$XtY - if (sum(xtx) == 1){ + if (length(xtx) == 1){ cb_model[[i]]$res <- xty - scaling_factor * beta } else { cb_model[[i]]$res <- xty - scaling_factor * xtx %*% beta @@ -129,7 +129,7 @@ colocboost_update <- function(cb_model, cb_model_para, cb_data) { # - profile-loglikelihood yty <- cb_data$data[[i]]$YtY / scaling_factor xty <- xty / scaling_factor - if (sum(xtx) == 1){ + if (length(xtx) == 1){ profile_log <- (yty - 2 * sum(beta * xty) + sum(beta^2)) * adj_dep } else { profile_log <- (yty - 2 * sum(beta * xty) + sum((xtx %*% as.matrix(beta)) * beta)) * adj_dep @@ -152,7 +152,7 @@ get_LD_jk <- function(jk1, X = NULL, XtX = NULL, N = NULL, remain_idx = NULL, P } else if (!is.null(XtX)) { jk1.remain <- which(remain_idx == jk1) corr <- rep(0, P) - if (sum(XtX) == 1 | length(jk1.remain)==0){ + if (length(XtX) == 1 | length(jk1.remain)==0){ corr[remain_idx] <- 1 } else { corr[remain_idx] <- XtX[, jk1.remain]