refactor: production-grade code quality pass — dedup, vectorize, C++17, RAII, safety

.gitignore

@@ -46,3 +46,6 @@ vignettes/*.pdf
 *.knit.md
 .Rproj.user
 packages
+
+# R Files for comparisions
+compare_perf.R

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<T> 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.

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 <xll198708@gmail.com>
 Author: Lilin Yin [aut],
   Haohao Zhang [aut],

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]]

R/MVP.Constants.R

@@ -0,0 +1,25 @@
+# Package-level constants for rMVP statistical methods.
+# Naming convention: .<METHOD>_<DESCRIPTION> (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

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) {

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)