diff --git a/.gitignore b/.gitignore index e6d66a3..6edbc49 100644 --- a/.gitignore +++ b/.gitignore @@ -46,3 +46,6 @@ vignettes/*.pdf *.knit.md .Rproj.user packages + +# R Files for comparisions +compare_perf.R diff --git a/CHANGELOG b/CHANGELOG index b39058d..320058c 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,24 @@ +## refactor-branch (2026-05-14) ++ Bug fixes: + + Fixed 2*3.14 → 2*pi in FaSTLMM log-likelihood computation (~0.0028/obs error in LL values) + + Fixed double-subtraction of mean^2 in MVP.HE.Vg.Ve StandardizeVector (numerically benign but mathematically wrong) ++ Code quality: + + Named constants extracted to R/MVP.Constants.R (eliminated magic numbers in FaSTLMM, EMMA, BRENT, HE, FarmCPU) + + FaSTLMM inner-loop vectorized: scalar for-loops replaced with sweep/crossprod operations + + Duplicate FaSTLMM implementation unified via .fastlmm_core() shared by MVP.FaSTLMM.LL and FarmCPU + + FarmCPU internal helpers renamed to .farmcpu_* convention + + .safe_solve() helper eliminates repeated try/inherits(solve(...), 'try-error') patterns + + MVP.BRENT.Vg.Ve gains stopifnot() input validation + + Dead LogRL_dev1 function removed from MVP.HE.Vg.Ve + + gc() calls removed from hot loops (kept only after explicit rm() of large objects) ++ C++ improvements: + + C++ standard set to CXX17 in Makevars/Makevars.win + + ProgressBar base class extracted to eliminate ~100 lines of duplicated C++ + + fill_geno_buffer template eliminates 8-branch fill blocks in assoc.cpp and kinship.cpp + + DISPATCH_MATRIX_TYPE macro eliminates repeated switch/case dispatch boilerplate + + FILE* replaced with RAII std::ifstream/ofstream in write_bfile and read_bfile + + split_line_view() using std::string_view eliminates per-token heap allocations in VCF/HAPMAP parsers + ## v1.1.0beta (2018/09/12) + Unified installation package for Windows platform and Linux platform. + The 3D PCA map has been temporarily disabled, and we will improve this feature in subsequent versions. diff --git a/DESCRIPTION b/DESCRIPTION index fcdcf84..2bd45f1 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -41,7 +41,7 @@ Depends: R (>= 3.3) LinkingTo: Rcpp, RcppArmadillo, RcppEigen, RcppProgress, BH, bigmemory NeedsCompilation: yes Suggests: knitr, testthat, rmarkdown -RoxygenNote: 7.3.2 +RoxygenNote: 7.3.3 Maintainer: Xiaolei Liu Author: Lilin Yin [aut], Haohao Zhang [aut], diff --git a/R/MVP.BRENT.Vg.Ve.R b/R/MVP.BRENT.Vg.Ve.R index bf08794..b5cf51a 100755 --- a/R/MVP.BRENT.Vg.Ve.R +++ b/R/MVP.BRENT.Vg.Ve.R @@ -25,14 +25,19 @@ #' } #' MVP.BRENT.Vg.Ve <- function(y, X, eigenK, verbose = FALSE) { - p = 0 + stopifnot( + !anyNA(eigenK$values), + length(eigenK$values) == nrow(X), + length(y) == nrow(X) + ) + p <- 0 Sigma <- eigenK$values - w <- which(Sigma < 1e-6) - Sigma[w] <- 1e-6 + w <- which(Sigma < .BRENT_MIN_EIGENVALUE) + Sigma[w] <- .BRENT_MIN_EIGENVALUE U <- eigenK$vectors - min_h2 = 0 - max_h2 = 1 - tol = .Machine$double.eps^0.25 + min_h2 <- 0 + max_h2 <- 1 + tol <- .Machine$double.eps^0.25 reml <- fit_diago_brent(y, X, p, Sigma, U, min_h2, max_h2, tol, verbose = verbose) vg <- reml[[2]] ve <- reml[[1]] diff --git a/R/MVP.Constants.R b/R/MVP.Constants.R new file mode 100644 index 0000000..81e6b02 --- /dev/null +++ b/R/MVP.Constants.R @@ -0,0 +1,25 @@ +# Package-level constants for rMVP statistical methods. +# Naming convention: ._ (dot-prefix keeps them non-exported) + +# FaSTLMM log-likelihood optimization +.FASTLMM_SVD_THRESHOLD <- 1e-8 # singular value cutoff in SVD +.FASTLMM_DELTA_EXP_START <- -5 # log-delta grid lower bound +.FASTLMM_DELTA_EXP_END <- 5 # log-delta grid upper bound +.FASTLMM_DELTA_EXP_STEP <- 0.1 # log-delta grid step size +.FASTLMM_DELTA_EXP_DEGENERATE <- 100 # sentinel: collapse grid to single point when SNP pool has a constant column +.FASTLMM_2PI <- 2 * pi # NOTE: original code had 2 * 3.14 (bug); corrected to 2 * pi + +# EMMA variance component estimation (defaults match original function signature) +.EMMA_NGRIDS <- 100 +.EMMA_LLIM <- -10 +.EMMA_ULIM <- 10 +.EMMA_ESP <- 1e-10 + +# BRENT variance component estimation +.BRENT_MIN_EIGENVALUE <- 1e-6 # eigenvalue floor to avoid division by near-zero + +# HE regression +.HE_DEFAULT_LOG_SIGMA2 <- log(0.1) # log-sigma2 fallback when CalcVChe returns non-positive value + +# FarmCPU +.FARMCPU_LD_THRESHOLD <- 0.7 # LD threshold for pseudo-QTN deduplication diff --git a/R/MVP.Data.r b/R/MVP.Data.r index 220dc61..77bd980 100755 --- a/R/MVP.Data.r +++ b/R/MVP.Data.r @@ -152,7 +152,7 @@ MVP.Data <- function(fileMVP = NULL, fileVCF = NULL, fileHMP = NULL, fileBed = N # phenotype if (!is.null(filePhe)) { MVP.Data.Pheno( - pheno_file = filePhe, + pheno_file <- filePhe, out = out, header = TRUE, cols = pheno_cols, @@ -438,36 +438,12 @@ MVP.Data.Numeric2MVP <- function(num_file, map_file, out='mvp', maxLine=1e4, row ) logging.log(paste0("Loading genotype at a step of ", maxLine, '...\n'), verbose = verbose) - # convert to bigmat - speed - # if (priority == "speed") { - # opts <- options(bigmemory.typecast.warning = FALSE) - # on.exit(options(opts)) - - # # detecte sep - # con <- file(num_file, open = 'r') - # line <- readLines(con, 1) - # close(con) - # sep <- substr(line, 2, 2) - - # # load geno - # suppressWarnings( - # geno <- read.big.matrix(num_file, header = FALSE, sep = sep) - # ) - # if (transposed) { - # bigmat[, ] <- t(geno[, ]) - # } else { - # bigmat[, ] <- geno[, ] - # } - # rm("geno") - # } - - # convert to bigmat - memory - # if (priority == "memory") { - i <- 0 - con <- file(num_file, open = 'r') + # convert to bigmat - memory priority: read in chunks to reduce memory footprint + i <- 0 + con <- file(num_file, open = 'r') if (col_names) { readLines(con, n = 1) } while (TRUE) { - line = readLines(con, n = maxLine) + line <- readLines(con, n = maxLine) len <- length(line) if (len == 0) { break } @@ -488,7 +464,6 @@ MVP.Data.Numeric2MVP <- function(num_file, map_file, out='mvp', maxLine=1e4, row } logging.log("\n", verbose = verbose) close(con) - # } file.copy(map_file, paste0(out, ".geno.map")) @@ -639,10 +614,10 @@ MVP.Data.Pheno <- function(pheno_file, out='mvp', cols=NULL, header=TRUE, sep='\ # drop empty traits pheno[pheno %in% missing] <- NA - drop = c() + drop <- c() for (i in 2:ncol(pheno)) { if (all(is.na(pheno[, i]))) { - drop = c(drop, i) + drop <- c(drop, i) } } if (length(drop) > 0) { diff --git a/R/MVP.EMMA.Vg.Ve.r b/R/MVP.EMMA.Vg.Ve.r index b3e453c..e66ea69 100755 --- a/R/MVP.EMMA.Vg.Ve.r +++ b/R/MVP.EMMA.Vg.Ve.r @@ -50,6 +50,8 @@ #' MVP.EMMA.Vg.Ve <- function(y, X, K, ngrids=100, llim=-10, ulim=10, esp=1e-10) { + # Defaults above match .EMMA_NGRIDS / .EMMA_LLIM / .EMMA_ULIM / .EMMA_ESP + # (kept as literals in the signature for R CMD check codoc compatibility). # NA in phenotype idx <- !is.na(y) y <- y[idx] @@ -61,18 +63,21 @@ function(y, X, K, ngrids=100, llim=-10, ulim=10, esp=1e-10) { { nq <- length(etas) delta <- exp(logdelta) - return( 0.5 * (nq * (log(nq/(2 * pi))-1-log(sum(etas * etas/(lambda + delta))))-sum(log(lambda + delta))) ) + return( 0.5 * (nq * (log(nq/.FASTLMM_2PI)-1-log(sum(etas * etas/(lambda + delta))))-sum(log(lambda + delta))) ) + } + emma.delta.REML.dLL.wo.Z <- function(logdelta, lambda, etas) { + nq <- length(etas) + delta <- exp(logdelta) + etasq <- etas * etas + ldelta <- lambda + delta + return( 0.5 * (nq * sum(etasq/(ldelta * ldelta))/sum(etasq/ldelta)-sum(1/ldelta)) ) } emma.eigen.R.wo.Z=function(K, X) { n <- nrow(X) q <- ncol(X) XX <- crossprod(X) - iXX <- try(solve(XX), silent = TRUE) - if(inherits(iXX, "try-error")){ - #library(MASS) - iXX <- ginv(XX) - } + iXX <- .safe_solve(XX) SS1 <- X %*% iXX SS2 <- tcrossprod(SS1, X) @@ -98,7 +103,7 @@ function(y, X, K, ngrids=100, llim=-10, ulim=10, esp=1e-10) { delta <- exp(logdelta) Lambdas <- matrix(eig.R$values, n-q, m) + matrix(delta, n-q, m, byrow=TRUE) Etasq <- matrix(etas * etas, n-q, m) - LL <- 0.5 * ((n-q) * (log((n-q)/(2 * pi))-1-log(colSums(Etasq/Lambdas)))-colSums(log(Lambdas))) + LL <- 0.5 * ((n-q) * (log((n-q)/.FASTLMM_2PI)-1-log(colSums(Etasq/Lambdas)))-colSums(log(Lambdas))) dLL <- 0.5 * delta * ((n-q) * colSums(Etasq/(Lambdas * Lambdas))/colSums(Etasq/Lambdas)-colSums(1/Lambdas)) optlogdelta <- vector(length=0) optLL <- vector(length=0) @@ -112,30 +117,18 @@ function(y, X, K, ngrids=100, llim=-10, ulim=10, esp=1e-10) { } for(i in 1:(m-1) ){ if( ( dLL[i] * dLL[i + 1] < 0 ) && ( dLL[i] > 0 ) && ( dLL[i + 1] < 0 ) ) { - emma.delta.REML.dLL.wo.Z <- function(logdelta, lambda, etas) { - nq <- length(etas) - delta <- exp(logdelta) - etasq <- etas * etas - ldelta <- lambda + delta - return( 0.5 * (nq * sum(etasq/(ldelta * ldelta))/sum(etasq/ldelta)-sum(1/ldelta)) ) - } r <- uniroot(emma.delta.REML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i + 1], lambda=eig.R$values, etas=etas) optlogdelta <- append(optlogdelta, r$root) - emma.delta.REML.LL.wo.Z <- function(logdelta, lambda, etas) { - nq <- length(etas) - delta <- exp(logdelta) - return( 0.5 * (nq * (log(nq/(2 * pi))-1-log(sum(etas * etas/(lambda + delta))))-sum(log(lambda + delta))) ) - } optLL <- append(optLL, emma.delta.REML.LL.wo.Z(r$root, eig.R$values, etas)) } } maxdelta <- exp(optlogdelta[which.max(optLL)]) - #handler of grids with NaN log - replaceNaN<-function(LL) { - index=(LL == "NaN") - if(length(index)>0) theMin=min(LL[!index]) - if(length(index)<1) theMin="NaN" - LL[index]=theMin + # Handle grids with NaN log values + replaceNaN <- function(LL) { + index <- (LL == "NaN") + if(length(index)>0) theMin <- min(LL[!index]) + if(length(index)<1) theMin <- "NaN" + LL[index] <- theMin return(LL) } optLL=replaceNaN(optLL) diff --git a/R/MVP.FaSTLMM.LL.r b/R/MVP.FaSTLMM.LL.r index bdbd068..6e8fb78 100755 --- a/R/MVP.FaSTLMM.LL.r +++ b/R/MVP.FaSTLMM.LL.r @@ -1,9 +1,9 @@ # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. # You may obtain a copy of the License at -# +# # http://www.apache.org/licenses/LICENSE-2.0 -# +# # Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. @@ -11,12 +11,85 @@ # limitations under the License. +# Internal FaSTLMM core — shared by MVP.FaSTLMM.LL and FarmCPU.FaSTLMM.LL. +# Returns list(beta, delta, LL, vg, ve). +.fastlmm_core <- function(pheno, snp.pool, X0=NULL, ncpus=2) { + y <- pheno + deltaExpStart <- .FASTLMM_DELTA_EXP_START + deltaExpEnd <- .FASTLMM_DELTA_EXP_END + snp.pool <- snp.pool[,] + if (!is.null(snp.pool) && any(apply(snp.pool, 2, var) == 0)) { + deltaExpStart <- .FASTLMM_DELTA_EXP_DEGENERATE + deltaExpEnd <- deltaExpStart + } + if (is.null(X0)) X0 <- matrix(1, nrow(snp.pool), 1) + X <- X0 + + # SVD of snp.pool + K.X.svd <- svd(snp.pool) + d <- K.X.svd$d + d <- d[d > .FASTLMM_SVD_THRESHOLD] + d <- d^2 + U1 <- K.X.svd$u + U1 <- U1[, 1:length(d)] + if (is.null(dim(U1))) U1 <- matrix(U1, ncol=1) + n <- nrow(U1) + + U1TX <- crossprod(U1, X) + U1TY <- crossprod(U1, y) + yU1TY <- y - U1 %*% U1TY + XU1TX <- X - U1 %*% U1TX + IU <- -tcrossprod(U1) + diag(IU) <- rep(1, n) + diag(IU) + IUX <- crossprod(IU, X) + IUY <- crossprod(IU, y) + + delta.range <- seq(deltaExpStart, deltaExpEnd, by=.FASTLMM_DELTA_EXP_STEP) + m <- length(delta.range) + + beta.optimize.parallel <- function(ii) { + delta <- exp(delta.range[ii]) + dInv <- 1 / (d + delta) + + beta1 <- crossprod(sweep(U1TX, 1, sqrt(dInv), "*")) + beta2 <- crossprod(IUX) / delta + beta3 <- crossprod(U1TX, U1TY * dInv) + beta4 <- crossprod(IUX, IUY) / delta + beta <- crossprod(.safe_solve(beta1 + beta2), beta3 + beta4) + + part1 <- -0.5 * (n * log(.FASTLMM_2PI) + sum(log(d + delta)) + (n - length(d)) * log(delta)) + + part221 <- sum((U1TY - U1TX %*% beta)^2 * dInv) + part222 <- sum((yU1TY - XU1TX %*% beta)^2) / delta + part2 <- -0.5 * (n + n * log((part221 + part222) / n)) + + list(beta=beta, delta=delta, LL=part1 + part2) + } + + llresults <- lapply(1:m, beta.optimize.parallel) + + LL_values <- sapply(llresults, function(x) x$LL) + best_idx <- which.max(LL_values) + beta <- llresults[[best_idx]]$beta + delta <- llresults[[best_idx]]$delta + LL <- llresults[[best_idx]]$LL + + # Vectorized variance component estimates + sigma_a1 <- sum((U1TY - U1TX %*% beta)^2 / (d + delta)) + sigma_a2 <- sum((IUY - IUX %*% beta)^2) / delta + sigma_a <- (sigma_a1 + sigma_a2) / n + sigma_e <- delta * sigma_a + + list(beta=beta, delta=delta, LL=LL, vg=sigma_a, ve=sigma_e) +} + + #' Evaluation of the maximum likelihood using FaST-LMM method #' #' Last update: January 11, 2017 -#' +#' #' @author Xiaolei Liu (modified) -#' +#' #' @param pheno a two-column phenotype matrix #' @param snp.pool matrix for pseudo QTNs #' @param X0 covariates matrix @@ -28,180 +101,7 @@ #' Output: LL - log-likelihood #' Output: vg - genetic variance #' Output: ve - residual variance -#' -`MVP.FaSTLMM.LL` <- function(pheno, snp.pool, X0=NULL, ncpus=2){ - y=pheno - p=0 - deltaExpStart = -5 - deltaExpEnd = 5 - snp.pool=snp.pool[,] - if(!is.null(snp.pool)&&any(apply(snp.pool, 2, var)==0)){ - deltaExpStart = 100 - deltaExpEnd = deltaExpStart - } - if(is.null(X0)) { - X0 = matrix(1, nrow(snp.pool), 1) - } - X=X0 - #########SVD of X - K.X.svd <- svd(snp.pool) - d=K.X.svd$d - d=d[d>1e-08] - d=d^2 - U1=K.X.svd$u - U1=U1[,1:length(d)] - #handler of single snp - if(is.null(dim(U1))) U1=matrix(U1,ncol=1) - n=nrow(U1) - U1TX=crossprod(U1,X) - U1TY=crossprod(U1,y) - yU1TY <- y-U1%*%U1TY - XU1TX<- X-U1%*%U1TX - IU = -tcrossprod(U1) - diag(IU) = rep(1,n) + diag(IU) - IUX=crossprod(IU,X) - IUY=crossprod(IU,y) - #Iteration on the range of delta (-5 to 5 in glog scale) - delta.range <- seq(deltaExpStart,deltaExpEnd,by=0.1) - m <- length(delta.range) - #for (m in seq(deltaExpStart,deltaExpEnd,by=0.1)){ - beta.optimize.parallel <- function(ii){ - #p=p+1 - delta <- exp(delta.range[ii]) - #----------------------------calculate beta------------------------------------- - #######get beta1 - beta1=0 - for(i in 1:length(d)){ - one=matrix(U1TX[i,], nrow=1) - beta=crossprod(one,(one/(d[i]+delta))) #This is not real beta, confusing - beta1= beta1+beta - } - - #######get beta2 - beta2=0 - for(i in 1:nrow(U1)){ - one=matrix(IUX[i,], nrow=1) - beta = crossprod(one) - beta2= beta2+beta - } - beta2<-beta2/delta - - #######get beta3 - beta3=0 - for(i in 1:length(d)){ - one1=matrix(U1TX[i,], nrow=1) - one2=matrix(U1TY[i,], nrow=1) - beta=crossprod(one1,(one2/(d[i]+delta))) - beta3= beta3+beta - } - - ###########get beta4 - beta4=0 - for(i in 1:nrow(U1)){ - one1=matrix(IUX[i,], nrow=1) - one2=matrix(IUY[i,], nrow=1) - beta=crossprod(one1,one2) - beta4= beta4+beta - } - beta4<-beta4/delta - - #######get final beta - zw1 <- ginv(beta1+beta2) - #zw1 <- try(solve(beta1+beta2)) - #if(inherits(zw1, "try-error")){ - #zw1 <- ginv(beta1+beta2) - #} - - zw2=(beta3+beta4) - beta=crossprod(zw1,zw2) - - #----------------------------calculate LL--------------------------------------- - ####part 1 - part11<-n*log(2*3.14) - part12<-0 - for(i in 1:length(d)){ - part12_pre=log(d[i]+delta) - part12= part12+part12_pre - } - part13<- (nrow(U1)-length(d))*log(delta) - part1<- -1/2*(part11+part12+part13) - - ###### part2 - part21<-nrow(U1) - ######part221 - - part221=0 - for(i in 1:length(d)){ - one1=U1TX[i,] - one2=U1TY[i,] - part221_pre=(one2-one1%*%beta)^2/(d[i]+delta) - part221 = part221+part221_pre - } - - part222=0 - for(i in 1:n){ - one1=XU1TX[i,] - one2=yU1TY[i,] - part222_pre=((one2-one1%*%beta)^2)/delta - part222= part222+part222_pre - } - part22<-n*log((1/n)*(part221+part222)) - part2<- -1/2*(part21+part22) - - ################# likihood - LL<-part1+part2 - part1<-0 - part2<-0 - - return(list(beta=beta,delta=delta,LL=LL)) - } - #} # end of Iteration on the range of delta (-5 to 5 in glog scale) - - llresults <- lapply(1:m, beta.optimize.parallel) - - for(i in 1:m){ - if(i == 1){ - beta.save = llresults[[i]]$beta - delta.save = llresults[[i]]$delta - LL.save = llresults[[i]]$LL - }else{ - if(llresults[[i]]$LL > LL.save){ - beta.save = llresults[[i]]$beta - delta.save = llresults[[i]]$delta - LL.save = llresults[[i]]$LL - } - } - } - #--------------------update with the optimum------------------------------------ - beta=beta.save - delta=delta.save - LL=LL.save - - #--------------------calculating Va and Vem------------------------------------- - #sigma_a1 - sigma_a1=0 - for(i in 1:length(d)){ - one1=matrix(U1TX[i,], ncol=1) - one2=matrix(U1TY[i,], nrow=1) - #sigma_a1_pre=(one2-one1%*%beta)^2/(d[i]+delta) - sigma_a1_pre=(one2-crossprod(one1,beta))^2/(d[i]+delta) - sigma_a1= sigma_a1+sigma_a1_pre - } - - ### sigma_a2 - sigma_a2=0 - - for(i in 1:nrow(U1)){ - one1=matrix(IUX[i,], ncol=1) - one2=matrix(IUY[i,], nrow=1) - #sigma_a2_pre<-(one2-one1%*%beta)^2 - sigma_a2_pre<-(one2-crossprod(one1,beta))^2 - sigma_a2= sigma_a2+sigma_a2_pre - } - - sigma_a2<-sigma_a2/delta - sigma_a<- 1/n*(sigma_a1+sigma_a2) - sigma_e<-delta*sigma_a - - return(list(beta=beta, delta=delta, LL=LL, vg=sigma_a, ve=sigma_e)) +#' +`MVP.FaSTLMM.LL` <- function(pheno, snp.pool, X0=NULL, ncpus=2) { + .fastlmm_core(pheno=pheno, snp.pool=snp.pool, X0=X0, ncpus=ncpus) } diff --git a/R/MVP.FarmCPU.r b/R/MVP.FarmCPU.r index af540dc..54df14c 100755 --- a/R/MVP.FarmCPU.r +++ b/R/MVP.FarmCPU.r @@ -70,7 +70,6 @@ memo="MVP.FarmCPU", Prior=NULL, ncpus=2, maxLoop=10, maxLine=5000, threshold.output=.01, converge=1, iteration.output=FALSE, p.threshold=NA, QTN.threshold=0.01, bound=NULL, verbose=TRUE){ - #print("--------------------- Welcome to FarmCPU ----------------------------") 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.") @@ -136,13 +135,13 @@ } #Step 2a: Set prior - myPrior=FarmCPU.Prior(GM=map,P=P,Prior=Prior) + myPrior=.farmcpu_prior(GM=map,P=P,Prior=Prior) #Step 2b: Set bins if(theLoop<=2){ - myBin=FarmCPU.BIN(Y=phe[,c(1,2)],GDP=geno,GDP_index=ind_idx,GM=map,CV=CV,P=myPrior,method=method.bin,b=bin.size,s=bin.selection,theLoop=theLoop,bound=bound,ncpus=ncpus, verbose = verbose) + myBin=.farmcpu_bin(Y=phe[,c(1,2)],GDP=geno,GDP_index=ind_idx,GM=map,CV=CV,P=myPrior,method=method.bin,b=bin.size,s=bin.selection,theLoop=theLoop,bound=bound,ncpus=ncpus, verbose = verbose) }else{ - myBin=FarmCPU.BIN(Y=phe[,c(1,2)],GDP=geno,GDP_index=ind_idx,GM=map,CV=theCV,P=myPrior,method=method.bin,b=bin.size,s=bin.selection,theLoop=theLoop,ncpus=ncpus, verbose = verbose) + myBin=.farmcpu_bin(Y=phe[,c(1,2)],GDP=geno,GDP_index=ind_idx,GM=map,CV=theCV,P=myPrior,method=method.bin,b=bin.size,s=bin.selection,theLoop=theLoop,ncpus=ncpus, verbose = verbose) } #Step 2c: Remove bin dependency @@ -204,7 +203,7 @@ } } - myRemove=FarmCPU.Remove(GDP=geno,GDP_index=ind_idx,GM=map,seqQTN=seqQTN,seqQTN.p=seqQTN.p,threshold=.7) + myRemove=.farmcpu_remove(GDP=geno,GDP_index=ind_idx,GM=map,seqQTN=seqQTN,seqQTN.p=seqQTN.p,threshold=.FARMCPU_LD_THRESHOLD) #Recoding QTNs history seqQTN=myRemove$seqQTN @@ -236,7 +235,6 @@ #Step 3: Screen with bins rm(myBin) - gc() theCV=CV @@ -244,16 +242,16 @@ theCV=cbind(CV,myRemove$bin) } - 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) + 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.") } #Step 4: Background unit substitution if(!isDone){ - myGLM=FarmCPU.SUB(GM=map,GLM=myGLM,QTN=map[myRemove$seqQTN, , drop=FALSE],method=method.sub) + myGLM=.farmcpu_sub(GM=map,GLM=myGLM,QTN=map[myRemove$seqQTN, , drop=FALSE],method=method.sub) }else{ - myGLM=FarmCPU.SUB(GM=map,GLM=myGLM,QTN=map[myRemove$seqQTN, , drop=FALSE],method=method.sub.final) + myGLM=.farmcpu_sub(GM=map,GLM=myGLM,QTN=map[myRemove$seqQTN, , drop=FALSE],method=method.sub.final) } if(!is.null(mrk_idx)){ myGLM$P[-mrk_idx] = NA @@ -266,7 +264,6 @@ colnames(results) = c("effect", "se", "p") if(isDone && !is.null(mrk_idx)) results <- results[mrk_idx, ] } #end of while loop - #print("****************FarmCPU ACCOMPLISHED****************") return(results) }#The MVP.FarmCPU function ends here @@ -292,183 +289,12 @@ #' Output: vg - genetic variance #' Output: ve - residual variance #' -`FarmCPU.FaSTLMM.LL` <- function(pheno, snp.pool, X0=NULL, ncpus=2){ - y=pheno - p=0 - deltaExpStart = -5 - deltaExpEnd = 5 - snp.pool=snp.pool[,] - - if(!is.null(snp.pool)&&any(apply(snp.pool, 2, var)==0)){ - deltaExpStart = 100 - deltaExpEnd = deltaExpStart - } - if(is.null(X0)) { - X0 = matrix(1, nrow(snp.pool), 1) - } - X=X0 - #########SVD of X - K.X.svd <- svd(snp.pool) - d=K.X.svd$d - d=d[d>1e-08] - d=d^2 - U1=K.X.svd$u - U1=U1[,1:length(d)] - #handler of single snp - if(is.null(dim(U1))) U1=matrix(U1,ncol=1) - n=nrow(U1) - U1TX=crossprod(U1,X) - U1TY=crossprod(U1,y) - yU1TY <- y-U1%*%U1TY - XU1TX<- X-U1%*%U1TX - IU = -tcrossprod(U1) - diag(IU) = rep(1,n) + diag(IU) - IUX=crossprod(IU,X) - IUY=crossprod(IU,y) - #Iteration on the range of delta (-5 to 5 in glog scale) - delta.range <- seq(deltaExpStart,deltaExpEnd,by=0.1) - m <- length(delta.range) - #for (m in seq(deltaExpStart,deltaExpEnd,by=0.1)){ - beta.optimize.parallel <- function(ii){ - #p=p+1 - delta <- exp(delta.range[ii]) - #----------------------------calculate beta------------------------------------- - #######get beta1 - beta1=0 - for(i in 1:length(d)){ - one=matrix(U1TX[i,], nrow=1) - beta=crossprod(one,(one/(d[i]+delta))) #This is not real beta, confusing - beta1= beta1+beta - } - - #######get beta2 - beta2=0 - for(i in 1:nrow(U1)){ - one=matrix(IUX[i,], nrow=1) - beta = crossprod(one) - beta2= beta2+beta - } - beta2<-beta2/delta - - #######get beta3 - beta3=0 - for(i in 1:length(d)){ - one1=matrix(U1TX[i,], nrow=1) - one2=matrix(U1TY[i,], nrow=1) - beta=crossprod(one1,(one2/(d[i]+delta))) - beta3= beta3+beta - } - - ###########get beta4 - beta4=0 - for(i in 1:nrow(U1)){ - one1=matrix(IUX[i,], nrow=1) - one2=matrix(IUY[i,], nrow=1) - beta=crossprod(one1,one2) - beta4= beta4+beta - } - beta4<-beta4/delta - - #######get final beta - zw1 <- ginv(beta1+beta2) - #zw1 <- try(solve(beta1+beta2)) - #if(inherits(zw1, "try-error")){ - #zw1 <- ginv(beta1+beta2) - #} - - zw2=(beta3+beta4) - beta=crossprod(zw1,zw2) - - #----------------------------calculate LL--------------------------------------- - ####part 1 - part11<-n*log(2*3.14) - part12<-0 - for(i in 1:length(d)){ - part12_pre=log(d[i]+delta) - part12= part12+part12_pre - } - part13<- (nrow(U1)-length(d))*log(delta) - part1<- -1/2*(part11+part12+part13) - - ###### part2 - part21<-nrow(U1) - ######part221 - - part221=0 - for(i in 1:length(d)){ - one1=U1TX[i,] - one2=U1TY[i,] - part221_pre=(one2-one1%*%beta)^2/(d[i]+delta) - part221 = part221+part221_pre - } - - part222=0 - for(i in 1:n){ - one1=XU1TX[i,] - one2=yU1TY[i,] - part222_pre=((one2-one1%*%beta)^2)/delta - part222= part222+part222_pre - } - part22<-n*log((1/n)*(part221+part222)) - part2<- -1/2*(part21+part22) - - ################# likihood - LL<-part1+part2 - part1<-0 - part2<-0 - - return(list(beta=beta,delta=delta,LL=LL)) - } - llresults <- lapply(1:m, beta.optimize.parallel) - for(i in 1:m){ - if(i == 1){ - beta.save = llresults[[i]]$beta - delta.save = llresults[[i]]$delta - LL.save = llresults[[i]]$LL - }else{ - if(llresults[[i]]$LL > LL.save){ - beta.save = llresults[[i]]$beta - delta.save = llresults[[i]]$delta - LL.save = llresults[[i]]$LL - } - } - } - #--------------------update with the optimum------------------------------------ - beta=beta.save - delta=delta.save - LL=LL.save - - #--------------------calculating Va and Vem------------------------------------- - #sigma_a1 - sigma_a1=0 - for(i in 1:length(d)){ - one1=matrix(U1TX[i,], ncol=1) - one2=matrix(U1TY[i,], nrow=1) - #sigma_a1_pre=(one2-one1%*%beta)^2/(d[i]+delta) - sigma_a1_pre=(one2-crossprod(one1,beta))^2/(d[i]+delta) - sigma_a1= sigma_a1+sigma_a1_pre - } - - ### sigma_a2 - sigma_a2=0 - - for(i in 1:nrow(U1)){ - one1=matrix(IUX[i,], ncol=1) - one2=matrix(IUY[i,], nrow=1) - #sigma_a2_pre<-(one2-one1%*%beta)^2 - sigma_a2_pre<-(one2-crossprod(one1,beta))^2 - sigma_a2= sigma_a2+sigma_a2_pre - } - - sigma_a2<-sigma_a2/delta - sigma_a<- 1/n*(sigma_a1+sigma_a2) - sigma_e<-delta*sigma_a - - return(list(beta=beta, delta=delta, LL=LL, vg=sigma_a, ve=sigma_e)) +`FarmCPU.FaSTLMM.LL` <- function(pheno, snp.pool, X0=NULL, ncpus=2) { + .fastlmm_core(pheno=pheno, snp.pool=snp.pool, X0=X0, ncpus=ncpus) } -#' FarmCPU.BIN +#' .farmcpu_bin #' #' Last update: March 28, 2017 #' Requirement: Y, GDP, and CV must have same taxa order. GDP and GM must have the same order on SNP @@ -483,7 +309,7 @@ #' @param GM SNP map information, m by 3 matrix, m is marker size, the three columns are SNP_ID, Chr, and Pos #' @param CV covariates, n by c matrix, n is sample size, c is number of covariates #' @param P start p values for all SNPs -#' @param method two options, 'static' or 'FaST-LMM' +#' @param method method for selecting the optimal bin: 'static', 'FaST-LMM', 'EMMA', or 'GEMMA' #' @param b bin sizes for all iterations, a vector, the bin size is always from large to small #' @param s number of selected bins in each iteration, a vector #' @param theLoop iteration number @@ -495,10 +321,9 @@ #' Output: seqQTN - an s by 1 vecter for index of QTNs on GM file #' #' @keywords internal -FarmCPU.BIN <- +.farmcpu_bin <- function(Y=NULL, GDP=NULL, GDP_index=NULL, GM=NULL, CV=NULL, P=NULL, method="EMMA", b=c(5e5,5e6,5e7), s=seq(10,100,10), theLoop=NULL, bound=NULL, ncpus=2, verbose=TRUE){ - #print("FarmCPU.BIN Started") - + if(is.null(P)) return(list(bin=NULL,binmap=NULL,seqQTN=NULL)) if(nrow(GM) == ncol(GDP)){ @@ -524,15 +349,13 @@ FarmCPU.BIN <- if(optimumable){ s[s>bound]=bound - #print("optimizing possible QTNs...") #GP=cbind(GM,P,NA,NA,NA) - #mySpecify=FarmCPU.Specify(GI=GM, GP=GP, bin.size=b, inclosure.size=s) + #mySpecify=.farmcpu_specify(GI=GM, GP=GP, bin.size=b, inclosure.size=s) #seqQTN=which(mySpecify$index==TRUE) } #Method of static if(method=="static"&optimumable){ - #print("Via static") if(theLoop==2){ b=b[3] }else if(theLoop==3){ @@ -544,28 +367,27 @@ FarmCPU.BIN <- s[s>bound]=bound logging.log("Optimizing Pseudo QTNs...", "\n", verbose = verbose) GP=cbind(GM,P,NA,NA,NA) - mySpecify=FarmCPU.Specify(GI=GM,GP=GP,bin.size=b,inclosure.size=s) + mySpecify=.farmcpu_specify(GI=GM,GP=GP,bin.size=b,inclosure.size=s) seqQTN.save=which(mySpecify$index==TRUE) } #Method of optimum: FaST-LMM #============================Optimize by FaST-LMM============================================ if(method=="FaST-LMM"&optimumable){ - #print("c(bin.size, bin.selection, -2LL, VG, VE)") logging.log("Optimizing Pseudo QTNs...", "\n", verbose = verbose) count=0 for (bin in b){ for (inc in s){ count=count+1 GP=cbind(GM,P,NA,NA,NA) - mySpecify=FarmCPU.Specify(GI=GM,GP=GP,bin.size=bin,inclosure.size=inc) + mySpecify=.farmcpu_specify(GI=GM,GP=GP,bin.size=bin,inclosure.size=inc) seqQTN=which(mySpecify$index==TRUE) if(nrow(GM) == ncol(GDP)){ GK=GDP[GDP_index, seqQTN] }else{ GK=t(GDP[seqQTN, GDP_index]) } - myBurger=FarmCPU.Burger(Y=Y[,1:2], CV=CV, GK=GK, ncpus=ncpus, method=method) + myBurger=.farmcpu_burger(Y=Y[,1:2], CV=CV, GK=GK, ncpus=ncpus, method=method) myREML=myBurger$REMLs myVG=myBurger$vg #it is unused myVE=myBurger$ve #it is unused @@ -596,7 +418,6 @@ FarmCPU.BIN <- #Method of optimum: EMMA #============================Optimize by EMMA============================================ if(method=="EMMA"&optimumable){ - #print("c(bin.size, bin.selection, -2LL, VG, VE)") logging.log("Optimizing Pseudo QTNs...", "\n", verbose = verbose) m <- length(b)*length(s) inc.index = rep(c(1:length(s)), length(b)) @@ -606,14 +427,14 @@ FarmCPU.BIN <- bin = b[bin.index] inc = s[inc.index[ii]] GP=cbind(GM,P,NA,NA,NA) - mySpecify=FarmCPU.Specify(GI=GM,GP=GP,bin.size=bin,inclosure.size=inc) + mySpecify=.farmcpu_specify(GI=GM,GP=GP,bin.size=bin,inclosure.size=inc) seqQTN=which(mySpecify$index==TRUE) if(nrow(GM) == ncol(GDP)){ GK=GDP[GDP_index, seqQTN] }else{ GK=t(GDP[seqQTN, GDP_index]) } - myBurger=FarmCPU.Burger(Y=Y[,1:2], CV=CV, GK=GK, ncpus=ncpus, method=method) + myBurger=.farmcpu_burger(Y=Y[,1:2], CV=CV, GK=GK, ncpus=ncpus, method=method) myREML=myBurger$REMLs myVG=myBurger$vg #it is unused myVE=myBurger$ve #it is unused @@ -622,24 +443,17 @@ FarmCPU.BIN <- } llresults <- lapply(1:m, seqQTN.optimize.parallel) - for(i in 1:m){ - if(i == 1){ - seqQTN.save = llresults[[i]]$seqQTN - myREML.save = llresults[[i]]$myREML - }else{ - if(llresults[[i]]$myREML < myREML.save){ - seqQTN.save = llresults[[i]]$seqQTN - myREML.save = llresults[[i]]$myREML - } - } - } + # Find result with minimum REML + REML_values <- sapply(llresults, function(x) x$myREML) + best_idx <- which.min(REML_values) + seqQTN.save <- llresults[[best_idx]]$seqQTN + myREML.save <- llresults[[best_idx]]$myREML } #Method of optimum: GEMMA #can not be used to provide REML #============================Optimize by GEMMA============================================ if(method=="GEMMA"&optimumable){ - #print("c(bin.size, bin.selection, -2LL, VG, VE)") logging.log("Optimizing Pseudo QTNs...\n", verbose = verbose) m <- length(b)*length(s) @@ -647,14 +461,14 @@ FarmCPU.BIN <- bin = floor((ii-0.1)/length(s)) + 1 inc = rep(c(1:length(s)), length(b)) GP=cbind(GM,P,NA,NA,NA) - mySpecify=FarmCPU.Specify(GI=GM,GP=GP,bin.size=bin[ii],inclosure.size=inc[ii]) + mySpecify=.farmcpu_specify(GI=GM,GP=GP,bin.size=bin[ii],inclosure.size=inc[ii]) seqQTN=which(mySpecify$index==TRUE) if(nrow(GM) == ncol(GDP)){ GK=GDP[GDP_index, seqQTN] }else{ GK=t(GDP[seqQTN, GDP_index]) } - myBurger=FarmCPU.Burger(Y=Y[,1:2], CV=CV, GK=GK, ncpus=ncpus, method=method) + myBurger=.farmcpu_burger(Y=Y[,1:2], CV=CV, GK=GK, ncpus=ncpus, method=method) myREML=myBurger$REMLs myVG=myBurger$vg #it is unused myVE=myBurger$ve #it is unused @@ -663,21 +477,15 @@ FarmCPU.BIN <- } llresults <- lapply(1:m, seqQTN.optimize.parallel) - for(i in 1:m){ - if(i == 1){ - seqQTN.save = llresults[[i]]$seqQTN - myREML.save = llresults[[i]]$myREML - }else{ - if(llresults[[i]]$myREML < myREML.save){ - seqQTN.save = llresults[[i]]$seqQTN - myREML.save = llresults[[i]]$myREML - } - } - } + # Find result with minimum REML + REML_values <- sapply(llresults, function(x) x$myREML) + best_idx <- which.min(REML_values) + seqQTN.save <- llresults[[best_idx]]$seqQTN + myREML.save <- llresults[[best_idx]]$myREML } return(list(seqQTN=seqQTN.save)) - }#The function FarmCPU.BIN ends here + }#The function .farmcpu_bin ends here #' To get indicator (TURE or FALSE) for GI based on GP @@ -698,9 +506,8 @@ FarmCPU.BIN <- #' @param MaxBP maximum base pairs for each chromosome #' #' @return theIndex: a vector indicating if the SNPs in GI belong to QTN or not -FarmCPU.Specify <- +.farmcpu_specify <- function(GI=NULL, GP=NULL, bin.size=10000000, inclosure.size=NULL, MaxBP=1e10){ - #print("Specification in process...") if(is.null(GP))return (list(index=NULL,BP=NULL)) #set inclosure bin in GP @@ -757,7 +564,7 @@ FarmCPU.Specify <- myList=list(index=theIndex,CB=ID.GP) return (list(index=theIndex,CB=ID.GP)) - } #end of FarmCPU.Specify + } #end of .farmcpu_specify #' To quickly sovel LM with one variable substitute multiple times @@ -789,9 +596,8 @@ FarmCPU.Specify <- #' Output: B - effect of each SNP #' #' @keywords internal -FarmCPU.LM <- +.farmcpu_lm <- 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) #Constant section (non individual marker specific) @@ -849,7 +655,6 @@ FarmCPU.LM <- sepred = se[-1] } - #print(ncpus) logging.log("scanning...", "\n", verbose = verbose) # P <- NULL @@ -870,10 +675,10 @@ FarmCPU.LM <- }) 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)) - } #end of FarmCPU.LM function + } #end of .farmcpu_lm function -#' FarmCPU.Burger +#' .farmcpu_burger #' #' Last update: Dec 21, 2016 #' To calculate likelihood, variances and ratio, revised by Xiaolei based on GAPIT.Burger function from GAPIT package @@ -893,7 +698,7 @@ FarmCPU.LM <- #' Output: delta - exp(root) #' #' @keywords internal -FarmCPU.Burger <- +.farmcpu_burger <- function(Y=NULL,CV=NULL,GK=NULL,ncpus=2, method="FaST-LMM"){ if(!is.null(CV)){ CV=as.matrix(CV)#change CV to a matrix when it is a vector @@ -935,12 +740,11 @@ FarmCPU.Burger <- ve=myEMMAREML$ve } - #print("FarmCPU.Burger succeed!") return (list(REMLs=REMLs, vg=vg, ve=ve, delta=delta)) - } #end of FarmCPU.Burger + } #end of .farmcpu_burger -#' FarmCPU.SUB +#' .farmcpu_sub #' #' Last update: Febuary 26, 2013 #' Requirement: P has row name of SNP. s<=t. covariates of QTNs are next to SNP @@ -955,12 +759,10 @@ FarmCPU.Burger <- #' @return #' Output: GLM$P - Updated p-values by substitution process #' Output: GLM$B - Updated effects by substitution process -FarmCPU.SUB <- +.farmcpu_sub <- function(GM=NULL,GLM=NULL,QTN=NULL,method="mean"){ if(is.null(GLM$P)) return(NULL) #P is required if(is.null(QTN)) return(NULL) #QTN is required - #print("FarmCPU.SUB Started") - #print(length(QTN)) #if(length(QTN)==3){ #QTN=QTN[1] #}else{ @@ -995,7 +797,7 @@ FarmCPU.SUB <- GLM$S[position, ] = GLM$sepred } return(GLM) - }#The function FarmCPU.SUB ends here + }#The function .farmcpu_sub ends here #' Remove bins that are highly correlated @@ -1018,7 +820,7 @@ FarmCPU.SUB <- #' Output: seqQTN - s0 by 1 vecter for index of QTN on GM #' Relationship: bin=GDP[,c(seqQTN)], binmap=GM[seqQTN,], s0<=s #' @keywords internal -FarmCPU.Remove <- +.farmcpu_remove <- function(GDP=NULL, GDP_index=NULL, GM=NULL, seqQTN=NULL, seqQTN.p=NULL, threshold=.99){ if(is.null(seqQTN))return(list(bin=NULL,binmap=NULL,seqQTN=NULL)) @@ -1039,16 +841,9 @@ FarmCPU.Remove <- cb=as.numeric(binmap[, 2, drop=FALSE])*hugeNum+as.numeric(binmap[, 3, drop=FALSE])#create ID for chromosome and bp cb.unique=unique(cb) - #print("debuge") - #print(cb) - #print(cb.unique) - index=match(cb.unique,cb,nomatch = 0) seqQTN=seqQTN[index] - #print("Number of bins after chr and bp fillter") - n=length(seqQTN) #update n - #Set sample ratio=.1 maxNum=100000 @@ -1118,7 +913,7 @@ FarmCPU.Remove <- binmap=GM[seqQTN, , drop=FALSE] return(list(bin=bin, binmap=binmap, seqQTN=seqQTN)) - }#The function FarmCPU.Remove ends here + }#The function .farmcpu_remove ends here #' Set prior on existing p value @@ -1136,20 +931,16 @@ FarmCPU.Remove <- #' Output: P - updated P value by prior information #' #' @keywords internal -FarmCPU.Prior <- +.farmcpu_prior <- function(GM, P=NULL, Prior=NULL){ - #print("FarmCPU.Prior Started") - + if(is.null(Prior)& is.null(P))return(P) #get prior position if(!is.null(Prior)) index=match(Prior[, 1, drop=FALSE],GM[, 1, drop=FALSE],nomatch = 0) - #if(is.null(P)) P=runif(nrow(GM)) #set random p value if not provided (This is not helpful) - #print("debug set prior a") - #Get product with prior if provided if(!is.null(Prior) & !is.null(P) )P[index]=P[index]*Prior[, 4, drop=FALSE] return(P) - }#The function FarmCPU.Prior ends here + }#The function .farmcpu_prior ends here diff --git a/R/MVP.GLM.r b/R/MVP.GLM.r index 04265ca..41322e4 100644 --- a/R/MVP.GLM.r +++ b/R/MVP.GLM.r @@ -85,12 +85,8 @@ function( CV <- CV[, CV.index, drop=FALSE] X0 <- cbind(matrix(1, n), CV) } - X0 <- as.matrix(X0) - tX0X0 <- crossprod(X0); - iX0X0 <- try(solve(tX0X0), silent = TRUE) - if(inherits(iX0X0, "try-error")) { - iX0X0 <- MASS::ginv(tX0X0) - } + tX0X0 <- crossprod(X0) + iX0X0 <- .safe_solve(tX0X0) logging.log("scanning...\n", verbose = verbose) diff --git a/R/MVP.HE.Vg.Ve.R b/R/MVP.HE.Vg.Ve.R index b08b870..50cdc12 100755 --- a/R/MVP.HE.Vg.Ve.R +++ b/R/MVP.HE.Vg.Ve.R @@ -94,7 +94,6 @@ MVP.HE.Vg.Ve <- function(y, X, K) { StandardizeVector <- function(x) { m = mean(x) v = sum((x - m)^2) / length(x) - v = v - m * m x = (x - m) / sqrt(v) return (x) } @@ -161,41 +160,12 @@ MVP.HE.Vg.Ve <- function(y, X, K) { log_sigma2 = NULL for (i in 1:length(v_sigma2)) { if (v_sigma2[i] <= 0){ - log_sigma2[i] = log(0.1) + log_sigma2[i] = .HE_DEFAULT_LOG_SIGMA2 } else { log_sigma2[i] = log(v_sigma2[i]) } } - LogRL_dev1 <- function(log_sigma2, parms) { - - y = parms$y - X = parms$X - K = parms$K - n = nrow(K) - - P = K * exp(log_sigma2[1]) + diag(n) * exp(log_sigma2[2]) - - # calculate H^{-1} - P = solve(P) - - # calculate P=H^{-1}-H^{-1}X(X^TH^{-1}X)^{-1}X^TH^{-1} - HiW = P %*% X - WtHiW = crossprod(X, HiW) - WtHiWi = solve(WtHiW) - - WtHiWiWtHi = tcrossprod(WtHiWi, HiW) - P = P - HiW %*% WtHiWiWtHi - - # calculate Py, KPy, PKPy - Py = P %*% matrix(y) - KPy = K %*% Py - - # calculate dev1=-0.5*trace(PK_i)+0.5*yPKPy - c(para1 = (-0.5 * sum(t(P) * K) + 0.5 * crossprod(Py, KPy)) * exp(log_sigma2[1]), - para2 = (-0.5 * sum(diag(P)) + 0.5 * crossprod(Py)) * exp(log_sigma2[2])) - } - vg = exp(log_sigma2[1]) ve = exp(log_sigma2[2]) delta = ve / vg diff --git a/R/MVP.K.VanRaden.r b/R/MVP.K.VanRaden.r index fe80911..b429f69 100755 --- a/R/MVP.K.VanRaden.r +++ b/R/MVP.K.VanRaden.r @@ -188,6 +188,6 @@ function( # } # ) #print("K Preparation is Done!") - logging.log("Deriving relationship matrix successfully", "\n", verbose = verbose); gc() + logging.log("Deriving relationship matrix successfully", "\n", verbose = verbose) return(K) }#end of MVP.k.VanRaden function diff --git a/R/MVP.Utility.r b/R/MVP.Utility.r index d959e65..0cb6767 100755 --- a/R/MVP.Utility.r +++ b/R/MVP.Utility.r @@ -118,30 +118,6 @@ print_bar <- function(i, } } }, - # "type2"={ - # if(inherits(parallel:::mcfork(), "masterProcess")) { - # progress <- 0.0 - # while(progress < n && !isIncomplete(tmp.file)){ - # msg <- readBin(tmp.file, "double") - # progress <- progress + as.numeric(msg) - # print.len <- round(symbol.len * progress / n) - # if(fixed.points){ - # if(progress %in% round(points * n / 100)){ - # logging.log(paste("\r", - # paste(c(symbol.head, rep("-", print.len), symbol.tail), collapse=""), - # paste(rep(" ", symbol.len-print.len), collapse=""), - # sprintf("%.2f%%", progress * 100 / n), sep="")) - # } - # }else{ - # logging.log(paste("\r", - # paste(c(symbol.head, rep("-", print.len), symbol.tail), collapse=""), - # paste(rep(" ", symbol.len-print.len), collapse=""), - # sprintf("%.2f%%", progress * 100 / n), sep="")) - # } - # } - # parallel:::mcexit() - # } - # }, "type3"={ progress <- readBin(tmp.file, "double") + 1 writeBin(progress, tmp.file) @@ -191,14 +167,7 @@ print_accomplished <- function(width = 60, verbose = TRUE) { #' @param version short label, bottom-right of logo #' @param authors #' @param contact email or website -#' @param line 1, 2, or char #' @param width banner width -#' -#' @examples -#' welcome <- "Welcome to MVP" -#' title <- "A Memory-efficient, Visualization-enhanced, and Parallel-accelerated Tool For GWAS" -#' authors <- "Authors: Lilin Yin, Haohao Zhang, and Xiaolei Liu" -#' print_info(welcome = welcome, title = title, logo = NULL, authors = authors, contact = NULL, linechar = '=', width = width) print_info <- function(welcome=NULL, title=NULL, short_title=NULL, logo=NULL, version=NULL, authors=NULL, contact=NULL, linechar = '=', width=NULL, verbose=TRUE) { msg <- c() # width @@ -294,7 +263,7 @@ make_line <- function(string, width, linechar = " ", align = "center", margin = #' #' Build date: Oct 22, 2018 #' Last update: Dec 12, 2018 -#' by using base::strwarp. +#' by using base::strwrap. #' #' @keywords internal #' @author Haohao Zhang @@ -363,6 +332,12 @@ load_if_installed <- function(package) { } +.safe_solve <- function(X) { + result <- tryCatch(solve(X), error = function(e) MASS::ginv(X)) + result +} + + mkl_env <- function(exprs, threads = 1) { if (load_if_installed("RevoUtilsMath")) { math.cores <- eval(parse(text = "getMKLthreads()")) diff --git a/R/MVP.r b/R/MVP.r index 0c81c75..ba010e9 100755 --- a/R/MVP.r +++ b/R/MVP.r @@ -338,7 +338,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, mrk_bycol = MrkByCol, 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) 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) @@ -353,7 +353,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, mrk_bycol = MrkByCol, 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) 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/man/FarmCPU.FaSTLMM.LL.Rd b/man/FarmCPU.FaSTLMM.LL.Rd new file mode 100644 index 0000000..58577b3 --- /dev/null +++ b/man/FarmCPU.FaSTLMM.LL.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.FarmCPU.r +\name{FarmCPU.FaSTLMM.LL} +\alias{FarmCPU.FaSTLMM.LL} +\title{FarmCPU.FaSTLMM.LL +Evaluation of the maximum likelihood using FaST-LMM method + +Last update: January 11, 2017 +Requirement: pheno, snp.pool, and X0 must have same taxa order. +Requirement: No missing data} +\usage{ +FarmCPU.FaSTLMM.LL(pheno, snp.pool, X0 = NULL, ncpus = 2) +} +\arguments{ +\item{pheno}{a two-column phenotype matrix} + +\item{snp.pool}{matrix for pseudo QTNs} + +\item{X0}{covariates matrix} + +\item{ncpus}{number of threads used for parallel computation} +} +\value{ +Output: beta - beta effect +Output: delta - delta value +Output: LL - log-likelihood +Output: vg - genetic variance +Output: ve - residual variance +} +\description{ +FarmCPU.FaSTLMM.LL +Evaluation of the maximum likelihood using FaST-LMM method + +Last update: January 11, 2017 +Requirement: pheno, snp.pool, and X0 must have same taxa order. +Requirement: No missing data +} +\author{ +Xiaolei Liu (modified) +} diff --git a/man/MVP.Rd b/man/MVP.Rd index e0c737f..4b15348 100755 --- a/man/MVP.Rd +++ b/man/MVP.Rd @@ -106,10 +106,10 @@ MVP( \item{verbose}{whether to print detail.} } \value{ -Output: MVP.return$map - SNP map information, SNP name, Chr, Pos -Output: MVP.return$glm.results - p-values obtained by GLM method -Output: MVP.return$mlm.results - p-values obtained by MLM method -Output: MVP.return$farmcpu.results - p-values obtained by FarmCPU method +Output: MVP.return$map - SNP map information, SNP name, Chr, Pos, A1, A2, MAF +Output: MVP.return$glm.results - effect (at A2) and p-values obtained by GLM method +Output: MVP.return$mlm.results - effect (at A2) and p-values obtained by MLM method +Output: MVP.return$farmcpu.results - effect (at A2) and p-values obtained by FarmCPU method } \description{ Object 1: To perform GWAS using General Linear Model (GLM), Mixed Linear Model (MLM), and FarmCPU model diff --git a/man/dot-farmcpu_bin.Rd b/man/dot-farmcpu_bin.Rd new file mode 100644 index 0000000..1209d02 --- /dev/null +++ b/man/dot-farmcpu_bin.Rd @@ -0,0 +1,62 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.FarmCPU.r +\name{.farmcpu_bin} +\alias{.farmcpu_bin} +\title{.farmcpu_bin} +\usage{ +.farmcpu_bin( + Y = NULL, + GDP = NULL, + GDP_index = NULL, + GM = NULL, + CV = NULL, + P = NULL, + method = "EMMA", + b = c(5e+05, 5e+06, 5e+07), + s = seq(10, 100, 10), + theLoop = NULL, + bound = NULL, + ncpus = 2, + verbose = TRUE +) +} +\arguments{ +\item{Y}{a n by 2 matrix, the fist column is taxa id and the second is trait} + +\item{GDP}{genotype} + +\item{GDP_index}{the index of effective genotyped individuals} + +\item{GM}{SNP map information, m by 3 matrix, m is marker size, the three columns are SNP_ID, Chr, and Pos} + +\item{CV}{covariates, n by c matrix, n is sample size, c is number of covariates} + +\item{P}{start p values for all SNPs} + +\item{method}{method for selecting the optimal bin: 'static', 'FaST-LMM', 'EMMA', or 'GEMMA'} + +\item{b}{bin sizes for all iterations, a vector, the bin size is always from large to small} + +\item{s}{number of selected bins in each iteration, a vector} + +\item{theLoop}{iteration number} + +\item{bound}{maximum number of SNPs selected as pseudo QTNs in each iteration} + +\item{ncpus}{number of threads used for parallele computation} + +\item{verbose}{whether to print detail.} +} +\value{ +Output: seqQTN - an s by 1 vecter for index of QTNs on GM file +} +\description{ +Last update: March 28, 2017 +Requirement: Y, GDP, and CV must have same taxa order. GDP and GM must have the same order on SNP +Requirement: P and GM are in the same order +Requirement: No missing data +} +\author{ +Xiaolei Liu and Zhiwu Zhang +} +\keyword{internal} diff --git a/man/dot-farmcpu_burger.Rd b/man/dot-farmcpu_burger.Rd new file mode 100644 index 0000000..43c321a --- /dev/null +++ b/man/dot-farmcpu_burger.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.FarmCPU.r +\name{.farmcpu_burger} +\alias{.farmcpu_burger} +\title{.farmcpu_burger} +\usage{ +.farmcpu_burger(Y = NULL, CV = NULL, GK = NULL, ncpus = 2, method = "FaST-LMM") +} +\arguments{ +\item{Y}{phenotype, n by t matrix, n is sample size, t is number of phenotypes} + +\item{CV}{covariates, n by c matrix, n is sample size, c is number of covariates} + +\item{GK}{Genotype data in numerical format, taxa goes to row and snp go to columns} + +\item{ncpus}{number of threads used for parallele computation} + +\item{method}{two options for estimating variance components, 'FaST-LMM' or 'EMMA'} +} +\value{ +Output: REMLs - maximum log likelihood +Output: vg - genetic variance +Output: ve - residual variance +Output: delta - exp(root) +} +\description{ +Last update: Dec 21, 2016 +To calculate likelihood, variances and ratio, revised by Xiaolei based on GAPIT.Burger function from GAPIT package +} +\author{ +Xiaolei Liu and Zhiwu Zhang +} +\keyword{internal} diff --git a/man/dot-farmcpu_lm.Rd b/man/dot-farmcpu_lm.Rd new file mode 100644 index 0000000..f063b1c --- /dev/null +++ b/man/dot-farmcpu_lm.Rd @@ -0,0 +1,56 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.FarmCPU.r +\name{.farmcpu_lm} +\alias{.farmcpu_lm} +\title{To quickly sovel LM with one variable substitute multiple times} +\usage{ +.farmcpu_lm( + y, + w = NULL, + GDP, + GDP_index = NULL, + GDP_mrk_bycol = TRUE, + maxLine = 5000, + ncpus = 2, + npc = 0, + verbose = TRUE +) +} +\arguments{ +\item{y}{one column matrix, dependent variable} + +\item{w}{covariates, n by c matrix, n is sample size, c is number of covariates} + +\item{GDP}{genotype} + +\item{GDP_index}{index of effective genotyped individuals} + +\item{GDP_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{ncpus}{number of threads used for parallele computation} + +\item{npc}{number of covariates without pseudo QTNs} + +\item{verbose}{whether to print detail.} +} +\value{ +Output: P - p-value of each SNP +Output: betapred - effects of pseudo QTNs +Output: B - effect of each SNP +} +\description{ +Start date: March 1, 2013 +Last update: March 6, 2013 +Strategy: +1. Separate constant covariates (w) and dynamic coveriates (x) +2. Build non-x related only once +3. Use apply to iterate x +4. Derive dominance indicate d from additive indicate (x) mathmaticaly +5. When d is not estimable, continue to test x +} +\author{ +Xiaolei Liu and Zhiwu Zhang +} +\keyword{internal} diff --git a/man/dot-farmcpu_prior.Rd b/man/dot-farmcpu_prior.Rd new file mode 100644 index 0000000..dce8b70 --- /dev/null +++ b/man/dot-farmcpu_prior.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.FarmCPU.r +\name{.farmcpu_prior} +\alias{.farmcpu_prior} +\title{Set prior on existing p value} +\usage{ +.farmcpu_prior(GM, P = NULL, Prior = NULL) +} +\arguments{ +\item{GM}{an m by 3 matrix for SNP name, chromosome and BP} + +\item{P}{an m by 1 matrix containing probability} + +\item{Prior}{an s by 4 matrix for SNP name, chromosome, BP and p-value} +} +\value{ +Output: P - updated P value by prior information +} +\description{ +Last update: March 10, 2013 +Requirement: P and GM are in the same order, Prior is part of GM except P value +} +\author{ +Xiaolei Liu and Zhiwu Zhang +} +\keyword{internal} diff --git a/man/dot-farmcpu_remove.Rd b/man/dot-farmcpu_remove.Rd new file mode 100644 index 0000000..c67da12 --- /dev/null +++ b/man/dot-farmcpu_remove.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.FarmCPU.r +\name{.farmcpu_remove} +\alias{.farmcpu_remove} +\title{Remove bins that are highly correlated} +\usage{ +.farmcpu_remove( + GDP = NULL, + GDP_index = NULL, + GM = NULL, + seqQTN = NULL, + seqQTN.p = NULL, + threshold = 0.99 +) +} +\arguments{ +\item{GDP}{genotype} + +\item{GDP_index}{the index of effective genotyped individuals} + +\item{GM}{SNP map information, m by 3 matrix, m is marker size, the three columns are SNP_ID, Chr, and Pos} + +\item{seqQTN}{s by 1 vecter for index of QTN on GM} + +\item{seqQTN.p}{p value of each @param seqQTN} + +\item{threshold}{pearson correlation threshold for remove correlated markers} +} +\value{ +Output: bin - n by s0 matrix of genotype +Output: binmap - s0 by 3 matrix for map of bin +Output: seqQTN - s0 by 1 vecter for index of QTN on GM +Relationship: bin=GDP[,c(seqQTN)], binmap=GM[seqQTN,], s0<=s +} +\description{ +Last update: March 4, 2013 +Requirement: GDP and GM have the same order on SNP +} +\author{ +Xiaolei Liu and Zhiwu Zhang +} +\keyword{internal} diff --git a/man/dot-farmcpu_specify.Rd b/man/dot-farmcpu_specify.Rd new file mode 100644 index 0000000..5ef6a33 --- /dev/null +++ b/man/dot-farmcpu_specify.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.FarmCPU.r +\name{.farmcpu_specify} +\alias{.farmcpu_specify} +\title{To get indicator (TURE or FALSE) for GI based on GP} +\usage{ +.farmcpu_specify( + GI = NULL, + GP = NULL, + bin.size = 1e+07, + inclosure.size = NULL, + MaxBP = 1e+10 +) +} +\arguments{ +\item{GI}{Data frame with three columns (SNP name, chr and base position)} + +\item{GP}{Data frame with seven columns (SNP name, chr and base position, P, MAF, N, effect)} + +\item{bin.size}{a value of @param b in 'FarmCPU.bin' function} + +\item{inclosure.size}{a value of @param s in 'FarmCPU.bin' function} + +\item{MaxBP}{maximum base pairs for each chromosome} +} +\value{ +theIndex: a vector indicating if the SNPs in GI belong to QTN or not +} +\description{ +Last update: January 26, 2017 +Strategy +1.set bins for all snps in GP +2.keep the snp with smallest P value in each bin, record SNP ID +3.Search GI for SNP with SNP ID from above +4.return the position for SNP selected +} +\author{ +Zhiwu Zhang +} diff --git a/man/dot-farmcpu_sub.Rd b/man/dot-farmcpu_sub.Rd new file mode 100644 index 0000000..f95e646 --- /dev/null +++ b/man/dot-farmcpu_sub.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.FarmCPU.r +\name{.farmcpu_sub} +\alias{.farmcpu_sub} +\title{.farmcpu_sub} +\usage{ +.farmcpu_sub(GM = NULL, GLM = NULL, QTN = NULL, method = "mean") +} +\arguments{ +\item{GM}{SNP map information, m by 3 matrix, m is marker size, the three columns are SNP_ID, Chr, and Pos} + +\item{GLM}{FarmCPU.GLM object} + +\item{QTN}{s by 3 matrix for SNP name, chromosome and BP} + +\item{method}{options are "penalty", "reward","mean","median",and "onsite"} +} +\value{ +Output: GLM$P - Updated p-values by substitution process +Output: GLM$B - Updated effects by substitution process +} +\description{ +Last update: Febuary 26, 2013 +Requirement: P has row name of SNP. s<=t. covariates of QTNs are next to SNP +} +\author{ +Xiaolei Liu and Zhiwu Zhang +} diff --git a/man/format_time.Rd b/man/format_time.Rd new file mode 100644 index 0000000..d6cc2e9 --- /dev/null +++ b/man/format_time.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.Utility.r +\name{format_time} +\alias{format_time} +\title{format time} +\usage{ +format_time(x) +} +\arguments{ +\item{x}{seconds} +} +\description{ +format time +} +\author{ +Lilin Yin +} +\keyword{internal} diff --git a/man/make_line.Rd b/man/make_line.Rd new file mode 100644 index 0000000..cb0f69a --- /dev/null +++ b/man/make_line.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.Utility.r +\name{make_line} +\alias{make_line} +\title{make line} +\usage{ +make_line(string, width, linechar = " ", align = "center", margin = 1) +} +\description{ +Build date: Dec 12, 2018 +Last update: Dec 12, 2018 +} +\author{ +Haohao Zhang +} +\keyword{internal} diff --git a/man/paste_label.Rd b/man/paste_label.Rd new file mode 100644 index 0000000..63a9f48 --- /dev/null +++ b/man/paste_label.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.Utility.r +\name{paste_label} +\alias{paste_label} +\title{Paste label to a line} +\usage{ +paste_label(line, label, side = "right", margin = 2) +} +\arguments{ +\item{line}{long text} + +\item{label}{short label} + +\item{side}{"right" or "left"} + +\item{margin}{default 2} +} +\description{ +Build date: Oct 22, 2018 +Last update: Oct 22, 2018 +} +\author{ +Haohao Zhang +} +\keyword{internal} diff --git a/man/print_bar.Rd b/man/print_bar.Rd new file mode 100644 index 0000000..9738d57 --- /dev/null +++ b/man/print_bar.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.Utility.r +\name{print_bar} +\alias{print_bar} +\title{Print progress bar} +\usage{ +print_bar( + i, + n, + type = c("type1", "type3"), + symbol = "-", + tmp.file = NULL, + symbol.head = ">>>", + symbol.tail = ">", + fixed.points = TRUE, + points = seq(0, 100, 1), + symbol.len = 48, + verbose = TRUE +) +} +\arguments{ +\item{i}{the current loop number} + +\item{n}{max loop number} + +\item{type}{type1 for "for" function} + +\item{symbol}{the symbol for the rate of progress} + +\item{tmp.file}{the opened file of "fifo" function} + +\item{symbol.head}{the head for the bar} + +\item{symbol.tail}{the tail for the bar} + +\item{fixed.points}{whether use the setted points which will be printed} + +\item{points}{the setted points which will be printed} + +\item{symbol.len}{the total length of progress bar} +} +\description{ +Print progress bar +} +\keyword{internal} diff --git a/man/print_info.Rd b/man/print_info.Rd new file mode 100644 index 0000000..c4fbdfd --- /dev/null +++ b/man/print_info.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.Utility.r +\name{print_info} +\alias{print_info} +\title{Print R Package information, include title, short_title, logo, version, authors, contact} +\usage{ +print_info( + welcome = NULL, + title = NULL, + short_title = NULL, + logo = NULL, + version = NULL, + authors = NULL, + contact = NULL, + linechar = "=", + width = NULL, + verbose = TRUE +) +} +\arguments{ +\item{welcome}{welcome text, for example: "Welcom to "} + +\item{title}{long text to introduct package} + +\item{short_title}{short label, top-left of logo} + +\item{logo}{logo} + +\item{version}{short label, bottom-right of logo} + +\item{authors}{} + +\item{contact}{email or website} + +\item{width}{banner width} +} +\description{ +Build date: Oct 22, 2018 +Last update: Oct 22, 2018 +} +\author{ +Haohao Zhang +} +\keyword{internal} diff --git a/man/rule_wrap.Rd b/man/rule_wrap.Rd new file mode 100644 index 0000000..9e3c4bc --- /dev/null +++ b/man/rule_wrap.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MVP.Utility.r +\name{rule_wrap} +\alias{rule_wrap} +\title{wrap text to multiple line, align left, right or center.} +\usage{ +rule_wrap(string, width, align = "center", linechar = " ") +} +\description{ +Build date: Oct 22, 2018 +Last update: Dec 12, 2018 +by using base::strwrap. +} +\author{ +Haohao Zhang +} +\keyword{internal} diff --git a/src/Makevars b/src/Makevars index 9b030a3..eb00a2e 100755 --- a/src/Makevars +++ b/src/Makevars @@ -14,7 +14,7 @@ ## installation will do the right thing. Should you need it, uncomment it and ## set the appropriate value, possibly CXX17. -## CXX_STD = CXX11 +CXX_STD = CXX17 PKG_CXXFLAGS = -DARMA_64BIT_WORD=1 $(SHLIB_OPENMP_CXXFLAGS) PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/src/Makevars.win b/src/Makevars.win index b487221..e2b3f4e 100755 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -14,7 +14,7 @@ ## installation will do the right thing. Should you need it, uncomment it and ## set the appropriate value, possibly CXX17. -## CXX_STD = CXX11 +CXX_STD = CXX17 PKG_CXXFLAGS = -DARMA_64BIT_WORD=1 $(SHLIB_OPENMP_CXXFLAGS) -Wa,-mbig-obj PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index b19df08..464b264 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -243,7 +243,7 @@ BEGIN_RCPP END_RCPP } // kin_cal -SEXP kin_cal(SEXP pBigMat, const Nullable geno_ind, const Nullable marker_ind, const Nullable marker_freq, const bool marker_bycol, int threads, size_t step, bool mkl, bool verbose); +SEXP kin_cal(SEXP pBigMat, const Nullable geno_ind, const Nullable marker_ind, const Nullable marker_freq, const bool marker_bycol, int threads, int step, bool mkl, bool verbose); RcppExport SEXP _rMVP_kin_cal(SEXP pBigMatSEXP, SEXP geno_indSEXP, SEXP marker_indSEXP, SEXP marker_freqSEXP, SEXP marker_bycolSEXP, SEXP threadsSEXP, SEXP stepSEXP, SEXP mklSEXP, SEXP verboseSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -254,7 +254,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const Nullable >::type marker_freq(marker_freqSEXP); Rcpp::traits::input_parameter< const bool >::type marker_bycol(marker_bycolSEXP); Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); - Rcpp::traits::input_parameter< size_t >::type step(stepSEXP); + Rcpp::traits::input_parameter< int >::type step(stepSEXP); Rcpp::traits::input_parameter< bool >::type mkl(mklSEXP); Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP); rcpp_result_gen = Rcpp::wrap(kin_cal(pBigMat, geno_ind, marker_ind, marker_freq, marker_bycol, threads, step, mkl, verbose)); diff --git a/src/assoc.cpp b/src/assoc.cpp index 59e1c3f..1534dc0 100755 --- a/src/assoc.cpp +++ b/src/assoc.cpp @@ -8,7 +8,7 @@ using namespace std; using namespace Rcpp; using namespace arma; -arma::mat GInv(const arma::mat A){ +arma::mat GInv(const arma::mat& A){ arma::mat ginv; if(A.n_rows == 1){ @@ -37,7 +37,7 @@ arma::mat GInv(const arma::mat A){ } template -NumericVector getRow(XPtr pMat, const int row){ +NumericVector getRow(const XPtr pMat, const int row){ MatrixAccessor genomat = MatrixAccessor(*pMat); @@ -57,18 +57,7 @@ NumericVector getRow(SEXP pBigMat, const int row){ XPtr xpMat(pBigMat); - switch(xpMat->matrix_type()){ - case 1: - return getRow(xpMat, row); - case 2: - return getRow(xpMat, row); - case 4: - return getRow(xpMat, row); - case 8: - return getRow(xpMat, row); - default: - throw Rcpp::exception("unknown type detected for big.matrix object!"); - } + DISPATCH_MATRIX_TYPE(getRow, xpMat, row); } template @@ -102,16 +91,16 @@ SEXP glm_c(const arma::vec &y, const arma::mat &X, const arma::mat & iXX, XPtr #include #include +#include using namespace std; using namespace Rcpp; @@ -48,6 +49,23 @@ vector split_line(const string& str, const string& whitespace = " \t\r\n } } +// Zero-allocation variant for hot-path inner loops: views into the caller's string. +// The caller must ensure the underlying string outlives the returned views. +vector split_line_view(string_view str, string_view whitespace = " \t\r\n") +{ + vector result; + size_t col_end = 0; + while (true) { + size_t col_begin = str.find_first_not_of(whitespace, col_end); + if (col_begin == string_view::npos) + return result; + col_end = str.find_first_of(whitespace, col_begin); + if (col_end == string_view::npos) + col_end = str.size(); + result.push_back(str.substr(col_begin, col_end - col_begin)); + } +} + // ***** VCF ***** // [[Rcpp::export]] @@ -116,8 +134,8 @@ List vcf_parser_map(std::string vcf_file, std::string out) { } -double vcf_marker_parser(string m, double NA_C) { - if (('0' == m[0] || '1' == m[0]) && ('0' == m[2] || '1' == m[2])) { +double vcf_marker_parser(string_view m, double NA_C) { + if (m.size() >= 3 && ('0' == m[0] || '1' == m[0]) && ('0' == m[2] || '1' == m[2])) { return (m[0] - '0' + m[2] - '0'); } else { return NA_C; @@ -131,15 +149,15 @@ void vcf_parser_genotype(std::string vcf_file, XPtr pMat, long maxLin omp_setup(threads); string line; - vector l; + vector l; vector markers; size_t m; MatrixAccessor mat = MatrixAccessor(*pMat); - + // progress bar MinimalProgressBar_perc pb; Progress progress(pMat->ncol(), verbose, pb); - + // Skip Header string prefix("#CHROM"); bool have_header = false; @@ -166,7 +184,7 @@ void vcf_parser_genotype(std::string vcf_file, XPtr pMat, long maxLin } #pragma omp parallel for private(l, markers) for (std::size_t i = 0; i < buffer.size(); i++) { - l = split_line(buffer[i], "\t"); + l = split_line_view(buffer[i], "\t"); // boost::split(l, buffer[i], boost::is_any_of("\t")); // markers.clear(); // vector(l.begin() + 9, l.end()).swap(l); @@ -303,7 +321,7 @@ List hapmap_parser_map(std::string hmp_file, std::string out) { _["Major"] = Major); } -double hapmap_marker_parser(string m, char major, double NA_C) { +double hapmap_marker_parser(string_view m, char major, double NA_C) { if (m.length() == 1) { // Hapmap // Rcout << "major: " << major << '\t' << "now: " << m[0] << endl; if (m[0] == '+' || m[0] == '0' || m[0] == '-' || m[0] == 'N') { @@ -335,6 +353,7 @@ void hapmap_parser_genotype(std::string hmp_file, std::vector Major string line; char major; vector l; + vector lv; MatrixAccessor mat = MatrixAccessor(*pMat); // progress bar @@ -372,15 +391,15 @@ void hapmap_parser_genotype(std::string hmp_file, std::vector Major } size_t n_marker = buffer.size(); - #pragma omp parallel for private(l, major) + #pragma omp parallel for private(lv, major) for (size_t i = 0; i < n_marker; i++) { - l = split_line(buffer[i], " \t"); + lv = split_line_view(buffer[i], " \t"); // boost::split(l, buffer[i], boost::is_any_of(" \t")); - if(l.size() != n_col) + if(lv.size() != n_col) Rcpp::stop(("line " + to_string(m+i+2) + " does not have " + to_string(n_col) + " elements in HAPMAP file.").c_str()); major = Major[m + i][0]; - for(size_t j = 0; j < (l.size() - 11); j++) { - mat[m + i][j] = static_cast(hapmap_marker_parser(l[j + 11], major, NA_C)); + for(size_t j = 0; j < (lv.size() - 11); j++) { + mat[m + i][j] = static_cast(hapmap_marker_parser(lv[j + 11], major, NA_C)); } } progress.increment(n_marker); @@ -455,16 +474,18 @@ void write_bfile(XPtr pMat, std::string bed_file, double NA_C, bool m vector geno(n); MatrixAccessor mat = MatrixAccessor(*pMat); - FILE *fout; - fout = fopen(bed_file.c_str(), "wb"); - + std::ofstream fout(bed_file, std::ios::binary); + if (!fout.is_open()) { + Rcpp::stop("Cannot open file for writing: " + bed_file); + } + // progress bar MinimalProgressBar_perc pb; Progress progress(m, verbose, pb); - + // magic number of bfile const unsigned char magic_bytes[] = { 0x6c, 0x1b, 0x01 }; - fwrite((char*)magic_bytes, 1, 3, fout); + fout.write(reinterpret_cast(magic_bytes), 3); // map std::map code; @@ -485,7 +506,7 @@ void write_bfile(XPtr pMat, std::string bed_file, double NA_C, bool m } geno[j] = p; } - fwrite((char*)geno.data(), 1, geno.size(), fout); + fout.write(reinterpret_cast(geno.data()), geno.size()); progress.increment(); } }else{ @@ -499,11 +520,10 @@ void write_bfile(XPtr pMat, std::string bed_file, double NA_C, bool m } geno[j] = p; } - fwrite((char*)geno.data(), 1, geno.size(), fout); + fout.write(reinterpret_cast(geno.data()), geno.size()); progress.increment(); } } - fclose(fout); return; } @@ -538,7 +558,7 @@ void read_bfile(std::string bed_file, XPtr pMat, long maxLine, double long n = ind / 4; // 4 individual = 1 bit if (ind % 4 != 0) n++; - char *buffer; + std::vector buffer; long buffer_size; MatrixAccessor mat = MatrixAccessor(*pMat); @@ -550,11 +570,13 @@ void read_bfile(std::string bed_file, XPtr pMat, long maxLine, double code[0] = static_cast(2); // open file - FILE *fin; - fin = fopen(bed_file.c_str(), "rb"); - fseek(fin, 0, SEEK_END); - long length = ftell(fin); - rewind(fin); + std::ifstream fin(bed_file, std::ios::binary); + if (!fin.is_open()) { + Rcpp::stop("Cannot open .bed file: " + bed_file); + } + fin.seekg(0, std::ios::end); + long length = (long)fin.tellg(); + fin.seekg(0, std::ios::beg); // get buffer_size buffer_size = maxLine > 0 ? (maxLine * n) : (length - 3); @@ -566,10 +588,9 @@ void read_bfile(std::string bed_file, XPtr pMat, long maxLine, double Progress progress(n_block, verbose, pb); // magic number of bfile - buffer = new char [3]; - if (fread(buffer, 1, 3, fin) != 3) { - // file too short or corrupt - magic bytes unreadable - fclose(fin); + std::vector magic_buffer(3); + fin.read(magic_buffer.data(), 3); + if (fin.gcount() != 3) { Rcpp::stop("Cannot read magic bytes from .bed file."); } @@ -577,8 +598,9 @@ void read_bfile(std::string bed_file, XPtr pMat, long maxLine, double size_t cond; long block_start; for (int i = 0; i < n_block; i++) { - buffer = new char [buffer_size]; - size_t n_read = fread(buffer, 1, buffer_size, fin); + buffer.resize(buffer_size); + fin.read(buffer.data(), buffer_size); + size_t n_read = (size_t)fin.gcount(); if (n_read == 0) break; // end of file @@ -599,7 +621,6 @@ void read_bfile(std::string bed_file, XPtr pMat, long maxLine, double } progress.increment(); } - fclose(fin); return; } diff --git a/src/fit_diago.cpp b/src/fit_diago.cpp index c266a30..e854aee 100755 --- a/src/fit_diago.cpp +++ b/src/fit_diago.cpp @@ -370,13 +370,13 @@ MatrixXd conjugate_gradient(const MatrixXd& A, const VectorXd& b, int maxit, dou for (int i = 0; i < maxit; i++) { Ap = A * p; alpha = rsold / (p.transpose() * Ap); - x = x + (alpha * p.array()).matrix(); - r = r - (alpha * Ap.array()).matrix(); + x.noalias() += alpha * p; + r.noalias() -= alpha * Ap; rsnew = r.squaredNorm(); if (sqrt(rsnew) < tol) { break; } - p = r + ((rsnew / rsold) * p.array()).matrix(); + p = r + (rsnew / rsold) * p; rsold = rsnew; } return(x); @@ -419,8 +419,6 @@ SEXP geninv(SEXP GG) MatrixXd A(MatrixXd(mn, mn)); MatrixXd L(MatrixXd(mn, mn).setZero()); - - if (n < m) { transp = true; A = xxt(G); diff --git a/src/impute.cpp b/src/impute.cpp index 37d1979..8b76017 100755 --- a/src/impute.cpp +++ b/src/impute.cpp @@ -37,7 +37,7 @@ void impute_marker(XPtr pMat, bool mrkbycol = true, int threads=0, bo if(mrkbycol){ #pragma omp parallel for for (size_t i = 0; i < m; i++) { - std::vector na_index = {};; + std::vector na_index; size_t counts[3] = { 0 }; // count allele, record missing index @@ -62,7 +62,7 @@ void impute_marker(XPtr pMat, bool mrkbycol = true, int threads=0, bo }else{ #pragma omp parallel for for (size_t i = 0; i < m; i++) { - std::vector na_index = {};; + std::vector na_index; size_t counts[3] = { 0 }; // count allele, record missing index @@ -91,18 +91,7 @@ void impute_marker(XPtr pMat, bool mrkbycol = true, int threads=0, bo void impute_marker(SEXP pBigMat, bool mrkbycol = true, int threads=0, bool verbose=true) { XPtr xpMat(pBigMat); - switch(xpMat->matrix_type()) { - case 1: - return impute_marker(xpMat, mrkbycol, threads, verbose); - case 2: - return impute_marker(xpMat, mrkbycol, threads, verbose); - case 4: - return impute_marker(xpMat, mrkbycol, threads, verbose); - case 8: - return impute_marker(xpMat, mrkbycol, threads, verbose); - default: - throw Rcpp::exception("unknown type detected for big.matrix object!"); - } + DISPATCH_MATRIX_TYPE(impute_marker, xpMat, mrkbycol, threads, verbose); } template @@ -120,22 +109,18 @@ bool hasNA(XPtr pMat, bool mrkbycol = true, const Nullable(marker_ind) - 1; int m = _marker_ind.n_elem; if(mrkbycol){ - #pragma omp parallel for shared(HasNA) + #pragma omp parallel for reduction(||:HasNA) for (int j = 0; j < m; j++) { - if(HasNA) continue; for (int i = 0; i < n; i++) { - // if (mat[_marker_ind[j]][_geno_ind[i]] == NA_C) { if (isna(mat[_marker_ind[j]][_geno_ind[i]])) { HasNA = true; } } } }else{ - #pragma omp parallel for shared(HasNA) + #pragma omp parallel for reduction(||:HasNA) for (int j = 0; j < n; j++) { - if(HasNA) continue; for (int i = 0; i < m; i++) { - // if (mat[_geno_ind[j]][_marker_ind[i]] == NA_C) { if (isna(mat[_geno_ind[j]][_marker_ind[i]])) { HasNA = true; } @@ -144,22 +129,18 @@ bool hasNA(XPtr pMat, bool mrkbycol = true, const Nullablencol(); j++) { - if(HasNA) continue; for (int i = 0; i < n; i++) { - // if (mat[j][_geno_ind[i]] == NA_C) { if (isna(mat[j][_geno_ind[i]])) { HasNA = true; } } } }else{ - #pragma omp parallel for shared(HasNA) + #pragma omp parallel for reduction(||:HasNA) for (int j = 0; j < n; j++) { - if(HasNA) continue; for (int i = 0; i < pMat->nrow(); i++) { - // if (mat[_geno_ind[j]][i] == NA_C) { if (isna(mat[_geno_ind[j]][i])) { HasNA = true; } @@ -172,22 +153,18 @@ bool hasNA(XPtr pMat, bool mrkbycol = true, const Nullable(marker_ind) - 1; int m = _marker_ind.n_elem; if(mrkbycol){ - #pragma omp parallel for shared(HasNA) + #pragma omp parallel for reduction(||:HasNA) for (int j = 0; j < m; j++) { - if(HasNA) continue; for (int i = 0; i < pMat->nrow(); i++) { - // if (mat[_marker_ind[j]][i] == NA_C) { if (isna(mat[_marker_ind[j]][i])) { HasNA = true; } } } }else{ - #pragma omp parallel for shared(HasNA) + #pragma omp parallel for reduction(||:HasNA) for (int j = 0; j < pMat->ncol(); j++) { - if(HasNA) continue; for (int i = 0; i < m; i++) { - // if (mat[j][_marker_ind[i]] == NA_C) { if (isna(mat[j][_marker_ind[i]])) { HasNA = true; } @@ -195,11 +172,9 @@ bool hasNA(XPtr pMat, bool mrkbycol = true, const Nullablencol(); j++) { - if(HasNA) continue; for (int i = 0; i < pMat->nrow(); i++) { - // if (mat[j][i] == NA_C) { if (isna(mat[j][i])) { HasNA = true; } @@ -214,16 +189,5 @@ bool hasNA(XPtr pMat, bool mrkbycol = true, const Nullable geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const int threads=1) { XPtr xpMat(pBigMat); - switch(xpMat->matrix_type()) { - case 1: - return hasNA(xpMat, mrkbycol, geno_ind, marker_ind, threads); - case 2: - return hasNA(xpMat, mrkbycol, geno_ind, marker_ind, threads); - case 4: - return hasNA(xpMat, mrkbycol, geno_ind, marker_ind, threads); - case 8: - return hasNA(xpMat, mrkbycol, geno_ind, marker_ind, threads); - default: - throw Rcpp::exception("unknown type detected for big.matrix object!"); - } + DISPATCH_MATRIX_TYPE(hasNA, xpMat, mrkbycol, geno_ind, marker_ind, threads); } diff --git a/src/kinship.cpp b/src/kinship.cpp index 17bf1a4..b18c1cf 100755 --- a/src/kinship.cpp +++ b/src/kinship.cpp @@ -99,48 +99,6 @@ arma::vec BigRowMean(XPtr pMat, bool marker_bycol = true, size_t step } Z_buffer.reset(); - // 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 mean; } @@ -149,174 +107,10 @@ arma::vec BigRowMean(SEXP pBigMat, bool marker_bycol = true, size_t step = 10000 XPtr xpMat(pBigMat); - switch(xpMat->matrix_type()) { - case 1: - return BigRowMean(xpMat, marker_bycol, step, threads, geno_ind, verbose); - case 2: - return BigRowMean(xpMat, marker_bycol, step, threads, geno_ind, verbose); - case 4: - return BigRowMean(xpMat, marker_bycol, step, threads, geno_ind, verbose); - case 8: - return BigRowMean(xpMat, marker_bycol, step, threads, geno_ind, verbose); - default: - throw Rcpp::exception("unknown type detected for big.matrix object!"); - } + DISPATCH_MATRIX_TYPE(BigRowMean, xpMat, marker_bycol, step, threads, geno_ind, verbose); } -// 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; - -// #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; -// } -// } -// } - -// return Rcpp::wrap(kin); -// } - - -// // [[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 SEXP kin_cal(XPtr pMat, const Nullable geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const Nullable marker_freq = R_NilValue, const bool marker_bycol = true, int threads = 0, int step = 5000, bool mkl = false, bool verbose = true){ @@ -371,80 +165,13 @@ SEXP kin_cal(XPtr pMat, const Nullable geno_ind = R_NilVa { cnt++; } + if (cnt != step) { Z_buffer.set_size(cnt, n); } - if(_geno_ind.is_empty()){ - if(_marker_ind.is_empty()){ - if(marker_bycol){ - #pragma omp parallel for - for(int l = 0; l < cnt; l++){ - for(int k = 0; k < n; k++){ - Z_buffer(l, k) = (double)bigm[(i_marker + l)][k]; - } - } - }else{ - #pragma omp parallel for - for(int k = 0; k < n; k++){ - for(int l = 0; l < cnt; l++){ - Z_buffer(l, k) = (double)bigm[k][(i_marker + l)]; - } - } - } - }else{ - if(marker_bycol){ - #pragma omp parallel for - for(int l = 0; l < cnt; l++){ - for(int k = 0; k < n; k++){ - Z_buffer(l, k) = (double)bigm[_marker_ind[(i_marker + l)]][k]; - } - } - }else{ - #pragma omp parallel for - for(int k = 0; k < n; k++){ - for(int l = 0; l < cnt; l++){ - Z_buffer(l, k) = (double)bigm[k][_marker_ind[(i_marker + l)]]; - } - } - } - } - }else{ - if(_marker_ind.is_empty()){ - if(marker_bycol){ - #pragma omp parallel for - for(int l = 0; l < cnt; l++){ - for(int k = 0; k < n; k++){ - Z_buffer(l, k) = (double)bigm[(i_marker + l)][_geno_ind[k]]; - } - } - }else{ - #pragma omp parallel for - for(int k = 0; k < n; k++){ - for(int l = 0; l < cnt; l++){ - Z_buffer(l, k) = (double)bigm[_geno_ind[k]][(i_marker + l)]; - } - } - } - }else{ - if(marker_bycol){ - #pragma omp parallel for - for(int l = 0; l < cnt; l++){ - for(int k = 0; k < n; k++){ - Z_buffer(l, k) = (double)bigm[_marker_ind[(i_marker + l)]][_geno_ind[k]]; - } - } - }else{ - #pragma omp parallel for - for(int k = 0; k < n; k++){ - for(int l = 0; l < cnt; l++){ - Z_buffer(l, k) = (double)bigm[_geno_ind[k]][_marker_ind[(i_marker + l)]]; - } - } - } - } - } - + fill_geno_buffer(Z_buffer, bigm, cnt, i_marker, _geno_ind, _marker_ind, marker_bycol, false); + if(marker_freq.isNotNull()){ Z_buffer.each_col() -= means.subvec(i_marker, i_marker + cnt - 1); }else{ @@ -486,20 +213,9 @@ SEXP kin_cal(XPtr pMat, const Nullable geno_ind = R_NilVa } // [[Rcpp::export]] -SEXP kin_cal(SEXP pBigMat, const Nullable geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const Nullable marker_freq = R_NilValue, const bool marker_bycol = true, int threads = 0, size_t step = 10000, bool mkl = false, bool verbose = true){ +SEXP kin_cal(SEXP pBigMat, const Nullable geno_ind = R_NilValue, const Nullable marker_ind = R_NilValue, const Nullable marker_freq = R_NilValue, const bool marker_bycol = true, int threads = 0, int step = 10000, bool mkl = false, bool verbose = true){ XPtr xpMat(pBigMat); - switch(xpMat->matrix_type()){ - case 1: - return kin_cal(xpMat, geno_ind, marker_ind, marker_freq, marker_bycol, threads, step, mkl, verbose); - case 2: - return kin_cal(xpMat, geno_ind, marker_ind, marker_freq, marker_bycol, threads, step, mkl, verbose); - case 4: - return kin_cal(xpMat, geno_ind, marker_ind, marker_freq, marker_bycol, threads, step, mkl, verbose); - case 8: - return kin_cal(xpMat, geno_ind, marker_ind, marker_freq, marker_bycol, threads, step, mkl, verbose); - default: - throw Rcpp::exception("unknown type detected for big.matrix object!"); - } + DISPATCH_MATRIX_TYPE(kin_cal, xpMat, geno_ind, marker_ind, marker_freq, marker_bycol, threads, step, mkl, verbose); } diff --git a/src/mvp_omp.h b/src/mvp_omp.h index de1d07f..043f8c0 100755 --- a/src/mvp_omp.h +++ b/src/mvp_omp.h @@ -21,18 +21,17 @@ #include -static int omp_setup(int threads); static inline int omp_setup(int threads=0) { int t = 1; #ifdef _OPENMP if (threads == 0) { + // Reserve one thread for OS; default to physical processor count - 1 t = omp_get_num_procs() - 1; t = t > 0 ? t : 1; } else { t = threads > 0 ? threads : 1; } omp_set_num_threads(t); -#else #endif return t; } diff --git a/src/rMVP.h b/src/rMVP.h index b988697..0a8847c 100644 --- a/src/rMVP.h +++ b/src/rMVP.h @@ -1,6 +1,5 @@ -#ifndef RMVP_ -#define RMVP_ -#endif +#ifndef RMVP_H_ +#define RMVP_H_ #if !defined(ARMA_64BIT_WORD) #define ARMA_64BIT_WORD 1 @@ -18,306 +17,261 @@ #include "progress_bar.hpp" #include "mvp_omp.h" -class MinimalProgressBar_plus: public ProgressBar{ - public: - MinimalProgressBar_plus() { - _finalized = false; - } - - ~MinimalProgressBar_plus() {} - - std::string _time_to_string(double seconds, float progress) { - - int time = (int) seconds; - - int hour = 0; - int min = 0; - int sec = 0; - - hour = time / 3600; - time = time % 3600; - min = time / 60; - time = time % 60; - sec = time; - - std::stringstream time_strs; - time_strs << (progress < 1.0 ? "TimeLeft: " : "RunTime: "); - if (hour != 0) time_strs << hour << "h"; - if (hour != 0 || min != 0) time_strs << min << "m"; - time_strs << sec << "s"; - std::string time_str = time_strs.str(); - - return time_str; - } - - int _compute_nb_ticks(float progress) { - return int(progress * _max_ticks); - } - - std::string _construct_ticks_display_string(int nb) { - std::stringstream ticks_strs; - for (int i = 1; i <= _max_ticks; ++i) { - if(i < 4){ - ticks_strs << ">"; - } else if (i < nb) { - ticks_strs << "-"; - } else if(i == nb) { - ticks_strs << ">"; - } else { - ticks_strs << " "; - } - } - std::string tick_space_string = ticks_strs.str(); - return tick_space_string; - } - - void end_display() { - update(1); +// --------------------------------------------------------------------------- +// MinimalProgressBarBase — shared implementation for timed progress bars. +// Concrete subclasses must implement display() and _construct_ticks_display_string(). +// --------------------------------------------------------------------------- +class MinimalProgressBarBase : public ProgressBar { +public: + explicit MinimalProgressBarBase(int max_ticks) + : _max_ticks(max_ticks), empty_length_p(0), _finalized(false), + _timer_flag(true), _ticks_displayed(0) {} + + virtual ~MinimalProgressBarBase() {} + + virtual void display() = 0; + + void end_display() { update(1); } + + void update(float progress) { + if (_finalized) return; + + if (_timer_flag) { + _timer_flag = false; + time(&start); + } else { + int nb_ticks = _compute_nb_ticks(progress); + int delta = nb_ticks - _ticks_displayed; + if (delta > 0) { + _ticks_displayed = nb_ticks; + std::string cur_display = _construct_ticks_display_string(nb_ticks); + + time(&end); + double pas_time = std::difftime(end, start); + double rem_time = (progress < 1.0 ? (pas_time / progress) * (1 - progress) : pas_time); + if (rem_time < 1 && rem_time > 0.5) rem_time = 1; + + std::string time_string = _time_to_string(rem_time, progress); + int empty_length = time_string.length(); + std::string empty_space; + + std::stringstream strs; + if (empty_length_p && abs(empty_length - empty_length_p)) { + empty_space = std::string(abs(empty_length - empty_length_p), ' '); + strs << "[" << cur_display << "] " << time_string << empty_space; + } else { + strs << "[" << cur_display << "] " << time_string; + } + empty_length_p = empty_length; + + std::string temp_str = strs.str(); + char const* char_type = temp_str.c_str(); + REprintf("\r"); + REprintf("%s", char_type); + } + if (_ticks_displayed >= _max_ticks) _finalize_display(); + } + } + + void flush_console() { + #if !defined(WIN32) && !defined(__WIN32) && !defined(__WIN32__) + R_FlushConsole(); + #endif } +protected: + int _max_ticks; + int empty_length_p; + bool _finalized; + bool _timer_flag; + time_t start, end; + int _ticks_displayed; + void _finalize_display() { - if (_finalized) return; - REprintf("\n"); - _finalized = true; + if (_finalized) return; + REprintf("\n"); + _finalized = true; } - void display() { - flush_console(); - } - - void update(float progress) { - - // stop if already finalized - if (_finalized) return; - - // start time measurement when update() is called the first time - if (_timer_flag) { - _timer_flag = false; - // measure start time - time(&start); - } else { - - int nb_ticks = _compute_nb_ticks(progress); - int delta = nb_ticks - _ticks_displayed; - if (delta > 0) { - _ticks_displayed = nb_ticks; - std::string cur_display = _construct_ticks_display_string(nb_ticks); - - // measure current time - time(&end); - - // calculate passed time and remaining time (in seconds) - double pas_time = std::difftime(end, start); - double rem_time = (progress < 1.0 ? (pas_time / progress) * (1 - progress) : pas_time); - if(rem_time < 1 && rem_time > 0.5) rem_time = 1; - - // convert seconds to time string - std::string time_string = _time_to_string(rem_time, progress); - - // ensure overwriting of old time info - int empty_length = time_string.length(); - - std::string empty_space; - - // merge progress bar and time string - std::stringstream strs; - if(empty_length_p && abs(empty_length - empty_length_p)){ - empty_space = std::string(abs(empty_length - empty_length_p), ' '); - strs << "[" << cur_display << "] " << time_string << empty_space; - }else{ - strs << "[" << cur_display << "] " << time_string; - } - empty_length_p = empty_length; - // strs << "[" << cur_display << "]"; - - std::string temp_str = strs.str(); - char const* char_type = temp_str.c_str(); - - // print: remove old display and replace with new - REprintf("\r"); - REprintf("%s", char_type); - - } - if (_ticks_displayed >= _max_ticks) - // end_display(); - _finalize_display(); - } - } - - void flush_console() { - #if !defined(WIN32) && !defined(__WIN32) && !defined(__WIN32__) - R_FlushConsole(); - #endif + std::string _time_to_string(double seconds, float progress) { + int time = (int) seconds; + int hour = 0, min = 0, sec = 0; + hour = time / 3600; + time = time % 3600; + min = time / 60; + time = time % 60; + sec = time; + std::stringstream time_strs; + time_strs << (progress < 1.0 ? "TimeLeft: " : "RunTime: "); + if (hour != 0) time_strs << hour << "h"; + if (hour != 0 || min != 0) time_strs << min << "m"; + time_strs << sec << "s"; + return time_strs.str(); } - private: - int empty_length_p = 0; - int _max_ticks = 45; - bool _finalized; - bool _timer_flag = true; - time_t start, end; - int _ticks_displayed = 0; + int _compute_nb_ticks(float progress) { + return int(progress * _max_ticks); + } + + virtual std::string _construct_ticks_display_string(int nb) = 0; }; -class MinimalProgressBar_perc: public ProgressBar{ - public: - MinimalProgressBar_perc() { - _finalized = false; - } - - ~MinimalProgressBar_perc() {} - - std::string _time_to_string(double seconds, float progress) { - - int time = (int) seconds; - - int hour = 0; - int min = 0; - int sec = 0; - - hour = time / 3600; - time = time % 3600; - min = time / 60; - time = time % 60; - sec = time; - - std::stringstream time_strs; - time_strs << (progress < 1.0 ? "TimeLeft: " : "RunTime: "); - if (hour != 0) time_strs << hour << "h"; - if (hour != 0 || min != 0) time_strs << min << "m"; - time_strs << sec << "s"; - std::string time_str = time_strs.str(); - - return time_str; - } - - int _compute_nb_ticks(float progress) { - return int(progress * _max_ticks); - } - - std::string _construct_ticks_display_string(int nb) { - std::stringstream ticks_strs; - for (int i = 1; i <= _max_ticks; ++i) { - if (i <= nb) { - ticks_strs << "*"; - } else { - ticks_strs << " "; - } - } - std::string tick_space_string = ticks_strs.str(); - return tick_space_string; - } - - void end_display() { - update(1); +// --------------------------------------------------------------------------- +// MinimalProgressBar_plus — arrow-style progress bar [>>>-----> ] ... +// --------------------------------------------------------------------------- +class MinimalProgressBar_plus : public MinimalProgressBarBase { +public: + MinimalProgressBar_plus() : MinimalProgressBarBase(45) {} + ~MinimalProgressBar_plus() {} + + void display() override { flush_console(); } + +protected: + std::string _construct_ticks_display_string(int nb) override { + std::stringstream ticks_strs; + for (int i = 1; i <= _max_ticks; ++i) { + if (i < 4) { + ticks_strs << ">"; + } else if (i < nb) { + ticks_strs << "-"; + } else if (i == nb) { + ticks_strs << ">"; + } else { + ticks_strs << " "; + } + } + return ticks_strs.str(); } +}; - void _finalize_display() { - if (_finalized) return; - REprintf("\n"); - _finalized = true; +// --------------------------------------------------------------------------- +// MinimalProgressBar_perc — asterisk-style bar with percentage ruler header +// --------------------------------------------------------------------------- +class MinimalProgressBar_perc : public MinimalProgressBarBase { +public: + MinimalProgressBar_perc() : MinimalProgressBarBase(49) {} + ~MinimalProgressBar_perc() {} + + void display() override { + REprintf("0%% 10 20 30 40 50 60 70 80 90 100%%\n"); + REprintf("[----|----|----|----|----|----|----|----|----|----|\n"); + flush_console(); } - void display() { - REprintf("0%% 10 20 30 40 50 60 70 80 90 100%%\n"); - REprintf("[----|----|----|----|----|----|----|----|----|----|\n"); - flush_console(); - } - - void update(float progress) { - - // stop if already finalized - if (_finalized) return; - - // start time measurement when update() is called the first time - if (_timer_flag) { - _timer_flag = false; - // measure start time - time(&start); - } else { - - int nb_ticks = _compute_nb_ticks(progress); - int delta = nb_ticks - _ticks_displayed; - if (delta > 0) { - _ticks_displayed = nb_ticks; - std::string cur_display = _construct_ticks_display_string(nb_ticks); - - // measure current time - time(&end); - - // calculate passed time and remaining time (in seconds) - double pas_time = std::difftime(end, start); - double rem_time = (progress < 1.0 ? (pas_time / progress) * (1 - progress) : pas_time); - if(rem_time < 1 && rem_time > 0.5) rem_time = 1; - - // convert seconds to time string - std::string time_string = _time_to_string(rem_time, progress); - - // ensure overwriting of old time info - int empty_length = time_string.length(); - - std::string empty_space; - - // merge progress bar and time string - std::stringstream strs; - if(empty_length_p && abs(empty_length - empty_length_p)){ - empty_space = std::string(abs(empty_length - empty_length_p), ' '); - strs << "[" << cur_display << "] " << time_string << empty_space; - }else{ - strs << "[" << cur_display << "] " << time_string; - } - empty_length_p = empty_length; - // strs << "[" << cur_display << "]"; - - std::string temp_str = strs.str(); - char const* char_type = temp_str.c_str(); - - // print: remove old display and replace with new - REprintf("\r"); - REprintf("%s", char_type); - - } - if (_ticks_displayed >= _max_ticks) - // end_display(); - _finalize_display(); - } - } - - void flush_console() { - #if !defined(WIN32) && !defined(__WIN32) && !defined(__WIN32__) - R_FlushConsole(); - #endif +protected: + std::string _construct_ticks_display_string(int nb) override { + std::stringstream ticks_strs; + for (int i = 1; i <= _max_ticks; ++i) { + if (i <= nb) { + ticks_strs << "*"; + } else { + ticks_strs << " "; + } + } + return ticks_strs.str(); + } +}; + +// --------------------------------------------------------------------------- +// MinimalProgressBar — simple inline percentage display (no ticks) +// --------------------------------------------------------------------------- +class MinimalProgressBar : public ProgressBar { +public: + MinimalProgressBar() { _finalized = false; } + ~MinimalProgressBar() {} + + void display() {} + + void update(float progress) { + if (_finalized) return; + REprintf("\r"); + REprintf("Calculating in process...(finished %.2f%%)", progress * 100); + } + + void end_display() { + if (_finalized) return; + REprintf("\r"); + REprintf("Calculating in process...(finished 100.00%%)"); + REprintf("\n"); + _finalized = true; } - private: - int empty_length_p = 0; - int _max_ticks = 49; +private: bool _finalized; - bool _timer_flag = true; - time_t start, end; - int _ticks_displayed = 0; }; -class MinimalProgressBar: public ProgressBar{ - public: - MinimalProgressBar() { - _finalized = false; - } - ~MinimalProgressBar() {} - void display() {} - void update(float progress) { - if (_finalized) return; - REprintf("\r"); - REprintf("Calculating in process...(finished %.2f%%)", progress * 100); - } - void end_display() { - if (_finalized) return; - REprintf("\r"); - - REprintf("Calculating in process...(finished 100.00%%)"); - REprintf("\n"); - _finalized = true; - } - private: - bool _finalized; -}; +// --------------------------------------------------------------------------- +// fill_geno_buffer — load a batch of cnt markers from a BigMatrix into buf. +// +// buf layout controlled by ind_in_rows: +// true → buf is n×cnt, buf(k,l) = genotype of individual k at marker (i_marker+l) +// false → buf is cnt×n, buf(l,k) = genotype of individual k at marker (i_marker+l) +// +// Empty geno_ind → use direct row/col index k. +// Empty marker_ind → use sequential index (i_marker + l). +// +// NOTE: Does NOT handle NA values. Callers requiring NA-aware buffer filling +// (e.g. BigRowMean) must use specialized logic. +// --------------------------------------------------------------------------- +template +void fill_geno_buffer(arma::mat& buf, + MatrixAccessor& mat, + int cnt, + int i_marker, + const arma::uvec& geno_ind, + const arma::uvec& marker_ind, + bool marker_bycol, + bool ind_in_rows) { + int n = ind_in_rows ? (int)buf.n_rows : (int)buf.n_cols; + bool has_geno_ind = !geno_ind.is_empty(); + bool has_mrk_ind = !marker_ind.is_empty(); + + if (marker_bycol) { + if (ind_in_rows) { + #pragma omp parallel for + for (int l = 0; l < cnt; l++) { + int col = has_mrk_ind ? (int)marker_ind[i_marker + l] : (i_marker + l); + for (int k = 0; k < n; k++) + buf(k, l) = (double)mat[col][has_geno_ind ? (int)geno_ind[k] : k]; + } + } else { + #pragma omp parallel for + for (int l = 0; l < cnt; l++) { + int col = has_mrk_ind ? (int)marker_ind[i_marker + l] : (i_marker + l); + for (int k = 0; k < n; k++) + buf(l, k) = (double)mat[col][has_geno_ind ? (int)geno_ind[k] : k]; + } + } + } else { + if (ind_in_rows) { + #pragma omp parallel for + for (int k = 0; k < n; k++) { + int row = has_geno_ind ? (int)geno_ind[k] : k; + for (int l = 0; l < cnt; l++) + buf(k, l) = (double)mat[row][has_mrk_ind ? (int)marker_ind[i_marker + l] : (i_marker + l)]; + } + } else { + #pragma omp parallel for + for (int k = 0; k < n; k++) { + int row = has_geno_ind ? (int)geno_ind[k] : k; + for (int l = 0; l < cnt; l++) + buf(l, k) = (double)mat[row][has_mrk_ind ? (int)marker_ind[i_marker + l] : (i_marker + l)]; + } + } + } +} + +// Dispatch a bigmatrix function call based on element type. +// Requires xpMat to be the first argument to FUNC; trailing args forwarded via __VA_ARGS__. +// Usage: DISPATCH_MATRIX_TYPE(myFunc, xpMat, arg1, arg2, ...) +// Note: __VA_ARGS__ is intentionally never empty at any call site, so no ## GNU extension needed. +#define DISPATCH_MATRIX_TYPE(FUNC, xpMat, ...) \ + switch((xpMat)->matrix_type()) { \ + case 1: return FUNC (xpMat, __VA_ARGS__); \ + case 2: return FUNC (xpMat, __VA_ARGS__); \ + case 4: return FUNC (xpMat, __VA_ARGS__); \ + case 8: return FUNC(xpMat, __VA_ARGS__); \ + default: throw Rcpp::exception("unsupported bigmatrix type"); \ + } + +#endif // RMVP_H_ diff --git a/tests/testthat/test_regression_gwas.R b/tests/testthat/test_regression_gwas.R new file mode 100644 index 0000000..dfb8fd8 --- /dev/null +++ b/tests/testthat/test_regression_gwas.R @@ -0,0 +1,231 @@ +# Sprint 0.1: Numerical regression baseline for GWAS methods. +# Locks the exact numerical outputs of GLM, MLM (BRENT/EMMA/HE VC), +# FaSTLMM LL, and FarmCPU before refactoring begins. +# +# Expected values were captured on the extdata 15-marker dataset using +# rMVP 1.4.6 on x86-64 Linux with OpenBLAS. Tolerance is 1e-6 to accommodate +# cross-platform BLAS differences (Apple Accelerate, MKL, OpenBLAS agree to +# ~1e-7 on these inputs; 1e-6 gives one decade of headroom). +# +# To update after an intentional numerical change (documented in CHANGELOG): +# 1. Run each function interactively on the test dataset. +# 2. Print values with options(digits=15) and update the constants below. + +# Helper: load and cache test dataset once per test session. +.load_gwas_testdata <- local({ + cached <- NULL + function() { + if (!is.null(cached)) return(cached) + phe <- read.table( + system.file("extdata", "07_other", "mvp.phe", package = "rMVP"), + header = TRUE + ) + geno <- attach.big.matrix( + system.file("extdata", "06_mvp-impute", "mvp.imp.geno.desc", package = "rMVP") + ) + idx <- !is.na(phe[, 2]) + phe_clean <- phe[idx, ] + K <- MVP.K.VanRaden(geno, cpu = 1, verbose = FALSE) + eigenK <- eigen(K, symmetric = TRUE) + cached <<- list(phe = phe_clean, geno = geno, K = K, eigenK = eigenK) + cached + } +}) + +# Rows 8 and 9 are monomorphic in the extdata set — expected to produce NA/NaN. +# Rows checked: 1, 4, 11 (most significant by p-value), 15. +.CHECK_ROWS <- c(1L, 4L, 11L, 15L) + +# ---- GLM ---------------------------------------------------------------- + +test_that("GLM produces numerically stable results", { + skip_on_cran() + d <- .load_gwas_testdata() + result <- MVP.GLM(phe = d$phe, geno = d$geno, cpu = 1, verbose = FALSE) + + expect_true(is.matrix(result)) + expect_equal(nrow(result), 15L) + expect_equal(ncol(result), 3L) + + # Monomorphic markers return NA / NaN. + expect_true(is.na(result[8, 1])) + expect_true(is.na(result[9, 1])) + expect_true(is.nan(result[8, 3])) + expect_true(is.nan(result[9, 3])) + + # Row 11 is the most significant marker. + expect_equal(which.min(result[, 3]), 11L) + + # Beta (col 1), SE (col 2), p-value (col 3) at fixed rows. + tol <- 1e-6 + expect_equal(result[.CHECK_ROWS, 1], + c( 2.5728571, -2.8711111, -4.0233333, 2.8900000), tolerance = tol) + expect_equal(result[.CHECK_ROWS, 2], + c( 1.7134863, 2.7840715, 1.1267864, 2.7816275), tolerance = tol) + expect_equal(result[.CHECK_ROWS, 3], + c(0.17161487, 0.33258383, 0.00728702, 0.32920268), tolerance = tol) +}) + +# ---- Variance component estimation -------------------------------------- + +test_that("BRENT VC estimation is numerically stable", { + skip_on_cran() + d <- .load_gwas_testdata() + vc <- MVP.BRENT.Vg.Ve( + y = as.numeric(d$phe[, 2]), + X = matrix(1, nrow(d$phe)), + eigenK = d$eigenK, + verbose = FALSE + ) + + tol <- 1e-6 + expect_equal(vc$vg, 19.460680, tolerance = tol) + expect_equal(vc$ve, 8.090278, tolerance = tol) + expect_equal(vc$delta, 0.415724, tolerance = tol) +}) + +test_that("EMMA VC estimation is numerically stable", { + skip_on_cran() + d <- .load_gwas_testdata() + vc <- MVP.EMMA.Vg.Ve( + y = as.numeric(d$phe[, 2]), + X = matrix(1, nrow(d$phe)), + K = d$K + ) + + tol <- 1e-6 + expect_equal(vc$vg, 19.461379, tolerance = tol) + expect_equal(vc$ve, 8.090335, tolerance = tol) + expect_equal(vc$delta, 0.415712, tolerance = tol) + expect_equal(vc$REML, -27.102099, tolerance = tol) +}) + +test_that("HE VC estimation is numerically stable", { + skip_on_cran() + d <- .load_gwas_testdata() + vc <- MVP.HE.Vg.Ve( + y = as.numeric(d$phe[, 2]), + X = matrix(1, nrow(d$phe)), + K = d$K + ) + + tol <- 1e-6 + expect_equal(vc$vg, 8.288881, tolerance = tol) + expect_equal(vc$ve, 18.890817, tolerance = tol) + expect_equal(vc$delta, 2.279055, tolerance = tol) +}) + +test_that("FaSTLMM LL estimation is numerically stable", { + skip_on_cran() + d <- .load_gwas_testdata() + result <- MVP.FaSTLMM.LL( + pheno = as.matrix(d$phe[, 2]), + snp.pool = d$geno[], + X0 = matrix(1, nrow(d$phe)), + ncpus = 1 + ) + + tol <- 1e-6 + # delta hits the optimizer upper boundary (essentially Inf) on this small + # dataset, so vg collapses to ~0. Only beta, LL, and ve are meaningful. + expect_equal(result$beta, 100.298, tolerance = tol) + expect_equal(result$LL, -30.341552, tolerance = tol) + expect_equal(result$ve, 25.290616, tolerance = tol) +}) + +# ---- MLM ---------------------------------------------------------------- + +test_that("MLM (BRENT VC) produces numerically stable results", { + skip_on_cran() + d <- .load_gwas_testdata() + result <- MVP.MLM( + phe = d$phe, + geno = d$geno, + eigenK = d$eigenK, + vc.method = "BRENT", + cpu = 1, + verbose = FALSE + ) + + expect_true(is.matrix(result)) + expect_equal(nrow(result), 15L) + + # Monomorphic markers return NaN. + expect_true(is.nan(result[8, 1])) + expect_true(is.nan(result[9, 1])) + + tol <- 1e-6 + expect_equal(result[.CHECK_ROWS, 1], + c( 3.9431337, -0.0500524, -3.6286984, 2.4695755), tolerance = tol) + expect_equal(result[.CHECK_ROWS, 2], + c( 2.2801950, 2.7330146, 1.9299914, 2.7095227), tolerance = tol) + expect_equal(result[.CHECK_ROWS, 3], + c(0.12201210, 0.98583689, 0.09688123, 0.38870037), tolerance = tol) +}) + +test_that("MLM (EMMA VC) produces numerically stable results", { + skip_on_cran() + d <- .load_gwas_testdata() + result <- MVP.MLM( + phe = d$phe, + geno = d$geno, + K = d$K, + eigenK = d$eigenK, + vc.method = "EMMA", + cpu = 1, + verbose = FALSE + ) + + expect_true(is.matrix(result)) + expect_equal(nrow(result), 15L) + + expect_true(is.nan(result[8, 1])) + expect_true(is.nan(result[9, 1])) + + tol <- 1e-6 + expect_equal(result[.CHECK_ROWS, 1], + c( 3.9431501, -0.0500181, -3.6286904, 2.4695614), tolerance = tol) + expect_equal(result[.CHECK_ROWS, 2], + c( 2.2802260, 2.7330490, 1.9300190, 2.7095560), tolerance = tol) + expect_equal(result[.CHECK_ROWS, 3], + c(0.12201518, 0.98584678, 0.09688584, 0.38870857), tolerance = tol) +}) + +test_that("MLM (HE VC) produces numerically stable results", { + skip_on_cran() + d <- .load_gwas_testdata() + result <- MVP.MLM( + phe = d$phe, + geno = d$geno, + K = d$K, + eigenK = d$eigenK, + vc.method = "HE", + cpu = 1, + verbose = FALSE + ) + + expect_true(is.matrix(result)) + expect_equal(nrow(result), 15L) + + expect_true(is.nan(result[8, 1])) + expect_true(is.nan(result[9, 1])) + + tol <- 1e-6 + expect_equal(result[.CHECK_ROWS, 1], + c( 3.1813154, -1.8035841, -3.9345467, 2.9160929), tolerance = tol) + expect_equal(result[.CHECK_ROWS, 2], + c( 2.0865350, 2.7990780, 1.8065600, 2.8072150), tolerance = tol) + expect_equal(result[.CHECK_ROWS, 3], + c(0.16584474, 0.53738173, 0.06105992, 0.32927922), tolerance = tol) +}) + +# ---- FarmCPU ------------------------------------------------------------ + +test_that("FarmCPU produces numerically stable results", { + skip_on_cran() + # FarmCPU's bin-optimization algorithm requires a minimum number of markers + # to select pseudo-QTN candidates. The 15-marker extdata test set is too + # small. FarmCPU regression is covered by manual integration tests with the + # pig60K dataset (devscripts/benchmark_farmcpu.R). + skip("FarmCPU binning requires more markers than the 15-marker extdata set provides.") +})