From 76a803a18ca6f1a19214a9c88b6a0a2da00ceb07 Mon Sep 17 00:00:00 2001 From: GAYNOR Chris Date: Mon, 19 Nov 2018 09:45:33 +0000 Subject: [PATCH] Added EM algorithm with 2 random effects --- DESCRIPTION | 14 +-- NAMESPACE | 2 + R/AlphaMME.R | 9 +- R/RcppExports.R | 53 +++++++++- man/AlphaMME-package.Rd | 9 -- man/AlphaMME.Rd | 13 +++ man/readMat.Rd | 7 +- man/solveMKM.Rd | 4 +- man/solveRRBLUP_EM.Rd | 31 ++++++ man/solveRRBLUP_EM2.Rd | 35 +++++++ src/MME.cpp | 213 +++++++++++++++++++++++++++++++++++++--- src/RcppExports.cpp | 45 ++++++++- src/init.c | 34 ++++--- src/optimize.cpp | 32 +++--- 14 files changed, 428 insertions(+), 73 deletions(-) delete mode 100644 man/AlphaMME-package.Rd create mode 100644 man/AlphaMME.Rd create mode 100644 man/solveRRBLUP_EM.Rd create mode 100644 man/solveRRBLUP_EM2.Rd diff --git a/DESCRIPTION b/DESCRIPTION index a299e80..ff9cfc7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,13 +1,15 @@ Package: AlphaMME Type: Package Title: Alpha Mixed Model Equations -Version: 0.1.0 -Date: 2018-01-25 -Authors@R: c(person("Chris","Gaynor",email="gaynor.robert@hotmail.com",role="cre")) -Description: This package contains mixed model equation solvers written in C++. +Version: 0.1.1 +Date: 2018-11-04 +Authors@R: c(person("Chris","Gaynor",email="gaynor.robert@hotmail.com",role=c("aut","cre"))) +Description: Fast mixed model solvers for genomic selection. All solvers are written in C++ + and based on solvers found in other R packages. License: MIT + file LICENSE URL: https://bitbucket.org/hickeyjohnteam/alphamme -Imports: Rcpp (>= 0.12.13), RcppArmadillo +Imports: Rcpp (>= 0.12.13) LinkingTo: Rcpp, RcppArmadillo -RoxygenNote: 6.0.1 +Encoding: UTF-8 +RoxygenNote: 6.1.0 NeedsCompilation: true diff --git a/NAMESPACE b/NAMESPACE index 1b8094d..0322f87 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,6 +13,8 @@ export(solveMVM) export(solveRRBLUP) export(solveRRBLUPMK) export(solveRRBLUPMV) +export(solveRRBLUP_EM) +export(solveRRBLUP_EM2) export(solveUVM) import(Rcpp) import(RcppArmadillo) diff --git a/R/AlphaMME.R b/R/AlphaMME.R index 90f5daa..13658f8 100644 --- a/R/AlphaMME.R +++ b/R/AlphaMME.R @@ -4,8 +4,11 @@ #' @title Alpha Mixed Model Equations #' #' @description -#' This package contains mixed model equation solvers written in C++. +#' This package is designed for genomic selection. It contains +#' several mixed model equation solvers from other R packages. +#' The solvers have been rewritten in C++ and streamlined for +#' increased speed and more efficient memory usage. #' #' @docType package -#' @name AlphaMME-package -NULL \ No newline at end of file +#' @name AlphaMME +NULL diff --git a/R/RcppExports.R b/R/RcppExports.R index 1499569..6900acb 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -8,10 +8,10 @@ #' file. Requires knowledge of the number of rows #' and columns in the file. #' -#' @param fileName path to the file to read +#' @param fileName path to the file being read #' @param rows number of rows to read in #' @param cols number of columns to read in -#' @param sep a single character seperating data entries +#' @param sep a single character delimiter seperating data entries #' @param skipRows number of rows to skip #' @param skipCols number of columns to skip #' @@ -109,10 +109,11 @@ solveRRBLUPMV <- function(Y, X, M, tol = 1e-6, maxIter = 1000L) { #' @param Zlist a list of Z matrices #' @param Klist a list of K matrices #' @param maxIter maximum number of iteration +#' @param tol tolerance for convergence #' #' @export -solveMKM <- function(y, X, Zlist, Klist, maxIter = 40L) { - .Call(`_AlphaMME_solveMKM`, y, X, Zlist, Klist, maxIter) +solveMKM <- function(y, X, Zlist, Klist, maxIter = 40L, tol = 1e-4) { + .Call(`_AlphaMME_solveMKM`, y, X, Zlist, Klist, maxIter, tol) } #' @title Solve Multikernel RR-BLUP @@ -130,6 +131,50 @@ solveRRBLUPMK <- function(y, X, Mlist, maxIter = 40L) { .Call(`_AlphaMME_solveRRBLUPMK`, y, X, Mlist, maxIter) } +#' @title Solve RR-BLUP with EM +#' +#' @description +#' Solves a univariate mixed model of form \eqn{y=X\beta+Mu+e} using +#' the Expectation-Maximization algorithm. +#' +#' @param Y a matrix with n rows and 1 column +#' @param X a matrix with n rows and x columns +#' @param M a matrix with n rows and m columns +#' @param Vu initial guess for variance of marker effects +#' @param Ve initial guess for error variance +#' @param tol tolerance for declaring convergence +#' @param maxIter maximum iteration for attempting convergence +#' @param useEM should EM algorithm be used. If false, no estimation of +#' variance components is performed. The initial values are treated as true. +#' +#' @export +solveRRBLUP_EM <- function(Y, X, M, Vu, Ve, tol = 1e-6, maxIter = 100L, useEM = TRUE) { + .Call(`_AlphaMME_solveRRBLUP_EM`, Y, X, M, Vu, Ve, tol, maxIter, useEM) +} + +#' @title Solve RR-BLUP with EM and 2 random effects +#' +#' @description +#' Solves a univariate mixed model of form \eqn{y=X\beta+M_1u_1+M_2u_2+e} using +#' the Expectation-Maximization algorithm. +#' +#' @param Y a matrix with n rows and 1 column +#' @param X a matrix with n rows and x columns +#' @param M1 a matrix with n rows and m1 columns +#' @param M2 a matrix with n rows and m2 columns +#' @param Vu1 initial guess for variance of the first marker effects +#' @param Vu2 initial guess for variance of the second marker effects +#' @param Ve initial guess for error variance +#' @param tol tolerance for declaring convergence +#' @param maxIter maximum iteration for attempting convergence +#' @param useEM should EM algorithm be used. If false, no estimation of +#' variance components is performed. The initial values are treated as true. +#' +#' @export +solveRRBLUP_EM2 <- function(Y, X, M1, M2, Vu1, Vu2, Ve, tol = 1e-6, maxIter = 100L, useEM = TRUE) { + .Call(`_AlphaMME_solveRRBLUP_EM2`, Y, X, M1, M2, Vu1, Vu2, Ve, tol, maxIter, useEM) +} + #' @title Calculate G Matrix #' #' @description diff --git a/man/AlphaMME-package.Rd b/man/AlphaMME-package.Rd deleted file mode 100644 index 489a7b3..0000000 --- a/man/AlphaMME-package.Rd +++ /dev/null @@ -1,9 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AlphaMME.R -\docType{package} -\name{AlphaMME-package} -\alias{AlphaMME-package} -\title{Alpha Mixed Model Equations} -\description{ -This package contains mixed model equation solvers written in C++. -} diff --git a/man/AlphaMME.Rd b/man/AlphaMME.Rd new file mode 100644 index 0000000..64b41f9 --- /dev/null +++ b/man/AlphaMME.Rd @@ -0,0 +1,13 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AlphaMME.R +\docType{package} +\name{AlphaMME} +\alias{AlphaMME} +\alias{AlphaMME-package} +\title{Alpha Mixed Model Equations} +\description{ +This package is designed for genomic selection. It contains +several mixed model equation solvers from other R packages. +The solvers have been rewritten in C++ and streamlined for +increased speed and more efficient memory usage. +} diff --git a/man/readMat.Rd b/man/readMat.Rd index 021c6f2..c65c60a 100644 --- a/man/readMat.Rd +++ b/man/readMat.Rd @@ -4,16 +4,17 @@ \alias{readMat} \title{Read Matrix} \usage{ -readMat(fileName, rows, cols, sep = " ", skipRows = 0L, skipCols = 0L) +readMat(fileName, rows, cols, sep = " ", skipRows = 0L, + skipCols = 0L) } \arguments{ -\item{fileName}{path to the file to read} +\item{fileName}{path to the file being read} \item{rows}{number of rows to read in} \item{cols}{number of columns to read in} -\item{sep}{a single character seperating data entries} +\item{sep}{a single character delimiter seperating data entries} \item{skipRows}{number of rows to skip} diff --git a/man/solveMKM.Rd b/man/solveMKM.Rd index c590f68..91d7c93 100644 --- a/man/solveMKM.Rd +++ b/man/solveMKM.Rd @@ -4,7 +4,7 @@ \alias{solveMKM} \title{Solve Multikernel Model} \usage{ -solveMKM(y, X, Zlist, Klist, maxIter = 40L) +solveMKM(y, X, Zlist, Klist, maxIter = 40L, tol = 1e-04) } \arguments{ \item{y}{a matrix with n rows and 1 column} @@ -16,6 +16,8 @@ solveMKM(y, X, Zlist, Klist, maxIter = 40L) \item{Klist}{a list of K matrices} \item{maxIter}{maximum number of iteration} + +\item{tol}{tolerance for convergence} } \description{ Solves a univariate mixed model with multiple random effects. diff --git a/man/solveRRBLUP_EM.Rd b/man/solveRRBLUP_EM.Rd new file mode 100644 index 0000000..83ea7aa --- /dev/null +++ b/man/solveRRBLUP_EM.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{solveRRBLUP_EM} +\alias{solveRRBLUP_EM} +\title{Solve RR-BLUP with EM} +\usage{ +solveRRBLUP_EM(Y, X, M, Vu, Ve, tol = 1e-06, maxIter = 100L, + useEM = TRUE) +} +\arguments{ +\item{Y}{a matrix with n rows and 1 column} + +\item{X}{a matrix with n rows and x columns} + +\item{M}{a matrix with n rows and m columns} + +\item{Vu}{initial guess for variance of marker effects} + +\item{Ve}{initial guess for error variance} + +\item{tol}{tolerance for declaring convergence} + +\item{maxIter}{maximum iteration for attempting convergence} + +\item{useEM}{should EM algorithm be used. If false, no estimation of +variance components is performed. The initial values are treated as true.} +} +\description{ +Solves a univariate mixed model of form \eqn{y=X\beta+Mu+e} using +the Expectation-Maximization algorithm. +} diff --git a/man/solveRRBLUP_EM2.Rd b/man/solveRRBLUP_EM2.Rd new file mode 100644 index 0000000..940a5d3 --- /dev/null +++ b/man/solveRRBLUP_EM2.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{solveRRBLUP_EM2} +\alias{solveRRBLUP_EM2} +\title{Solve RR-BLUP with EM and 2 random effects} +\usage{ +solveRRBLUP_EM2(Y, X, M1, M2, Vu1, Vu2, Ve, tol = 1e-06, + maxIter = 100L, useEM = TRUE) +} +\arguments{ +\item{Y}{a matrix with n rows and 1 column} + +\item{X}{a matrix with n rows and x columns} + +\item{M1}{a matrix with n rows and m1 columns} + +\item{M2}{a matrix with n rows and m2 columns} + +\item{Vu1}{initial guess for variance of the first marker effects} + +\item{Vu2}{initial guess for variance of the second marker effects} + +\item{Ve}{initial guess for error variance} + +\item{tol}{tolerance for declaring convergence} + +\item{maxIter}{maximum iteration for attempting convergence} + +\item{useEM}{should EM algorithm be used. If false, no estimation of +variance components is performed. The initial values are treated as true.} +} +\description{ +Solves a univariate mixed model of form \eqn{y=X\beta+M_1u_1+M_2u_2+e} using +the Expectation-Maximization algorithm. +} diff --git a/src/MME.cpp b/src/MME.cpp index 0ad6693..c18532a 100644 --- a/src/MME.cpp +++ b/src/MME.cpp @@ -3,13 +3,27 @@ #include "alphamme.h" #include #include +#ifdef ARMA_USE_LAPACK -// Note: Fortran compiler appends '_' to subroutine name -// See http://www.netlib.org/lapack/explore-html/ for description of args -extern "C" void dsyevr_(char* JOBZ, char* RANGE, char* UPLO, long long int* N, double* A, long long int* LDA, double* VL, - double* VU, long long int* IL, long long int* IU, double* ABSTOL, long long int* M, double* W, double* Z, - long long int* LDZ, long long int* ISUPPZ, double* WORK, long long int* LWORK, long long int* IWORK, - long long int* LIWORK, long long int* INFO); +#if !defined(ARMA_BLAS_CAPITALS) +#define arma_dsyevr dsyevr +#else +#define arma_dsyevr DSYEVR +#endif + +extern "C" +void arma_fortran(arma_dsyevr)(char* JOBZ, char* RANGE, char* UPLO, long long int* N, double* A, long long int* LDA, double* VL, + double* VU, long long int* IL, long long int* IU, double* ABSTOL, long long int* M, double* W, double* Z, + long long int* LDZ, long long int* ISUPPZ, double* WORK, long long int* LWORK, long long int* IWORK, + long long int* LIWORK, long long int* INFO); +#endif + +// // Note: Fortran compiler appends '_' to subroutine name +// // See http://www.netlib.org/lapack/explore-html/ for description of args +// extern "C" void dsyevr_(char* JOBZ, char* RANGE, char* UPLO, long long int* N, double* A, long long int* LDA, double* VL, +// double* VU, long long int* IL, long long int* IU, double* ABSTOL, long long int* M, double* W, double* Z, +// long long int* LDZ, long long int* ISUPPZ, double* WORK, long long int* LWORK, long long int* IWORK, +// long long int* LIWORK, long long int* INFO); // Replacement for Armadillo's eig_sym // Fixes an error with decompisition of large matrices on Eddie @@ -46,15 +60,15 @@ int eigen2(arma::vec& eigval, arma::mat& eigvec, arma::mat X, long long int LIWORK = -1; // To be calculated long long int INFO = 0; // Calculate LWORK and LIWORK - dsyevr_(&JOBZ,&RANGE,&UPLO,&N,&*X.begin(),&LDA,&VL,&VU,&IL,&IU,&ABSTOL,&M,&*eigval.begin(), - &*eigvec.begin(),&LDZ,&*ISUPPZ.begin(),&tmpWORK,&LWORK,&tmpIWORK,&LIWORK,&INFO); + F77_CALL(dsyevr)(&JOBZ,&RANGE,&UPLO,&N,&*X.begin(),&LDA,&VL,&VU,&IL,&IU,&ABSTOL,&M,&*eigval.begin(), + &*eigvec.begin(),&LDZ,&*ISUPPZ.begin(),&tmpWORK,&LWORK,&tmpIWORK,&LIWORK,&INFO); LWORK = (long long int) tmpWORK; LIWORK = tmpIWORK; // Allocate WORK and IWORK arma::vec WORK(LWORK); arma::Col IWORK(LIWORK); // Perform decomposition - dsyevr_(&JOBZ,&RANGE,&UPLO,&N,&*X.begin(),&LDA,&VL,&VU,&IL,&IU,&ABSTOL,&M,&*eigval.begin(), + F77_CALL(dsyevr)(&JOBZ,&RANGE,&UPLO,&N,&*X.begin(),&LDA,&VL,&VU,&IL,&IU,&ABSTOL,&M,&*eigval.begin(), &*eigvec.begin(),&LDZ,&*ISUPPZ.begin(),&*WORK.begin(),&LWORK,&*IWORK.begin(),&LIWORK,&INFO); return INFO; // Return error code } @@ -77,10 +91,10 @@ Rcpp::List objREML(double param, Rcpp::List args){ //' file. Requires knowledge of the number of rows //' and columns in the file. //' -//' @param fileName path to the file to read +//' @param fileName path to the file being read //' @param rows number of rows to read in //' @param cols number of columns to read in -//' @param sep a single character seperating data entries +//' @param sep a single character delimiter seperating data entries //' @param skipRows number of rows to skip //' @param skipCols number of columns to skip //' @@ -476,14 +490,14 @@ Rcpp::List solveRRBLUPMV(const arma::mat& Y, const arma::mat& X, //' @param Zlist a list of Z matrices //' @param Klist a list of K matrices //' @param maxIter maximum number of iteration +//' @param tol tolerance for convergence //' //' @export // [[Rcpp::export]] Rcpp::List solveMKM(arma::mat& y, arma::mat& X, arma::field& Zlist, arma::field& Klist, - int maxIter=40){ - double tol = 1e-4; + int maxIter=40, double tol=1e-4){ int k = Klist.n_elem; int n = y.n_rows; int q = X.n_cols; @@ -696,6 +710,179 @@ Rcpp::List solveRRBLUPMK(arma::mat& y, arma::mat& X, Rcpp::Named("iter")=iter); } +//' @title Solve RR-BLUP with EM +//' +//' @description +//' Solves a univariate mixed model of form \eqn{y=X\beta+Mu+e} using +//' the Expectation-Maximization algorithm. +//' +//' @param Y a matrix with n rows and 1 column +//' @param X a matrix with n rows and x columns +//' @param M a matrix with n rows and m columns +//' @param Vu initial guess for variance of marker effects +//' @param Ve initial guess for error variance +//' @param tol tolerance for declaring convergence +//' @param maxIter maximum iteration for attempting convergence +//' @param useEM should EM algorithm be used. If false, no estimation of +//' variance components is performed. The initial values are treated as true. +//' +//' @export +// [[Rcpp::export]] +Rcpp::List solveRRBLUP_EM(const arma::mat& Y, const arma::mat& X, + const arma::mat& M, double Vu, double Ve, + double tol=1e-6, int maxIter=100, + bool useEM=true){ + double lambda = Ve/Vu; + double delta=0; + int iter=0; + arma::uword n=Y.n_rows,m=M.n_cols,q=X.n_cols; + arma::mat RHS(q+m,q+m),LHS(q+m,1),Rvec(q+m,1); + RHS(arma::span(0,q-1),arma::span(0,q-1)) = X.t()*X; + RHS(arma::span(0,q-1),arma::span(q,q+m-1)) = X.t()*M; + RHS(arma::span(q,q+m-1),arma::span(0,q-1)) = M.t()*X; + RHS(arma::span(q,q+m-1),arma::span(q,q+m-1)) = M.t()*M; + RHS(arma::span(q,q+m-1),arma::span(q,q+m-1)).diag() += lambda; + Rvec(arma::span(0,q-1),0) = X.t()*Y; + Rvec(arma::span(q,q+m-1),0) = M.t()*Y; + arma::mat RHSinv = pinv(RHS); + LHS = RHSinv*Rvec; + if(useEM){ + Ve = as_scalar(Y.t()*Y-LHS.t()*Rvec)/(n-q); + Vu = as_scalar( + LHS(arma::span(q,q+m-1),0).t()*LHS(arma::span(q,q+m-1),0)+ + Ve*sum(RHSinv(arma::span(q,q+m-1),arma::span(q,q+m-1)).diag()) + )/m; + delta = Ve/Vu-lambda; + while(fabs(delta)>tol){ + RHS(arma::span(q,q+m-1),arma::span(q,q+m-1)).diag() += delta; + lambda += delta; + RHSinv = pinv(RHS); + LHS = RHSinv*Rvec; + iter++; + if(iter>=maxIter){ + Rcpp::Rcerr<<"Warning: did not converge, reached maxIter\n"; + break; + } + Ve = as_scalar(Y.t()*Y-LHS.t()*Rvec)/(n-q); + Vu = as_scalar( + LHS(arma::span(q,q+m-1),0).t()*LHS(arma::span(q,q+m-1),0)+ + Ve*sum(RHSinv(arma::span(q,q+m-1),arma::span(q,q+m-1)).diag()) + )/m; + delta = Ve/Vu-lambda; + } + } + return Rcpp::List::create(Rcpp::Named("Vu")=Vu, + Rcpp::Named("Ve")=Ve, + Rcpp::Named("beta")=LHS.rows(arma::span(0,q-1)), + Rcpp::Named("u")=LHS.rows(arma::span(q,q+m-1)), + Rcpp::Named("iter")=iter); +} + +//' @title Solve RR-BLUP with EM and 2 random effects +//' +//' @description +//' Solves a univariate mixed model of form \eqn{y=X\beta+M_1u_1+M_2u_2+e} using +//' the Expectation-Maximization algorithm. +//' +//' @param Y a matrix with n rows and 1 column +//' @param X a matrix with n rows and x columns +//' @param M1 a matrix with n rows and m1 columns +//' @param M2 a matrix with n rows and m2 columns +//' @param Vu1 initial guess for variance of the first marker effects +//' @param Vu2 initial guess for variance of the second marker effects +//' @param Ve initial guess for error variance +//' @param tol tolerance for declaring convergence +//' @param maxIter maximum iteration for attempting convergence +//' @param useEM should EM algorithm be used. If false, no estimation of +//' variance components is performed. The initial values are treated as true. +//' +//' @export +// [[Rcpp::export]] +Rcpp::List solveRRBLUP_EM2(const arma::mat& Y, const arma::mat& X, + const arma::mat& M1, const arma::mat& M2, + double Vu1, double Vu2, double Ve, + double tol=1e-6, int maxIter=100, + bool useEM=true){ + double lambda1 = Ve/Vu1; + double lambda2 = Ve/Vu2; + double delta1=0,delta2=0,VeN=0,Vu1N=0,Vu2N=0; + int iter=0; + arma::uword n=Y.n_rows,m=M1.n_cols,q=X.n_cols; + arma::mat RHS(q+2*m,q+2*m),LHS(q+2*m,1),Rvec(q+2*m,1); + // Top row + RHS(arma::span(0,q-1),arma::span(0,q-1)) = X.t()*X; + RHS(arma::span(0,q-1),arma::span(q,q+m-1)) = X.t()*M1; + RHS(arma::span(0,q-1),arma::span(q+m,q+2*m-1)) = X.t()*M2; + // Second row + RHS(arma::span(q,q+m-1),arma::span(0,q-1)) = M1.t()*X; + RHS(arma::span(q,q+m-1),arma::span(q,q+m-1)) = M1.t()*M1; + RHS(arma::span(q,q+m-1),arma::span(q+m,q+2*m-1)) = M1.t()*M2; + // Third row + RHS(arma::span(q+m,q+2*m-1),arma::span(0,q-1)) = M2.t()*X; + RHS(arma::span(q+m,q+2*m-1),arma::span(q,q+m-1)) = M2.t()*M1; + RHS(arma::span(q+m,q+2*m-1),arma::span(q+m,q+2*m-1)) = M2.t()*M2; + // Add to diagonal + RHS(arma::span(q,q+m-1),arma::span(q,q+m-1)).diag() += lambda1; + RHS(arma::span(q+m,q+2*m-1),arma::span(q+m,q+2*m-1)).diag() += lambda2; + Rvec(arma::span(0,q-1),0) = X.t()*Y; + Rvec(arma::span(q,q+m-1),0) = M1.t()*Y; + Rvec(arma::span(q+m,q+2*m-1),0) = M2.t()*Y; + arma::mat RHSinv = inv(RHS); + LHS = RHSinv*Rvec; + if(useEM){ + VeN = as_scalar(Y.t()*Y-LHS.t()*Rvec)/(n-q); + Vu1N = as_scalar( + LHS(arma::span(q,q+m-1),0).t()*LHS(arma::span(q,q+m-1),0)+ + Ve*sum(RHSinv(arma::span(q,q+m-1),arma::span(q,q+m-1)).diag()) + )/m; + Vu2N = as_scalar( + LHS(arma::span(q+m,q+2*m-1),0).t()*LHS(arma::span(q+m,q+2*m-1),0)+ + Ve*sum(RHSinv(arma::span(q+m,q+2*m-1),arma::span(q+m,q+2*m-1)).diag()) + )/m; + delta1 = VeN/Vu1N-lambda1; + delta2 = VeN/Vu2N-lambda2; + while((fabs(delta1)>tol) || (fabs(delta2)>tol)){ + Ve = VeN; + Vu1 = Vu1N; + Vu2 = Vu2N; + RHS(arma::span(q,q+m-1),arma::span(q,q+m-1)).diag() += delta1; + RHS(arma::span(q+m,q+2*m-1),arma::span(q+m,q+2*m-1)).diag() += delta2; + lambda1 += delta1; + lambda2 += delta2; + RHSinv = inv(RHS); + LHS = RHSinv*Rvec; + iter++; + if(iter>=maxIter){ + Rcpp::Rcerr<<"Warning: did not converge, reached maxIter\n"; + break; + } + VeN = as_scalar(Y.t()*Y-LHS.t()*Rvec)/(n-q); + Vu1N = as_scalar( + LHS(arma::span(q,q+m-1),0).t()*LHS(arma::span(q,q+m-1),0)+ + Ve*sum(RHSinv(arma::span(q,q+m-1),arma::span(q,q+m-1)).diag()) + )/m; + Vu2N = as_scalar( + LHS(arma::span(q+m,q+2*m-1),0).t()*LHS(arma::span(q+m,q+2*m-1),0)+ + Ve*sum(RHSinv(arma::span(q+m,q+2*m-1),arma::span(q+m,q+2*m-1)).diag()) + )/m; + delta1 = VeN/Vu1N-lambda1; + delta2 = VeN/Vu2N-lambda2; + } + } + arma::vec Vu(2); + Vu(0) = Vu1; + Vu(1) = Vu2; + arma::mat u(m,2); + u.col(0) = LHS.rows(arma::span(q,q+m-1)); + u.col(1) = LHS.rows(arma::span(q+m,q+2*m-1)); + return Rcpp::List::create(Rcpp::Named("Vu")=Vu, + Rcpp::Named("Ve")=Ve, + Rcpp::Named("beta")=LHS.rows(arma::span(0,q-1)), + Rcpp::Named("u")=u, + Rcpp::Named("iter")=iter); +} + + //' @title Calculate G Matrix //' //' @description diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 9ace088..43e7dde 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -94,8 +94,8 @@ BEGIN_RCPP END_RCPP } // solveMKM -Rcpp::List solveMKM(arma::mat& y, arma::mat& X, arma::field& Zlist, arma::field& Klist, int maxIter); -RcppExport SEXP _AlphaMME_solveMKM(SEXP ySEXP, SEXP XSEXP, SEXP ZlistSEXP, SEXP KlistSEXP, SEXP maxIterSEXP) { +Rcpp::List solveMKM(arma::mat& y, arma::mat& X, arma::field& Zlist, arma::field& Klist, int maxIter, double tol); +RcppExport SEXP _AlphaMME_solveMKM(SEXP ySEXP, SEXP XSEXP, SEXP ZlistSEXP, SEXP KlistSEXP, SEXP maxIterSEXP, SEXP tolSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -104,7 +104,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::field& >::type Zlist(ZlistSEXP); Rcpp::traits::input_parameter< arma::field& >::type Klist(KlistSEXP); Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP); - rcpp_result_gen = Rcpp::wrap(solveMKM(y, X, Zlist, Klist, maxIter)); + Rcpp::traits::input_parameter< double >::type tol(tolSEXP); + rcpp_result_gen = Rcpp::wrap(solveMKM(y, X, Zlist, Klist, maxIter, tol)); return rcpp_result_gen; END_RCPP } @@ -122,6 +123,44 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// solveRRBLUP_EM +Rcpp::List solveRRBLUP_EM(const arma::mat& Y, const arma::mat& X, const arma::mat& M, double Vu, double Ve, double tol, int maxIter, bool useEM); +RcppExport SEXP _AlphaMME_solveRRBLUP_EM(SEXP YSEXP, SEXP XSEXP, SEXP MSEXP, SEXP VuSEXP, SEXP VeSEXP, SEXP tolSEXP, SEXP maxIterSEXP, SEXP useEMSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type M(MSEXP); + Rcpp::traits::input_parameter< double >::type Vu(VuSEXP); + Rcpp::traits::input_parameter< double >::type Ve(VeSEXP); + Rcpp::traits::input_parameter< double >::type tol(tolSEXP); + Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP); + Rcpp::traits::input_parameter< bool >::type useEM(useEMSEXP); + rcpp_result_gen = Rcpp::wrap(solveRRBLUP_EM(Y, X, M, Vu, Ve, tol, maxIter, useEM)); + return rcpp_result_gen; +END_RCPP +} +// solveRRBLUP_EM2 +Rcpp::List solveRRBLUP_EM2(const arma::mat& Y, const arma::mat& X, const arma::mat& M1, const arma::mat& M2, double Vu1, double Vu2, double Ve, double tol, int maxIter, bool useEM); +RcppExport SEXP _AlphaMME_solveRRBLUP_EM2(SEXP YSEXP, SEXP XSEXP, SEXP M1SEXP, SEXP M2SEXP, SEXP Vu1SEXP, SEXP Vu2SEXP, SEXP VeSEXP, SEXP tolSEXP, SEXP maxIterSEXP, SEXP useEMSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type M1(M1SEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type M2(M2SEXP); + Rcpp::traits::input_parameter< double >::type Vu1(Vu1SEXP); + Rcpp::traits::input_parameter< double >::type Vu2(Vu2SEXP); + Rcpp::traits::input_parameter< double >::type Ve(VeSEXP); + Rcpp::traits::input_parameter< double >::type tol(tolSEXP); + Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP); + Rcpp::traits::input_parameter< bool >::type useEM(useEMSEXP); + rcpp_result_gen = Rcpp::wrap(solveRRBLUP_EM2(Y, X, M1, M2, Vu1, Vu2, Ve, tol, maxIter, useEM)); + return rcpp_result_gen; +END_RCPP +} // calcG arma::mat calcG(arma::mat X); RcppExport SEXP _AlphaMME_calcG(SEXP XSEXP) { diff --git a/src/init.c b/src/init.c index 901e43e..6d58200 100644 --- a/src/init.c +++ b/src/init.c @@ -16,28 +16,32 @@ extern SEXP _AlphaMME_fastPairDist(SEXP, SEXP); extern SEXP _AlphaMME_gaussKernel(SEXP, SEXP); extern SEXP _AlphaMME_readMat(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _AlphaMME_solveAniModel(SEXP, SEXP, SEXP); -extern SEXP _AlphaMME_solveMKM(SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _AlphaMME_solveMKM(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _AlphaMME_solveMVM(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _AlphaMME_solveRRBLUP(SEXP, SEXP, SEXP); +extern SEXP _AlphaMME_solveRRBLUP_EM(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _AlphaMME_solveRRBLUP_EM2(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _AlphaMME_solveRRBLUPMK(SEXP, SEXP, SEXP, SEXP); extern SEXP _AlphaMME_solveRRBLUPMV(SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _AlphaMME_solveUVM(SEXP, SEXP, SEXP, SEXP); static const R_CallMethodDef CallEntries[] = { - {"_AlphaMME_calcD", (DL_FUNC) &_AlphaMME_calcD, 1}, - {"_AlphaMME_calcG", (DL_FUNC) &_AlphaMME_calcG, 1}, - {"_AlphaMME_calcGIbs", (DL_FUNC) &_AlphaMME_calcGIbs, 1}, - {"_AlphaMME_fastDist", (DL_FUNC) &_AlphaMME_fastDist, 1}, - {"_AlphaMME_fastPairDist", (DL_FUNC) &_AlphaMME_fastPairDist, 2}, - {"_AlphaMME_gaussKernel", (DL_FUNC) &_AlphaMME_gaussKernel, 2}, - {"_AlphaMME_readMat", (DL_FUNC) &_AlphaMME_readMat, 6}, - {"_AlphaMME_solveAniModel", (DL_FUNC) &_AlphaMME_solveAniModel, 3}, - {"_AlphaMME_solveMKM", (DL_FUNC) &_AlphaMME_solveMKM, 5}, - {"_AlphaMME_solveMVM", (DL_FUNC) &_AlphaMME_solveMVM, 6}, - {"_AlphaMME_solveRRBLUP", (DL_FUNC) &_AlphaMME_solveRRBLUP, 3}, - {"_AlphaMME_solveRRBLUPMK", (DL_FUNC) &_AlphaMME_solveRRBLUPMK, 4}, - {"_AlphaMME_solveRRBLUPMV", (DL_FUNC) &_AlphaMME_solveRRBLUPMV, 5}, - {"_AlphaMME_solveUVM", (DL_FUNC) &_AlphaMME_solveUVM, 4}, + {"_AlphaMME_calcD", (DL_FUNC) &_AlphaMME_calcD, 1}, + {"_AlphaMME_calcG", (DL_FUNC) &_AlphaMME_calcG, 1}, + {"_AlphaMME_calcGIbs", (DL_FUNC) &_AlphaMME_calcGIbs, 1}, + {"_AlphaMME_fastDist", (DL_FUNC) &_AlphaMME_fastDist, 1}, + {"_AlphaMME_fastPairDist", (DL_FUNC) &_AlphaMME_fastPairDist, 2}, + {"_AlphaMME_gaussKernel", (DL_FUNC) &_AlphaMME_gaussKernel, 2}, + {"_AlphaMME_readMat", (DL_FUNC) &_AlphaMME_readMat, 6}, + {"_AlphaMME_solveAniModel", (DL_FUNC) &_AlphaMME_solveAniModel, 3}, + {"_AlphaMME_solveMKM", (DL_FUNC) &_AlphaMME_solveMKM, 6}, + {"_AlphaMME_solveMVM", (DL_FUNC) &_AlphaMME_solveMVM, 6}, + {"_AlphaMME_solveRRBLUP", (DL_FUNC) &_AlphaMME_solveRRBLUP, 3}, + {"_AlphaMME_solveRRBLUP_EM", (DL_FUNC) &_AlphaMME_solveRRBLUP_EM, 8}, + {"_AlphaMME_solveRRBLUP_EM2", (DL_FUNC) &_AlphaMME_solveRRBLUP_EM2, 10}, + {"_AlphaMME_solveRRBLUPMK", (DL_FUNC) &_AlphaMME_solveRRBLUPMK, 4}, + {"_AlphaMME_solveRRBLUPMV", (DL_FUNC) &_AlphaMME_solveRRBLUPMV, 5}, + {"_AlphaMME_solveUVM", (DL_FUNC) &_AlphaMME_solveUVM, 4}, {NULL, NULL, 0} }; diff --git a/src/optimize.cpp b/src/optimize.cpp index aafc0fd..449e43f 100644 --- a/src/optimize.cpp +++ b/src/optimize.cpp @@ -15,8 +15,8 @@ using namespace Rcpp; * maximize searches for a maximum value * eps controls the tolerance for convergence */ -List optimize(List (*objective)(double, List), List args, double l, - double u, int maxIter=1000, bool maximize=false, +List optimize(List (*objective)(double, List), List args, double l, + double u, int maxIter=1000, bool maximize=false, bool evalU=false, bool evalL=false, double eps=1.0e-9){ double MACHEPS_SQRT = sqrt(2.2204460492503131e-016); double c = (3-sqrt(5))/2; @@ -32,26 +32,26 @@ List optimize(List (*objective)(double, List), List args, double l, if(maximize) fx = -fx; double fv = fx; double fw = fx; - + int numiter = 0; - + while(true){ double m = 0.5*(l+u); double tol = MACHEPS_SQRT*std::abs(x)+eps; double tol2 = 2.0*tol; - + // Check the stopping criterion if(std::abs(x-m)<=(tol2-(0.5*(u-l)))){ break; } - + // Check maximum iterations if (++numiter>maxIter){ break; } - + double p=0,q=0,r=0,d=0,z=0; - + if(std::abs(e)>tol){ // Fit parabola r = (x-w)*(fx-fv); @@ -66,7 +66,7 @@ List optimize(List (*objective)(double, List), List args, double l, r = e; e = d; } - + if((std::abs(p)=tol){ z = x+d; @@ -94,7 +94,7 @@ List optimize(List (*objective)(double, List), List args, double l, fOut = objective(z,args); double funcu = fOut["objective"]; if(maximize) funcu = -funcu; - + // Update if(funcu<=fx){ if(z