diff --git a/DESCRIPTION b/DESCRIPTION index 40d065e..9b22a0d 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,8 +2,8 @@ Package: rMVP Type: Package Title: Memory-Efficient, Visualize-Enhanced, Parallel-Accelerated GWAS Tool -Version: 1.3.0 -Date: 2024-12-12 +Version: 1.3.5 +Date: 2024-12-30 Authors@R: c( person("Lilin", "Yin", role = "aut"), person("Haohao", "Zhang", role = "aut"), person("Zhenshuang", "Tang", role = "aut"), diff --git a/R/MVP.FarmCPU.r b/R/MVP.FarmCPU.r index ba28518..af540dc 100755 --- a/R/MVP.FarmCPU.r +++ b/R/MVP.FarmCPU.r @@ -74,6 +74,8 @@ if(!is.big.matrix(geno)) stop("genotype should be in 'big.matrix' format.") if(sum(is.na(phe[, 2])) != 0) stop("NAs are not allowed in phenotype.") + if(nrow(map) != ncol(geno) & nrow(map) != nrow(geno)) stop("The number of markers in genotype and map doesn't match!") + mrk_bycol <- ncol(geno) == nrow(map) if(is.null(ind_idx)){ if(nrow(phe) != ncol(geno) && nrow(phe) != nrow(geno)) stop("number of individuals does not match in phenotype and genotype.") n <- ifelse(nrow(phe) == ncol(geno), ncol(geno), nrow(geno)) @@ -242,7 +244,7 @@ theCV=cbind(CV,myRemove$bin) } - myGLM=FarmCPU.LM(y=phe[,2],GDP=geno,GDP_index=ind_idx,w=theCV,maxLine=maxLine,ncpus=ncpus,npc=npc, verbose=verbose) + myGLM=FarmCPU.LM(y=phe[,2],GDP=geno,GDP_index=ind_idx,GDP_mrk_bycol=mrk_bycol,w=theCV,maxLine=maxLine,ncpus=ncpus,npc=npc, verbose=verbose) if(!is.null(seqQTN)){ if(ncol(myGLM$P) != (npc + length(seqQTN) + 1)) stop("wrong dimensions.") } @@ -499,7 +501,7 @@ FarmCPU.BIN <- if(is.null(P)) return(list(bin=NULL,binmap=NULL,seqQTN=NULL)) - if(nrow(Y) == nrow(GDP)){ + if(nrow(GM) == ncol(GDP)){ if(is.null(GDP_index)) GDP_index=seq(1, nrow(GDP)) }else{ if(is.null(GDP_index)) GDP_index=seq(1, ncol(GDP)) @@ -558,7 +560,7 @@ FarmCPU.BIN <- GP=cbind(GM,P,NA,NA,NA) mySpecify=FarmCPU.Specify(GI=GM,GP=GP,bin.size=bin,inclosure.size=inc) seqQTN=which(mySpecify$index==TRUE) - if(nrow(Y) == nrow(GDP)){ + if(nrow(GM) == ncol(GDP)){ GK=GDP[GDP_index, seqQTN] }else{ GK=t(GDP[seqQTN, GDP_index]) @@ -606,7 +608,7 @@ FarmCPU.BIN <- GP=cbind(GM,P,NA,NA,NA) mySpecify=FarmCPU.Specify(GI=GM,GP=GP,bin.size=bin,inclosure.size=inc) seqQTN=which(mySpecify$index==TRUE) - if(nrow(Y) == nrow(GDP)){ + if(nrow(GM) == ncol(GDP)){ GK=GDP[GDP_index, seqQTN] }else{ GK=t(GDP[seqQTN, GDP_index]) @@ -647,7 +649,7 @@ FarmCPU.BIN <- GP=cbind(GM,P,NA,NA,NA) mySpecify=FarmCPU.Specify(GI=GM,GP=GP,bin.size=bin[ii],inclosure.size=inc[ii]) seqQTN=which(mySpecify$index==TRUE) - if(nrow(Y) == nrow(GDP)){ + if(nrow(GM) == ncol(GDP)){ GK=GDP[GDP_index, seqQTN] }else{ GK=t(GDP[seqQTN, GDP_index]) @@ -775,6 +777,7 @@ FarmCPU.Specify <- #' @param w covariates, n by c matrix, n is sample size, c is number of covariates #' @param GDP genotype #' @param GDP_index index of effective genotyped individuals +#' @param GDP_mrk_bycol whether the markers are stored by columns in genotype (i.e. M is a n by m matrix) #' @param maxLine the number of markers handled at a time, smaller value would reduce the memory cost #' @param ncpus number of threads used for parallele computation #' @param npc number of covariates without pseudo QTNs @@ -787,7 +790,7 @@ FarmCPU.Specify <- #' #' @keywords internal FarmCPU.LM <- - function(y, w=NULL, GDP, GDP_index=NULL, maxLine=5000, ncpus=2, npc=0, verbose=TRUE){ + function(y, w=NULL, GDP, GDP_index=NULL, GDP_mrk_bycol=TRUE, maxLine=5000, ncpus=2, npc=0, verbose=TRUE){ #print("FarmCPU.LM started") if(is.null(y)) return(NULL) if(is.null(GDP)) return(NULL) @@ -863,7 +866,7 @@ FarmCPU.LM <- # if(ncol(w) == 50) write.csv(P, "P.csv") mkl_env({ - results <- glm_c(y=y, X=w, iXX = wwi, GDP@address, geno_ind=GDP_index, step=maxLine, verbose=verbose, threads=ncpus) + results <- glm_c(y=y, X=w, iXX = wwi, GDP@address, geno_ind=GDP_index, marker_bycol=GDP_mrk_bycol, step=maxLine, verbose=verbose, threads=ncpus) }) return(list(P=results[ ,-c(1:3), drop=FALSE], betapred=betapred, sepred=sepred, B=results[ , 1, drop=FALSE], S=results[ , 2, drop=FALSE])) # return(list(P=P, betapred=betapred, B=as.matrix(B), S=S)) diff --git a/R/MVP.GLM.r b/R/MVP.GLM.r index 044e6da..b9ec4b3 100644 --- a/R/MVP.GLM.r +++ b/R/MVP.GLM.r @@ -23,6 +23,7 @@ #' @param CV Covariance, design matrix(n * x) for the fixed effects #' @param ind_idx the index of effective genotyped individuals #' @param mrk_idx the index of effective markers used in analysis +#' @param mrk_bycol whether the markers are stored by columns in genotype (i.e. M is a n by m matrix) #' @param maxLine the number of markers handled at a time, smaller value would reduce the memory cost #' @param cpu number of cpus used for parallel computation #' @param verbose whether to print detail. @@ -52,6 +53,7 @@ function( CV=NULL, ind_idx=NULL, mrk_idx=NULL, + mrk_bycol=TRUE, maxLine=5000, cpu=1, verbose=TRUE @@ -60,6 +62,11 @@ function( if(is.null(ind_idx)){ if(nrow(phe) != ncol(geno) && nrow(phe) != nrow(geno)) stop("number of individuals does not match in phenotype and genotype.") n <- ifelse(nrow(phe) == ncol(geno), ncol(geno), nrow(geno)) + if(mrk_bycol){ + if(nrow(phe) != nrow(geno)) stop("set the argument 'mrk_bycol=FALSE'?") + }else{ + if(nrow(phe) != ncol(geno)) stop("set the argument 'mrk_bycol=TRUE'?") + } }else{ n <- length(ind_idx) if(nrow(phe) != n) stop("number of individuals does not match in phenotype and genotype.") @@ -84,7 +91,7 @@ function( logging.log("scanning...\n", verbose = verbose) mkl_env({ - results <- glm_c(y = ys, X = X0, iXX = iX0X0, geno@address, geno_ind = ind_idx, marker_ind = mrk_idx, step = maxLine, verbose = verbose, threads = cpu) + results <- glm_c(y = ys, X = X0, iXX = iX0X0, geno@address, geno_ind = ind_idx, marker_ind = mrk_idx, marker_bycol = mrk_bycol, step = maxLine, verbose = verbose, threads = cpu) }) return(results[, c(1, 2, ncol(results))]) diff --git a/R/MVP.MLM.r b/R/MVP.MLM.r index 1c2ad02..a700027 100644 --- a/R/MVP.MLM.r +++ b/R/MVP.MLM.r @@ -22,6 +22,7 @@ #' @param CV covariates #' @param ind_idx the index of effective genotyped individuals #' @param mrk_idx the index of effective markers used in analysis +#' @param mrk_bycol whether the markers are stored by columns in genotype (i.e. M is a n by m matrix) #' @param REML a list that contains ve and vg #' @param maxLine the number of markers handled at a time, smaller value would reduce the memory cost #' @param cpu number of cpus used for parallel computation @@ -59,6 +60,7 @@ function( CV=NULL, ind_idx=NULL, mrk_idx=NULL, + mrk_bycol=TRUE, REML=NULL, maxLine=5000, cpu=1, @@ -69,6 +71,11 @@ function( if(is.null(ind_idx)){ if(nrow(phe) != ncol(geno) && nrow(phe) != nrow(geno)) stop("number of individuals does not match in phenotype and genotype.") n <- ifelse(nrow(phe) == ncol(geno), ncol(geno), nrow(geno)) + if(mrk_bycol){ + if(nrow(phe) != nrow(geno)) stop("set the argument 'mrk_bycol=FALSE'?") + }else{ + if(nrow(phe) != ncol(geno)) stop("set the argument 'mrk_bycol=TRUE'?") + } }else{ n <- length(ind_idx) if(nrow(phe) != n) stop("number of individuals does not match in phenotype and genotype.") @@ -125,7 +132,7 @@ function( logging.log("scanning...\n", verbose = verbose) mkl_env({ - results <- mlm_c(y = ys, X = X0, U = U, vgs = vgs, geno@address, ind_idx, mrk_idx, step = maxLine, verbose = verbose, threads = cpu) + results <- mlm_c(y = ys, X = X0, U = U, vgs = vgs, geno@address, geno_ind = ind_idx, marker_ind = mrk_idx, marker_bycol = mrk_bycol, step = maxLine, verbose = verbose, threads = cpu) }) return(results) diff --git a/R/MVP.r b/R/MVP.r index e453aa0..99d3b65 100755 --- a/R/MVP.r +++ b/R/MVP.r @@ -156,11 +156,11 @@ function(phe, geno, map, K=NULL, nPC.GLM=NULL, nPC.MLM=NULL, nPC.FarmCPU=NULL, logging.log(paste("Markers are detected to be stored by", ifelse(MrkByCol, "column", "row")), "\n", verbose = verbose) logging.log("Analyzed trait:", colnames(phe)[2], "\n", verbose = verbose) logging.log("Number of threads used:", ncpus, "\n", verbose = verbose) - hpclib <- grepl("mkl", sessionInfo()$LAPACK) | grepl("openblas", sessionInfo()$LAPACK) | !inherits(try(Revo.version, silent=TRUE),"try-error") + hpclib <- grepl("mkl", sessionInfo()$LAPACK) | grepl("openblas", sessionInfo()$LAPACK) | eval(parse(text = "!inherits(try(Revo.version,silent=TRUE),'try-error')")) if(!hpclib){ logging.log("No high performance math library detected! The computational efficiency would be greatly reduced\n", verbose = verbose) }else{ - if(grepl("mkl", sessionInfo()$LAPACK) | !inherits(try(Revo.version, silent=TRUE),"try-error")){ + if(grepl("mkl", sessionInfo()$LAPACK) | eval(parse(text = "!inherits(try(Revo.version,silent=TRUE),'try-error')"))){ logging.log("Math Kernel Library is detected, nice job!\n", verbose = verbose) }else{ logging.log("OpenBLAS Library is detected, nice job!\n", verbose = verbose) @@ -190,7 +190,7 @@ function(phe, geno, map, K=NULL, nPC.GLM=NULL, nPC.MLM=NULL, nPC.FarmCPU=NULL, } } - logging.log("Check if NAs exist in genotype...", "\n", verbose = verbose) + logging.log("Check if NAs exist in genotype", "\n", verbose = verbose) if(hasNA(geno@address, mrkbycol = MrkByCol, geno_ind = seqTaxa, threads = ncpus)) stop("NA is not allowed in genotype, use 'MVP.Data.impute' to impute.") logging.log("Calculate allele frequency...", "\n", verbose = verbose) @@ -340,7 +340,7 @@ function(phe, geno, map, K=NULL, nPC.GLM=NULL, nPC.MLM=NULL, nPC.FarmCPU=NULL, logging.log("-------------------------GWAS Start-------------------------", "\n", verbose = verbose) if (glm.run) { logging.log("General Linear Model (GLM) Start...", "\n", verbose = verbose) - glm.results <- MVP.GLM(phe=phe, geno=geno, CV=CV.GLM, ind_idx=seqTaxa, mrk_idx=geno_marker_index, maxLine = maxLine, cpu=ncpus, verbose = verbose);gc() + glm.results <- MVP.GLM(phe=phe, geno=geno, CV=CV.GLM, ind_idx=seqTaxa, mrk_idx=geno_marker_index, mrk_bycol = MrkByCol, maxLine = maxLine, cpu=ncpus, verbose = verbose);gc() colnames(glm.results) <- c("Effect", "SE", paste(colnames(phe)[2],"GLM",sep=".")) z = glm.results[, 1]/glm.results[, 2] lambda = median(z^2, na.rm=TRUE)/qchisq(1/2, df = 1,lower.tail=FALSE) @@ -355,7 +355,7 @@ function(phe, geno, map, K=NULL, nPC.GLM=NULL, nPC.MLM=NULL, nPC.FarmCPU=NULL, if (mlm.run) { logging.log("Mixed Linear Model (MLM) Start...", "\n", verbose = verbose) - mlm.results <- MVP.MLM(phe=phe, geno=geno, K=K, eigenK=eigenK, CV=CV.MLM, ind_idx=seqTaxa, mrk_idx=geno_marker_index, maxLine = maxLine, cpu=ncpus, vc.method=vc.method, verbose = verbose);gc() + mlm.results <- MVP.MLM(phe=phe, geno=geno, K=K, eigenK=eigenK, CV=CV.MLM, ind_idx=seqTaxa, mrk_idx=geno_marker_index, mrk_bycol = MrkByCol, maxLine = maxLine, cpu=ncpus, vc.method=vc.method, verbose = verbose);gc() colnames(mlm.results) <- c("Effect", "SE", paste(colnames(phe)[2],"MLM",sep=".")) z = mlm.results[, 1]/mlm.results[, 2] lambda = median(z^2, na.rm=TRUE)/qchisq(1/2, df = 1,lower.tail=FALSE) diff --git a/R/RcppExports.R b/R/RcppExports.R index 6cda37c..5d9fd19 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -5,12 +5,12 @@ getRow <- function(pBigMat, row) { .Call(`_rMVP_getRow`, pBigMat, row) } -glm_c <- function(y, X, iXX, pBigMat, geno_ind = NULL, marker_ind = NULL, step = 10000L, verbose = TRUE, threads = 0L) { - .Call(`_rMVP_glm_c`, y, X, iXX, pBigMat, geno_ind, marker_ind, step, verbose, threads) +glm_c <- function(y, X, iXX, pBigMat, geno_ind = NULL, marker_ind = NULL, marker_bycol = TRUE, step = 10000L, verbose = TRUE, threads = 0L) { + .Call(`_rMVP_glm_c`, y, X, iXX, pBigMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads) } -mlm_c <- function(y, X, U, vgs, pBigMat, geno_ind = NULL, marker_ind = NULL, step = 10000L, verbose = TRUE, threads = 0L) { - .Call(`_rMVP_mlm_c`, y, X, U, vgs, pBigMat, geno_ind, marker_ind, step, verbose, threads) +mlm_c <- function(y, X, U, vgs, pBigMat, geno_ind = NULL, marker_ind = NULL, marker_bycol = TRUE, step = 10000L, verbose = TRUE, threads = 0L) { + .Call(`_rMVP_mlm_c`, y, X, U, vgs, pBigMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads) } vcf_parser_map <- function(vcf_file, out) { @@ -61,16 +61,8 @@ hasNA <- function(pBigMat, mrkbycol = TRUE, geno_ind = NULL, marker_ind = NULL, .Call(`_rMVP_hasNA`, pBigMat, mrkbycol, geno_ind, marker_ind, threads) } -BigRowMean <- function(pBigMat, mrkbycol = TRUE, threads = 0L, geno_ind = NULL) { - .Call(`_rMVP_BigRowMean`, pBigMat, mrkbycol, threads, geno_ind) -} - -kin_cal_m <- function(pBigMat, threads = 0L, verbose = TRUE) { - .Call(`_rMVP_kin_cal_m`, pBigMat, threads, verbose) -} - -kin_cal_s <- function(pBigMat, threads = 0L, mkl = FALSE, verbose = TRUE) { - .Call(`_rMVP_kin_cal_s`, pBigMat, threads, mkl, verbose) +BigRowMean <- function(pBigMat, marker_bycol = TRUE, step = 10000L, threads = 0L, geno_ind = NULL, verbose = TRUE) { + .Call(`_rMVP_BigRowMean`, pBigMat, marker_bycol, step, threads, geno_ind, verbose) } kin_cal <- function(pBigMat, geno_ind = NULL, marker_ind = NULL, marker_freq = NULL, marker_bycol = TRUE, threads = 0L, step = 10000L, mkl = FALSE, verbose = TRUE) { diff --git a/README.md b/README.md index 5ddb3fd..6b60cda 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@

-### #----------------------***rMVP [v1.1.0]() is coming, and stronger again!***------------------------# +### #----------------------***rMVP [v1.3.5]() is coming, and stronger and faster again!***------------------------# ### Repos: **Github:** https://github.com/xiaolei-lab/rMVP diff --git a/man/MVP.GLM.Rd b/man/MVP.GLM.Rd index f09baed..4cf3d05 100755 --- a/man/MVP.GLM.Rd +++ b/man/MVP.GLM.Rd @@ -10,6 +10,7 @@ MVP.GLM( CV = NULL, ind_idx = NULL, mrk_idx = NULL, + mrk_bycol = TRUE, maxLine = 5000, cpu = 1, verbose = TRUE @@ -26,6 +27,8 @@ MVP.GLM( \item{mrk_idx}{the index of effective markers used in analysis} +\item{mrk_bycol}{whether the markers are stored by columns in genotype (i.e. M is a n by m matrix)} + \item{maxLine}{the number of markers handled at a time, smaller value would reduce the memory cost} \item{cpu}{number of cpus used for parallel computation} diff --git a/man/MVP.MLM.Rd b/man/MVP.MLM.Rd index 5943e99..b588c08 100755 --- a/man/MVP.MLM.Rd +++ b/man/MVP.MLM.Rd @@ -12,6 +12,7 @@ MVP.MLM( CV = NULL, ind_idx = NULL, mrk_idx = NULL, + mrk_bycol = TRUE, REML = NULL, maxLine = 5000, cpu = 1, @@ -34,6 +35,8 @@ MVP.MLM( \item{mrk_idx}{the index of effective markers used in analysis} +\item{mrk_bycol}{whether the markers are stored by columns in genotype (i.e. M is a n by m matrix)} + \item{REML}{a list that contains ve and vg} \item{maxLine}{the number of markers handled at a time, smaller value would reduce the memory cost} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index a741d05..b19df08 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -25,8 +25,8 @@ BEGIN_RCPP END_RCPP } // glm_c -SEXP glm_c(const arma::vec& y, const arma::mat& X, const arma::mat& iXX, SEXP pBigMat, const Nullable geno_ind, const Nullable marker_ind, const int step, const bool verbose, const int threads); -RcppExport SEXP _rMVP_glm_c(SEXP ySEXP, SEXP XSEXP, SEXP iXXSEXP, SEXP pBigMatSEXP, SEXP geno_indSEXP, SEXP marker_indSEXP, SEXP stepSEXP, SEXP verboseSEXP, SEXP threadsSEXP) { +SEXP glm_c(const arma::vec& y, const arma::mat& X, const arma::mat& iXX, SEXP pBigMat, const Nullable geno_ind, const Nullable marker_ind, const bool marker_bycol, const int step, const bool verbose, const int threads); +RcppExport SEXP _rMVP_glm_c(SEXP ySEXP, SEXP XSEXP, SEXP iXXSEXP, SEXP pBigMatSEXP, SEXP geno_indSEXP, SEXP marker_indSEXP, SEXP marker_bycolSEXP, SEXP stepSEXP, SEXP verboseSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -36,16 +36,17 @@ BEGIN_RCPP Rcpp::traits::input_parameter< SEXP >::type pBigMat(pBigMatSEXP); Rcpp::traits::input_parameter< const Nullable >::type geno_ind(geno_indSEXP); Rcpp::traits::input_parameter< const Nullable >::type marker_ind(marker_indSEXP); + Rcpp::traits::input_parameter< const bool >::type marker_bycol(marker_bycolSEXP); Rcpp::traits::input_parameter< const int >::type step(stepSEXP); Rcpp::traits::input_parameter< const bool >::type verbose(verboseSEXP); Rcpp::traits::input_parameter< const int >::type threads(threadsSEXP); - rcpp_result_gen = Rcpp::wrap(glm_c(y, X, iXX, pBigMat, geno_ind, marker_ind, step, verbose, threads)); + rcpp_result_gen = Rcpp::wrap(glm_c(y, X, iXX, pBigMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads)); return rcpp_result_gen; END_RCPP } // mlm_c -SEXP mlm_c(const arma::vec& y, const arma::mat& X, const arma::mat& U, const double vgs, SEXP pBigMat, const Nullable geno_ind, const Nullable marker_ind, const int step, const bool verbose, const int threads); -RcppExport SEXP _rMVP_mlm_c(SEXP ySEXP, SEXP XSEXP, SEXP USEXP, SEXP vgsSEXP, SEXP pBigMatSEXP, SEXP geno_indSEXP, SEXP marker_indSEXP, SEXP stepSEXP, SEXP verboseSEXP, SEXP threadsSEXP) { +SEXP mlm_c(const arma::vec& y, const arma::mat& X, const arma::mat& U, const double vgs, SEXP pBigMat, const Nullable geno_ind, const Nullable marker_ind, const bool marker_bycol, const int step, const bool verbose, const int threads); +RcppExport SEXP _rMVP_mlm_c(SEXP ySEXP, SEXP XSEXP, SEXP USEXP, SEXP vgsSEXP, SEXP pBigMatSEXP, SEXP geno_indSEXP, SEXP marker_indSEXP, SEXP marker_bycolSEXP, SEXP stepSEXP, SEXP verboseSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -56,10 +57,11 @@ BEGIN_RCPP Rcpp::traits::input_parameter< SEXP >::type pBigMat(pBigMatSEXP); Rcpp::traits::input_parameter< const Nullable >::type geno_ind(geno_indSEXP); Rcpp::traits::input_parameter< const Nullable >::type marker_ind(marker_indSEXP); + Rcpp::traits::input_parameter< const bool >::type marker_bycol(marker_bycolSEXP); Rcpp::traits::input_parameter< const int >::type step(stepSEXP); Rcpp::traits::input_parameter< const bool >::type verbose(verboseSEXP); Rcpp::traits::input_parameter< const int >::type threads(threadsSEXP); - rcpp_result_gen = Rcpp::wrap(mlm_c(y, X, U, vgs, pBigMat, geno_ind, marker_ind, step, verbose, threads)); + rcpp_result_gen = Rcpp::wrap(mlm_c(y, X, U, vgs, pBigMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads)); return rcpp_result_gen; END_RCPP } @@ -225,43 +227,18 @@ BEGIN_RCPP END_RCPP } // BigRowMean -arma::vec BigRowMean(SEXP pBigMat, bool mrkbycol, int threads, const Nullable geno_ind); -RcppExport SEXP _rMVP_BigRowMean(SEXP pBigMatSEXP, SEXP mrkbycolSEXP, SEXP threadsSEXP, SEXP geno_indSEXP) { +arma::vec BigRowMean(SEXP pBigMat, bool marker_bycol, size_t step, int threads, const Nullable geno_ind, const bool verbose); +RcppExport SEXP _rMVP_BigRowMean(SEXP pBigMatSEXP, SEXP marker_bycolSEXP, SEXP stepSEXP, SEXP threadsSEXP, SEXP geno_indSEXP, SEXP verboseSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< SEXP >::type pBigMat(pBigMatSEXP); - Rcpp::traits::input_parameter< bool >::type mrkbycol(mrkbycolSEXP); + Rcpp::traits::input_parameter< bool >::type marker_bycol(marker_bycolSEXP); + Rcpp::traits::input_parameter< size_t >::type step(stepSEXP); Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); Rcpp::traits::input_parameter< const Nullable >::type geno_ind(geno_indSEXP); - rcpp_result_gen = Rcpp::wrap(BigRowMean(pBigMat, mrkbycol, threads, geno_ind)); - return rcpp_result_gen; -END_RCPP -} -// kin_cal_m -SEXP kin_cal_m(SEXP pBigMat, int threads, bool verbose); -RcppExport SEXP _rMVP_kin_cal_m(SEXP pBigMatSEXP, SEXP threadsSEXP, SEXP verboseSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< SEXP >::type pBigMat(pBigMatSEXP); - Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); - Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP); - rcpp_result_gen = Rcpp::wrap(kin_cal_m(pBigMat, threads, verbose)); - return rcpp_result_gen; -END_RCPP -} -// kin_cal_s -SEXP kin_cal_s(SEXP pBigMat, int threads, bool mkl, bool verbose); -RcppExport SEXP _rMVP_kin_cal_s(SEXP pBigMatSEXP, SEXP threadsSEXP, SEXP mklSEXP, SEXP verboseSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< SEXP >::type pBigMat(pBigMatSEXP); - Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); - Rcpp::traits::input_parameter< bool >::type mkl(mklSEXP); - Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP); - rcpp_result_gen = Rcpp::wrap(kin_cal_s(pBigMat, threads, mkl, verbose)); + Rcpp::traits::input_parameter< const bool >::type verbose(verboseSEXP); + rcpp_result_gen = Rcpp::wrap(BigRowMean(pBigMat, marker_bycol, step, threads, geno_ind, verbose)); return rcpp_result_gen; END_RCPP } @@ -287,8 +264,8 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_rMVP_getRow", (DL_FUNC) &_rMVP_getRow, 2}, - {"_rMVP_glm_c", (DL_FUNC) &_rMVP_glm_c, 9}, - {"_rMVP_mlm_c", (DL_FUNC) &_rMVP_mlm_c, 10}, + {"_rMVP_glm_c", (DL_FUNC) &_rMVP_glm_c, 10}, + {"_rMVP_mlm_c", (DL_FUNC) &_rMVP_mlm_c, 11}, {"_rMVP_vcf_parser_map", (DL_FUNC) &_rMVP_vcf_parser_map, 2}, {"_rMVP_vcf_parser_genotype", (DL_FUNC) &_rMVP_vcf_parser_genotype, 5}, {"_rMVP_hapmap_parser_map", (DL_FUNC) &_rMVP_hapmap_parser_map, 2}, @@ -301,9 +278,7 @@ static const R_CallMethodDef CallEntries[] = { {"_rMVP_geninv", (DL_FUNC) &_rMVP_geninv, 1}, {"_rMVP_impute_marker", (DL_FUNC) &_rMVP_impute_marker, 4}, {"_rMVP_hasNA", (DL_FUNC) &_rMVP_hasNA, 5}, - {"_rMVP_BigRowMean", (DL_FUNC) &_rMVP_BigRowMean, 4}, - {"_rMVP_kin_cal_m", (DL_FUNC) &_rMVP_kin_cal_m, 3}, - {"_rMVP_kin_cal_s", (DL_FUNC) &_rMVP_kin_cal_s, 4}, + {"_rMVP_BigRowMean", (DL_FUNC) &_rMVP_BigRowMean, 6}, {"_rMVP_kin_cal", (DL_FUNC) &_rMVP_kin_cal, 9}, {NULL, NULL, 0} }; diff --git a/src/assoc.cpp b/src/assoc.cpp index eb34db0..fc12a1c 100755 --- a/src/assoc.cpp +++ b/src/assoc.cpp @@ -67,12 +67,11 @@ NumericVector getRow(SEXP pBigMat, const int row){ } template -SEXP glm_c(const arma::vec &y, const arma::mat &X, const arma::mat & iXX, XPtr pMat, const Nullable geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const int step = 10000, const bool verbose = true, const int threads = 0){ +SEXP glm_c(const arma::vec &y, const arma::mat &X, const arma::mat & iXX, XPtr pMat, const Nullable geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const bool marker_bycol = true, const int step = 10000, const bool verbose = true, const int threads = 0){ omp_setup(threads); MatrixAccessor genomat = MatrixAccessor(*pMat); - bool marker_bycol = (y.n_elem == pMat->nrow()); int n; uvec _geno_ind; @@ -245,32 +244,31 @@ SEXP glm_c(const arma::vec &y, const arma::mat &X, const arma::mat & iXX, XPtr geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const int step = 10000, const bool verbose = true, const int threads = 0){ +SEXP glm_c(const arma::vec & y, const arma::mat & X, const arma::mat & iXX, SEXP pBigMat, const Nullable geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const bool marker_bycol = true, const int step = 10000, const bool verbose = true, const int threads = 0){ XPtr xpMat(pBigMat); switch(xpMat->matrix_type()){ case 1: - return glm_c(y, X, iXX, xpMat, geno_ind, marker_ind, step, verbose, threads); + return glm_c(y, X, iXX, xpMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads); case 2: - return glm_c(y, X, iXX, xpMat, geno_ind, marker_ind, step, verbose, threads); + return glm_c(y, X, iXX, xpMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads); case 4: - return glm_c(y, X, iXX, xpMat, geno_ind, marker_ind, step, verbose, threads); + return glm_c(y, X, iXX, xpMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads); case 8: - return glm_c(y, X, iXX, xpMat, geno_ind, marker_ind, step, verbose, threads); + return glm_c(y, X, iXX, xpMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads); default: throw Rcpp::exception("unknown type detected for big.matrix object!"); } } template -SEXP mlm_c(const arma::vec & y, const arma::mat & X, const arma::mat & U, const double vgs, XPtr pMat, const Nullable geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const int step = 10000, const bool verbose = true, const int threads = 0){ +SEXP mlm_c(const arma::vec & y, const arma::mat & X, const arma::mat & U, const double vgs, XPtr pMat, const Nullable geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const bool marker_bycol = true, const int step = 10000, const bool verbose = true, const int threads = 0){ omp_setup(threads); MatrixAccessor genomat = MatrixAccessor(*pMat); - bool marker_bycol = (y.n_elem == pMat->nrow()); - + int n; uvec _geno_ind; if(geno_ind.isNotNull()){ @@ -303,7 +301,6 @@ SEXP mlm_c(const arma::vec & y, const arma::mat & X, const arma::mat & U, const arma::mat res(m, 3); arma::mat iXXs(q0 + 1, q0 + 1); - arma::mat Z_buffer(n, step, fill::none); int i = 0, j = 0; int i_marker = 0; @@ -426,19 +423,19 @@ SEXP mlm_c(const arma::vec & y, const arma::mat & X, const arma::mat & U, const } // [[Rcpp::export]] -SEXP mlm_c(const arma::vec & y, const arma::mat & X, const arma::mat & U, const double vgs, SEXP pBigMat, const Nullable geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const int step = 10000, const bool verbose = true, const int threads = 0){ +SEXP mlm_c(const arma::vec & y, const arma::mat & X, const arma::mat & U, const double vgs, SEXP pBigMat, const Nullable geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const bool marker_bycol = true, const int step = 10000, const bool verbose = true, const int threads = 0){ XPtr xpMat(pBigMat); switch(xpMat->matrix_type()){ case 1: - return mlm_c(y, X, U, vgs, xpMat, geno_ind, marker_ind, step, verbose, threads); + return mlm_c(y, X, U, vgs, xpMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads); case 2: - return mlm_c(y, X, U, vgs, xpMat, geno_ind, marker_ind, step, verbose, threads); + return mlm_c(y, X, U, vgs, xpMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads); case 4: - return mlm_c(y, X, U, vgs, xpMat, geno_ind, marker_ind, step, verbose, threads); + return mlm_c(y, X, U, vgs, xpMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads); case 8: - return mlm_c(y, X, U, vgs, xpMat, geno_ind, marker_ind, step, verbose, threads); + return mlm_c(y, X, U, vgs, xpMat, geno_ind, marker_ind, marker_bycol, step, verbose, threads); default: throw Rcpp::exception("unknown type detected for big.matrix object!"); } diff --git a/src/kinship.cpp b/src/kinship.cpp index f993274..e375665 100755 --- a/src/kinship.cpp +++ b/src/kinship.cpp @@ -14,13 +14,13 @@ using namespace Rcpp; using namespace arma; template -arma::vec BigRowMean(XPtr pMat, bool mrkbycol = true, int threads = 0, const Nullable geno_ind = R_NilValue){ +arma::vec BigRowMean(XPtr pMat, bool marker_bycol = true, size_t step = 10000, int threads = 0, const Nullable geno_ind = R_NilValue, const bool verbose = true){ omp_setup(threads); MatrixAccessor bigm = MatrixAccessor(*pMat); int n; - int m = mrkbycol ? pMat->ncol() : pMat->nrow(); + int m = marker_bycol ? pMat->ncol() : pMat->nrow(); arma::vec mean(m); uvec _geno_ind; @@ -28,225 +28,288 @@ arma::vec BigRowMean(XPtr pMat, bool mrkbycol = true, int threads = 0 _geno_ind = as(geno_ind) - 1; n = _geno_ind.n_elem; }else{ - n = mrkbycol ? pMat->nrow() : pMat->ncol(); + n = marker_bycol ? pMat->nrow() : pMat->ncol(); } - if(mrkbycol){ - if(_geno_ind.is_empty()){ - #pragma omp parallel for - for (int j = 0; j < m; j++){ - double p1 = 0.0; - for(int k = 0; k < n; k++){ - p1 += bigm[j][k]; - } - mean[j] = p1 / n; - } - }else{ - #pragma omp parallel for - for (int j = 0; j < m; j++){ - double p1 = 0.0; - for(int k = 0; k < n; k++){ - p1 += bigm[j][_geno_ind[k]]; - } - mean[j] = p1 / n; - } + MinimalProgressBar_plus pb; + Progress progress(m, verbose, pb); + + arma::mat Z_buffer(n, step, fill::none); + int i = 0, j = 0; + int i_marker = 0; + for (;i < m;) { + + int cnt = 0; + for (; j < m && cnt < step; j++) + { + cnt++; } - }else{ + + if (cnt != step) { + Z_buffer.set_size(n, cnt); + } + if(_geno_ind.is_empty()){ - #pragma omp parallel for - for (int j = 0; j < m; j++){ - double p1 = 0.0; + if(marker_bycol){ + #pragma omp parallel for + for(int l = 0; l < cnt; l++){ + for(int k = 0; k < n; k++){ + Z_buffer(k, l) = bigm[(i_marker + l)][k]; + } + } + }else{ + #pragma omp parallel for for(int k = 0; k < n; k++){ - p1 += bigm[k][j]; + for(int l = 0; l < cnt; l++){ + Z_buffer(k, l) = bigm[k][(i_marker + l)]; + } } - mean[j] = p1 / n; } }else{ - #pragma omp parallel for - for (int j = 0; j < m; j++){ - double p1 = 0.0; + if(marker_bycol){ + #pragma omp parallel for + for(int l = 0; l < cnt; l++){ + for(int k = 0; k < n; k++){ + Z_buffer(k, l) = bigm[(i_marker + l)][_geno_ind[k]]; + } + } + }else{ + #pragma omp parallel for for(int k = 0; k < n; k++){ - p1 += bigm[_geno_ind[k]][j]; + for(int l = 0; l < cnt; l++){ + Z_buffer(k, l) = bigm[_geno_ind[k]][(i_marker + l)]; + } } - mean[j] = p1 / n; } } - } - return mean; -} - -// [[Rcpp::export]] -arma::vec BigRowMean(SEXP pBigMat, bool mrkbycol = true, int threads = 0, const Nullable geno_ind = R_NilValue){ - - XPtr xpMat(pBigMat); + #pragma omp parallel for + for(int l = 0; l < cnt; l++){ + mean[i_marker + l] = arma::mean(Z_buffer.col(l)); + } - switch(xpMat->matrix_type()) { - case 1: - return BigRowMean(xpMat, mrkbycol, threads, geno_ind); - case 2: - return BigRowMean(xpMat, mrkbycol, threads, geno_ind); - case 4: - return BigRowMean(xpMat, mrkbycol, threads, geno_ind); - case 8: - return BigRowMean(xpMat, mrkbycol, threads, geno_ind); - default: - throw Rcpp::exception("unknown type detected for big.matrix object!"); + i = j; + i_marker += cnt; + if(!Progress::check_abort()) progress.increment(cnt); } -} - - -template -SEXP kin_cal_m(XPtr pMat, int threads = 0, bool verbose = true){ - - omp_setup(threads); - - if(verbose) - Rcout << "Computing GRM under mode: Memory" << endl; - - MatrixAccessor bigm = MatrixAccessor(*pMat); - - int n = pMat->ncol(); - int m = pMat->nrow(); - int i = 0, j = 0, k = 0; - MinimalProgressBar pb; - - arma::vec Mean = BigRowMean(pMat, threads); - double SUM = sum((0.5 * Mean) % (1 - 0.5 * Mean)); - - arma::mat kin(n, n); - arma::vec coli(m); - arma::vec colj(m); - - Progress p(n, verbose, pb); - - if(verbose) - Rcout << "Scale the genotype matrix and compute Z'Z" << endl; + Z_buffer.reset(); - #pragma omp parallel for schedule(dynamic) firstprivate(coli, colj) private(i, j, k) - for(i = 0; i < n; i++){ - for(k = 0; k < m; k++){ - coli[k] = bigm[i][k] - Mean[k]; - } - if ( ! Progress::check_abort() ) { - p.increment(); - for(j = i; j < n; j++){ - for(k = 0; k < m; k++){ - colj[k] = bigm[j][k] - Mean[k]; - } - kin(i, j) = kin(j, i) = 0.5 * sum(coli % colj) / SUM; - } - } - } + // if(marker_bycol){ + // if(_geno_ind.is_empty()){ + // #pragma omp parallel for + // for (int j = 0; j < m; j++){ + // double p1 = 0.0; + // for(int k = 0; k < n; k++){ + // p1 += bigm[j][k]; + // } + // mean[j] = p1 / n; + // } + // }else{ + // #pragma omp parallel for + // for (int j = 0; j < m; j++){ + // double p1 = 0.0; + // for(int k = 0; k < n; k++){ + // p1 += bigm[j][_geno_ind[k]]; + // } + // mean[j] = p1 / n; + // } + // } + // }else{ + // if(_geno_ind.is_empty()){ + // #pragma omp parallel for + // for (int j = 0; j < m; j++){ + // double p1 = 0.0; + // for(int k = 0; k < n; k++){ + // p1 += bigm[k][j]; + // } + // mean[j] = p1 / n; + // } + // }else{ + // #pragma omp parallel for + // for (int j = 0; j < m; j++){ + // double p1 = 0.0; + // for(int k = 0; k < n; k++){ + // p1 += bigm[_geno_ind[k]][j]; + // } + // mean[j] = p1 / n; + // } + // } + // } - return Rcpp::wrap(kin); + return mean; } - // [[Rcpp::export]] -SEXP kin_cal_m(SEXP pBigMat, int threads = 0, bool verbose = true){ - +arma::vec BigRowMean(SEXP pBigMat, bool marker_bycol = true, size_t step = 10000, int threads = 0, const Nullable geno_ind = R_NilValue, const bool verbose = true){ + XPtr xpMat(pBigMat); - switch(xpMat->matrix_type()){ + switch(xpMat->matrix_type()) { case 1: - return kin_cal_m(xpMat, threads, verbose); + return BigRowMean(xpMat, marker_bycol, step, threads, geno_ind, verbose); case 2: - return kin_cal_m(xpMat, threads, verbose); + return BigRowMean(xpMat, marker_bycol, step, threads, geno_ind, verbose); case 4: - return kin_cal_m(xpMat, threads, verbose); + return BigRowMean(xpMat, marker_bycol, step, threads, geno_ind, verbose); case 8: - return kin_cal_m(xpMat, threads, verbose); + return BigRowMean(xpMat, marker_bycol, step, threads, geno_ind, verbose); default: throw Rcpp::exception("unknown type detected for big.matrix object!"); } } -template -SEXP kin_cal_s(XPtr pMat, int threads = 0, bool mkl = false, bool verbose = true){ - - omp_setup(threads); - - if(verbose) - Rcout << "Computing GRM under mode: Speed" << endl; - - MatrixAccessor bigm = MatrixAccessor(*pMat); - - #ifdef _OPENMP - #else - if(!mkl) - mkl = true; - #endif - if(threads == 1) - mkl = true; - - int n = pMat->ncol(); - int m = pMat->nrow(); - int i = 0, j = 0; - MinimalProgressBar pb; +// template +// SEXP kin_cal_m(XPtr pMat, int threads = 0, bool verbose = true){ - arma::vec Mean = BigRowMean(pMat, threads); - double SUM = sum((0.5 * Mean) % (1 - 0.5 * Mean)); +// omp_setup(threads); - arma::mat kin(n, n); - arma::mat geno(m, n); +// if(verbose) +// Rcout << "Computing GRM under mode: Memory" << endl; - if(verbose) - Rcout << "Scale the genotype matrix" << endl; +// MatrixAccessor bigm = MatrixAccessor(*pMat); - #pragma omp parallel for schedule(dynamic) private(i, j) - for(i = 0; i < n; i++){ - for(j = 0; j < m; j++){ - geno(j, i) = bigm[i][j] - Mean[j]; - } - } +// int n = pMat->ncol(); +// int m = pMat->nrow(); +// int i = 0, j = 0, k = 0; +// MinimalProgressBar pb; - if(verbose) - Rcout << "Computing Z'Z" << endl; - - if(mkl){ - kin = geno.t() * geno / SUM / 2; - }else{ +// arma::vec Mean = BigRowMean(pMat, threads); +// double SUM = sum((0.5 * Mean) % (1 - 0.5 * Mean)); - Progress p(n, verbose, pb); - arma::colvec coli; - - #pragma omp parallel for schedule(dynamic) private(i, j, coli) - for(i = 0; i < n; i++){ - coli = geno.col(i); - if ( ! Progress::check_abort() ) { - p.increment(); - for(j = i; j < n; j++){ - kin(j, i) = kin(i, j) = 0.5 * sum(coli % geno.col(j)) / SUM; - } - } - } +// arma::mat kin(n, n); +// arma::vec coli(m); +// arma::vec colj(m); - } +// Progress p(n, verbose, pb); - return Rcpp::wrap(kin); -} +// if(verbose) +// Rcout << "Scale the genotype matrix and compute Z'Z" << endl; +// #pragma omp parallel for schedule(dynamic) firstprivate(coli, colj) private(i, j, k) +// for(i = 0; i < n; i++){ +// for(k = 0; k < m; k++){ +// coli[k] = bigm[i][k] - Mean[k]; +// } +// if ( ! Progress::check_abort() ) { +// p.increment(); +// for(j = i; j < n; j++){ +// for(k = 0; k < m; k++){ +// colj[k] = bigm[j][k] - Mean[k]; +// } +// kin(i, j) = kin(j, i) = 0.5 * sum(coli % colj) / SUM; +// } +// } +// } -// [[Rcpp::export]] -SEXP kin_cal_s(SEXP pBigMat, int threads = 0, bool mkl = false, bool verbose = true){ +// return Rcpp::wrap(kin); +// } - XPtr xpMat(pBigMat); - switch(xpMat->matrix_type()){ - case 1: - return kin_cal_s(xpMat, threads, mkl, verbose); - case 2: - return kin_cal_s(xpMat, threads, mkl, verbose); - case 4: - return kin_cal_s(xpMat, threads, mkl, verbose); - case 8: - return kin_cal_s(xpMat, threads, mkl, verbose); - default: - throw Rcpp::exception("unknown type detected for big.matrix object!"); - } -} +// // [[Rcpp::export]] +// SEXP kin_cal_m(SEXP pBigMat, int threads = 0, bool verbose = true){ + +// XPtr xpMat(pBigMat); + +// switch(xpMat->matrix_type()){ +// case 1: +// return kin_cal_m(xpMat, threads, verbose); +// case 2: +// return kin_cal_m(xpMat, threads, verbose); +// case 4: +// return kin_cal_m(xpMat, threads, verbose); +// case 8: +// return kin_cal_m(xpMat, threads, verbose); +// default: +// throw Rcpp::exception("unknown type detected for big.matrix object!"); +// } +// } + + +// template +// SEXP kin_cal_s(XPtr pMat, int threads = 0, bool mkl = false, bool verbose = true){ + +// omp_setup(threads); + +// if(verbose) +// Rcout << "Computing GRM under mode: Speed" << endl; + +// MatrixAccessor bigm = MatrixAccessor(*pMat); + +// #ifdef _OPENMP +// #else +// if(!mkl) +// mkl = true; +// #endif +// if(threads == 1) +// mkl = true; + +// int n = pMat->ncol(); +// int m = pMat->nrow(); +// int i = 0, j = 0; +// MinimalProgressBar pb; + +// arma::vec Mean = BigRowMean(pMat, threads); +// double SUM = sum((0.5 * Mean) % (1 - 0.5 * Mean)); + +// arma::mat kin(n, n); +// arma::mat geno(m, n); + +// if(verbose) +// Rcout << "Scale the genotype matrix" << endl; + +// #pragma omp parallel for schedule(dynamic) private(i, j) +// for(i = 0; i < n; i++){ +// for(j = 0; j < m; j++){ +// geno(j, i) = bigm[i][j] - Mean[j]; +// } +// } + +// if(verbose) +// Rcout << "Computing Z'Z" << endl; + +// if(mkl){ +// kin = geno.t() * geno / SUM / 2; +// }else{ + +// Progress p(n, verbose, pb); +// arma::colvec coli; + +// #pragma omp parallel for schedule(dynamic) private(i, j, coli) +// for(i = 0; i < n; i++){ +// coli = geno.col(i); +// if ( ! Progress::check_abort() ) { +// p.increment(); +// for(j = i; j < n; j++){ +// kin(j, i) = kin(i, j) = 0.5 * sum(coli % geno.col(j)) / SUM; +// } +// } +// } + +// } + +// return Rcpp::wrap(kin); +// } + + +// // [[Rcpp::export]] +// SEXP kin_cal_s(SEXP pBigMat, int threads = 0, bool mkl = false, bool verbose = true){ + +// XPtr xpMat(pBigMat); + +// switch(xpMat->matrix_type()){ +// case 1: +// return kin_cal_s(xpMat, threads, mkl, verbose); +// case 2: +// return kin_cal_s(xpMat, threads, mkl, verbose); +// case 4: +// return kin_cal_s(xpMat, threads, mkl, verbose); +// case 8: +// return kin_cal_s(xpMat, threads, mkl, verbose); +// default: +// throw Rcpp::exception("unknown type detected for big.matrix object!"); +// } +// } template