From 49194747f1af29ae6fd51b7d5d86525a7be9ccac Mon Sep 17 00:00:00 2001 From: GAYNOR Chris Date: Thu, 1 Feb 2018 13:07:30 +0000 Subject: [PATCH] Added code --- .Rbuildignore | 3 + .gitignore | 9 + AlphaMME.Rproj | 21 + DESCRIPTION | 13 + NAMESPACE | 19 + R/AlphaMME.R | 11 + R/RcppExports.R | 228 +++++++++++ man/AlphaMME-package.Rd | 9 + man/calcD.Rd | 20 + man/calcG.Rd | 17 + man/calcGIbs.Rd | 18 + man/fastDist.Rd | 19 + man/fastPairDist.Rd | 21 + man/gaussKernel.Rd | 21 + man/readMat.Rd | 29 ++ man/solveAniModel.Rd | 19 + man/solveMKM.Rd | 22 + man/solveMVM.Rd | 24 ++ man/solveRRBLUP.Rd | 18 + man/solveRRBLUPMK.Rd | 20 + man/solveRRBLUPMV.Rd | 22 + man/solveUVM.Rd | 20 + src/MME.cpp | 884 ++++++++++++++++++++++++++++++++++++++++ src/Makevars | 3 + src/Makevars.win | 3 + src/RcppExports.cpp | 192 +++++++++ src/alphamme.h | 7 + src/init.c | 48 +++ src/optimize.cpp | 155 +++++++ src/optimize.h | 9 + 30 files changed, 1904 insertions(+) create mode 100644 .Rbuildignore create mode 100644 .gitignore create mode 100644 AlphaMME.Rproj create mode 100644 DESCRIPTION create mode 100644 NAMESPACE create mode 100644 R/AlphaMME.R create mode 100644 R/RcppExports.R create mode 100644 man/AlphaMME-package.Rd create mode 100644 man/calcD.Rd create mode 100644 man/calcG.Rd create mode 100644 man/calcGIbs.Rd create mode 100644 man/fastDist.Rd create mode 100644 man/fastPairDist.Rd create mode 100644 man/gaussKernel.Rd create mode 100644 man/readMat.Rd create mode 100644 man/solveAniModel.Rd create mode 100644 man/solveMKM.Rd create mode 100644 man/solveMVM.Rd create mode 100644 man/solveRRBLUP.Rd create mode 100644 man/solveRRBLUPMK.Rd create mode 100644 man/solveRRBLUPMV.Rd create mode 100644 man/solveUVM.Rd create mode 100644 src/MME.cpp create mode 100644 src/Makevars create mode 100644 src/Makevars.win create mode 100644 src/RcppExports.cpp create mode 100644 src/alphamme.h create mode 100644 src/init.c create mode 100644 src/optimize.cpp create mode 100644 src/optimize.h diff --git a/.Rbuildignore b/.Rbuildignore new file mode 100644 index 0000000..17435f3 --- /dev/null +++ b/.Rbuildignore @@ -0,0 +1,3 @@ +^.*\.Rproj$ +^\.Rproj\.user$ +Notes.txt diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b2e2242 --- /dev/null +++ b/.gitignore @@ -0,0 +1,9 @@ +.Rproj.user +.Rhistory +.RData +.Ruserdata +src/*.o +src/*.so +src/*.dll +Notes.txt + diff --git a/AlphaMME.Rproj b/AlphaMME.Rproj new file mode 100644 index 0000000..270314b --- /dev/null +++ b/AlphaMME.Rproj @@ -0,0 +1,21 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX + +AutoAppendNewline: Yes +StripTrailingWhitespace: Yes + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source +PackageRoxygenize: rd,collate,namespace diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..a299e80 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,13 @@ +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++. +License: MIT + file LICENSE +URL: https://bitbucket.org/hickeyjohnteam/alphamme +Imports: Rcpp (>= 0.12.13), RcppArmadillo +LinkingTo: Rcpp, RcppArmadillo +RoxygenNote: 6.0.1 +NeedsCompilation: true diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..1b8094d --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,19 @@ +# Generated by roxygen2: do not edit by hand + +export(calcD) +export(calcG) +export(calcGIbs) +export(fastDist) +export(fastPairDist) +export(gaussKernel) +export(readMat) +export(solveAniModel) +export(solveMKM) +export(solveMVM) +export(solveRRBLUP) +export(solveRRBLUPMK) +export(solveRRBLUPMV) +export(solveUVM) +import(Rcpp) +import(RcppArmadillo) +useDynLib(AlphaMME, .registration = TRUE) diff --git a/R/AlphaMME.R b/R/AlphaMME.R new file mode 100644 index 0000000..90f5daa --- /dev/null +++ b/R/AlphaMME.R @@ -0,0 +1,11 @@ +#' @useDynLib AlphaMME, .registration = TRUE +#' @import Rcpp RcppArmadillo + +#' @title Alpha Mixed Model Equations +#' +#' @description +#' This package contains mixed model equation solvers written in C++. +#' +#' @docType package +#' @name AlphaMME-package +NULL \ No newline at end of file diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..1c7c1c0 --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,228 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#' @title Read Matrix +#' +#' @description +#' Uses C++ to quickly read a matrix from a text +#' file. Requires knowledge of the number of rows +#' and columns in the file. +#' +#' @param fileName path to the file to 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 skipRows number of rows to skip +#' @param skipCols number of columns to skip +#' +#' @return a numeric matrix +#' +#' @export +readMat <- function(fileName, rows, cols, sep = ' ', skipRows = 0L, skipCols = 0L) { + .Call(`_AlphaMME_readMat`, fileName, rows, cols, sep, skipRows, skipCols) +} + +#' @title Solve Univariate Model +#' +#' @description +#' Solves a univariate mixed model of form \eqn{y=X\beta+Zu+e} +#' +#' @param y a matrix with n rows and 1 column +#' @param X a matrix with n rows and x columns +#' @param Z a matrix with n rows and m columns +#' @param K a matrix with m rows and m columns +#' +#' @export +solveUVM <- function(y, X, Z, K) { + .Call(`_AlphaMME_solveUVM`, y, X, Z, K) +} + +#' @title Solve animal model +#' +#' @description +#' Solves a univariate mixed model of form \eqn{y=X\beta+u+e} +#' +#' @param y a matrix with n rows and 1 column +#' @param X a matrix with n rows and x columns +#' @param K the numeric relationship matrix +#' with n rows and n columns +#' +#' @export +solveAniModel <- function(y, X, K) { + .Call(`_AlphaMME_solveAniModel`, y, X, K) +} + +#' @title Solve RR-BLUP +#' +#' @description +#' Solves a univariate mixed model of form \eqn{y=X\beta+Mu+e} +#' +#' @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 +#' +#' @export +solveRRBLUP <- function(y, X, M) { + .Call(`_AlphaMME_solveRRBLUP`, y, X, M) +} + +#' @title Solve Multivariate Model +#' +#' @description +#' Solves a multivariate mixed model of form \eqn{Y=X\beta+Zu+e} +#' +#' @param Y a matrix with n rows and q columns +#' @param X a matrix with n rows and x columns +#' @param Z a matrix with n rows and m columns +#' @param K a matrix with m rows and m columns +#' @param tol tolerance for convergence +#' @param maxIter maximum number of iteration +#' +#' @export +solveMVM <- function(Y, X, Z, K, tol = 1e-6, maxIter = 1000L) { + .Call(`_AlphaMME_solveMVM`, Y, X, Z, K, tol, maxIter) +} + +#' @title Solve Multivariate RR-BLUP +#' +#' @description +#' Solves a multivariate mixed model of form \eqn{Y=X\beta+Mu+e} +#' +#' @param Y a matrix with n rows and q columns +#' @param X a matrix with n rows and x columns +#' @param M a matrix with n rows and m columns +#' @param tol tolerance for convergence +#' @param maxIter maximum number of iteration +#' +#' @export +solveRRBLUPMV <- function(Y, X, M, tol = 1e-6, maxIter = 1000L) { + .Call(`_AlphaMME_solveRRBLUPMV`, Y, X, M, tol, maxIter) +} + +#' @title Solve Multikernel Model +#' +#' @description +#' Solves a univariate mixed model with multiple random effects. +#' +#' @param y a matrix with n rows and 1 column +#' @param X a matrix with n rows and x columns +#' @param Zlist a list of Z matrices +#' @param Klist a list of K matrices +#' @param maxIter maximum number of iteration +#' +#' @export +solveMKM <- function(y, X, Zlist, Klist, maxIter = 40L) { + .Call(`_AlphaMME_solveMKM`, y, X, Zlist, Klist, maxIter) +} + +#' @title Solve Multikernel RR-BLUP +#' +#' @description +#' Solves a univariate mixed model with multiple random effects. +#' +#' @param y a matrix with n rows and 1 column +#' @param X a matrix with n rows and x columns +#' @param Mlist a list of M matrices +#' @param maxIter maximum number of iteration +#' +#' @export +solveRRBLUPMK <- function(y, X, Mlist, maxIter = 40L) { + .Call(`_AlphaMME_solveRRBLUPMK`, y, X, Mlist, maxIter) +} + +#' @title Calculate G Matrix +#' +#' @description +#' Calculates the genomic relationship matrix. +#' +#' @param X a matrix of marker genotypes scored as 0,1,2 +#' +#' @return a matrix of the realized genomic relationships +#' +#' @export +calcG <- function(X) { + .Call(`_AlphaMME_calcG`, X) +} + +#' @title Calculate Dominance Matrix +#' +#' @description +#' Calculates the dominance relationship matrix. +#' +#' @param X a matrix of marker genotypes scored as 0,1,2 +#' +#' @references +#' \cite{Nishio, M, and M. Satoh. 2014. Including Dominance Effects in the Genomic BLUP Method for Genomic Evaluation. PLOS ONE 9(1): e85792.} +#' +#' @return a matrix of the realized dominance relationships +#' +#' @export +calcD <- function(X) { + .Call(`_AlphaMME_calcD`, X) +} + +#' @title Calculate IBS G Matrix +#' +#' @description +#' Calculates an identity-by-state genomic relationship matrix +#' based on simple matching. +#' +#' @param X a matrix of marker genotypes scored as 0,1,2 +#' +#' @return a matrix of genomic relationships +#' +#' @export +calcGIbs <- function(X) { + .Call(`_AlphaMME_calcGIbs`, X) +} + +#' @title Calculate Euclidean distance +#' +#' @description +#' Calculates a Euclidean distance matrix using a binomial +#' theorem trick. Results in much faster computation than the +#' \code{dist} function in package \code{stats}. +#' +#' @param X a numeric matrix +#' +#' @return a matrix of columnwise distances +#' +#' @export +fastDist <- function(X) { + .Call(`_AlphaMME_fastDist`, X) +} + +#' @title Calculate Paired Euclidean distance +#' +#' @description +#' Calculates a Euclidean distance between two matrices using +#' a binomial theorem trick. +#' +#' @param X a numeric matrix +#' @param Y a numeric matrix +#' +#' @return a matrix of columnwise distances between matrices +#' X and Y +#' +#' @export +fastPairDist <- function(X, Y) { + .Call(`_AlphaMME_fastPairDist`, X, Y) +} + +#' @title Calculate Gaussian Kernel +#' +#' @description +#' Calculates a Gaussian kernel using a Euclidean distance +#' matrix. +#' +#' @param D a matrix of Euclidean distances, +#' see \code{\link{fastDist}} +#' @param theta the tuning parameter +#' +#' @return a numeric matrix +#' +#' @export +gaussKernel <- function(D, theta) { + .Call(`_AlphaMME_gaussKernel`, D, theta) +} + diff --git a/man/AlphaMME-package.Rd b/man/AlphaMME-package.Rd new file mode 100644 index 0000000..489a7b3 --- /dev/null +++ b/man/AlphaMME-package.Rd @@ -0,0 +1,9 @@ +% 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/calcD.Rd b/man/calcD.Rd new file mode 100644 index 0000000..dc6c790 --- /dev/null +++ b/man/calcD.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{calcD} +\alias{calcD} +\title{Calculate Dominance Matrix} +\usage{ +calcD(X) +} +\arguments{ +\item{X}{a matrix of marker genotypes scored as 0,1,2} +} +\value{ +a matrix of the realized dominance relationships +} +\description{ +Calculates the dominance relationship matrix. +} +\references{ +\cite{Nishio, M, and M. Satoh. 2014. Including Dominance Effects in the Genomic BLUP Method for Genomic Evaluation. PLOS ONE 9(1): e85792.} +} diff --git a/man/calcG.Rd b/man/calcG.Rd new file mode 100644 index 0000000..d6c5abf --- /dev/null +++ b/man/calcG.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{calcG} +\alias{calcG} +\title{Calculate G Matrix} +\usage{ +calcG(X) +} +\arguments{ +\item{X}{a matrix of marker genotypes scored as 0,1,2} +} +\value{ +a matrix of the realized genomic relationships +} +\description{ +Calculates the genomic relationship matrix. +} diff --git a/man/calcGIbs.Rd b/man/calcGIbs.Rd new file mode 100644 index 0000000..76b2fa8 --- /dev/null +++ b/man/calcGIbs.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{calcGIbs} +\alias{calcGIbs} +\title{Calculate IBS G Matrix} +\usage{ +calcGIbs(X) +} +\arguments{ +\item{X}{a matrix of marker genotypes scored as 0,1,2} +} +\value{ +a matrix of genomic relationships +} +\description{ +Calculates an identity-by-state genomic relationship matrix +based on simple matching. +} diff --git a/man/fastDist.Rd b/man/fastDist.Rd new file mode 100644 index 0000000..a182f43 --- /dev/null +++ b/man/fastDist.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{fastDist} +\alias{fastDist} +\title{Calculate Euclidean distance} +\usage{ +fastDist(X) +} +\arguments{ +\item{X}{a numeric matrix} +} +\value{ +a matrix of columnwise distances +} +\description{ +Calculates a Euclidean distance matrix using a binomial +theorem trick. Results in much faster computation than the +\code{dist} function in package \code{stats}. +} diff --git a/man/fastPairDist.Rd b/man/fastPairDist.Rd new file mode 100644 index 0000000..0dfa277 --- /dev/null +++ b/man/fastPairDist.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{fastPairDist} +\alias{fastPairDist} +\title{Calculate Paired Euclidean distance} +\usage{ +fastPairDist(X, Y) +} +\arguments{ +\item{X}{a numeric matrix} + +\item{Y}{a numeric matrix} +} +\value{ +a matrix of columnwise distances between matrices +X and Y +} +\description{ +Calculates a Euclidean distance between two matrices using +a binomial theorem trick. +} diff --git a/man/gaussKernel.Rd b/man/gaussKernel.Rd new file mode 100644 index 0000000..e28a83f --- /dev/null +++ b/man/gaussKernel.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{gaussKernel} +\alias{gaussKernel} +\title{Calculate Gaussian Kernel} +\usage{ +gaussKernel(D, theta) +} +\arguments{ +\item{D}{a matrix of Euclidean distances, +see \code{\link{fastDist}}} + +\item{theta}{the tuning parameter} +} +\value{ +a numeric matrix +} +\description{ +Calculates a Gaussian kernel using a Euclidean distance +matrix. +} diff --git a/man/readMat.Rd b/man/readMat.Rd new file mode 100644 index 0000000..021c6f2 --- /dev/null +++ b/man/readMat.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{readMat} +\alias{readMat} +\title{Read Matrix} +\usage{ +readMat(fileName, rows, cols, sep = " ", skipRows = 0L, skipCols = 0L) +} +\arguments{ +\item{fileName}{path to the file to 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{skipRows}{number of rows to skip} + +\item{skipCols}{number of columns to skip} +} +\value{ +a numeric matrix +} +\description{ +Uses C++ to quickly read a matrix from a text +file. Requires knowledge of the number of rows +and columns in the file. +} diff --git a/man/solveAniModel.Rd b/man/solveAniModel.Rd new file mode 100644 index 0000000..493877d --- /dev/null +++ b/man/solveAniModel.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{solveAniModel} +\alias{solveAniModel} +\title{Solve animal model} +\usage{ +solveAniModel(y, X, K) +} +\arguments{ +\item{y}{a matrix with n rows and 1 column} + +\item{X}{a matrix with n rows and x columns} + +\item{K}{the numeric relationship matrix +with n rows and n columns} +} +\description{ +Solves a univariate mixed model of form \eqn{y=X\beta+u+e} +} diff --git a/man/solveMKM.Rd b/man/solveMKM.Rd new file mode 100644 index 0000000..c590f68 --- /dev/null +++ b/man/solveMKM.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{solveMKM} +\alias{solveMKM} +\title{Solve Multikernel Model} +\usage{ +solveMKM(y, X, Zlist, Klist, maxIter = 40L) +} +\arguments{ +\item{y}{a matrix with n rows and 1 column} + +\item{X}{a matrix with n rows and x columns} + +\item{Zlist}{a list of Z matrices} + +\item{Klist}{a list of K matrices} + +\item{maxIter}{maximum number of iteration} +} +\description{ +Solves a univariate mixed model with multiple random effects. +} diff --git a/man/solveMVM.Rd b/man/solveMVM.Rd new file mode 100644 index 0000000..ab9c858 --- /dev/null +++ b/man/solveMVM.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{solveMVM} +\alias{solveMVM} +\title{Solve Multivariate Model} +\usage{ +solveMVM(Y, X, Z, K, tol = 1e-06, maxIter = 1000L) +} +\arguments{ +\item{Y}{a matrix with n rows and q columns} + +\item{X}{a matrix with n rows and x columns} + +\item{Z}{a matrix with n rows and m columns} + +\item{K}{a matrix with m rows and m columns} + +\item{tol}{tolerance for convergence} + +\item{maxIter}{maximum number of iteration} +} +\description{ +Solves a multivariate mixed model of form \eqn{Y=X\beta+Zu+e} +} diff --git a/man/solveRRBLUP.Rd b/man/solveRRBLUP.Rd new file mode 100644 index 0000000..a921eec --- /dev/null +++ b/man/solveRRBLUP.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{solveRRBLUP} +\alias{solveRRBLUP} +\title{Solve RR-BLUP} +\usage{ +solveRRBLUP(y, X, M) +} +\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} +} +\description{ +Solves a univariate mixed model of form \eqn{y=X\beta+Mu+e} +} diff --git a/man/solveRRBLUPMK.Rd b/man/solveRRBLUPMK.Rd new file mode 100644 index 0000000..393f1cf --- /dev/null +++ b/man/solveRRBLUPMK.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{solveRRBLUPMK} +\alias{solveRRBLUPMK} +\title{Solve Multikernel RR-BLUP} +\usage{ +solveRRBLUPMK(y, X, Mlist, maxIter = 40L) +} +\arguments{ +\item{y}{a matrix with n rows and 1 column} + +\item{X}{a matrix with n rows and x columns} + +\item{Mlist}{a list of M matrices} + +\item{maxIter}{maximum number of iteration} +} +\description{ +Solves a univariate mixed model with multiple random effects. +} diff --git a/man/solveRRBLUPMV.Rd b/man/solveRRBLUPMV.Rd new file mode 100644 index 0000000..3573924 --- /dev/null +++ b/man/solveRRBLUPMV.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{solveRRBLUPMV} +\alias{solveRRBLUPMV} +\title{Solve Multivariate RR-BLUP} +\usage{ +solveRRBLUPMV(Y, X, M, tol = 1e-06, maxIter = 1000L) +} +\arguments{ +\item{Y}{a matrix with n rows and q columns} + +\item{X}{a matrix with n rows and x columns} + +\item{M}{a matrix with n rows and m columns} + +\item{tol}{tolerance for convergence} + +\item{maxIter}{maximum number of iteration} +} +\description{ +Solves a multivariate mixed model of form \eqn{Y=X\beta+Mu+e} +} diff --git a/man/solveUVM.Rd b/man/solveUVM.Rd new file mode 100644 index 0000000..358a0b4 --- /dev/null +++ b/man/solveUVM.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{solveUVM} +\alias{solveUVM} +\title{Solve Univariate Model} +\usage{ +solveUVM(y, X, Z, K) +} +\arguments{ +\item{y}{a matrix with n rows and 1 column} + +\item{X}{a matrix with n rows and x columns} + +\item{Z}{a matrix with n rows and m columns} + +\item{K}{a matrix with m rows and m columns} +} +\description{ +Solves a univariate mixed model of form \eqn{y=X\beta+Zu+e} +} diff --git a/src/MME.cpp b/src/MME.cpp new file mode 100644 index 0000000..7113dcc --- /dev/null +++ b/src/MME.cpp @@ -0,0 +1,884 @@ +// solveUVM and solveMVM are based on R/EMMREML functions +// solveMKM is based on the mmer function in R/sommer +#include "alphamme.h" +#include +#include + +// 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 +// If calcVec = false, eigvec is not used +// It would be better to template this function +int eigen2(arma::vec& eigval, arma::mat& eigvec, arma::mat X, + bool calcVec = true){ + char JOBZ; + if(calcVec){ + JOBZ = 'V'; + }else{ + JOBZ = 'N'; + } + char RANGE = 'A'; + char UPLO = 'L'; + long long int N = X.n_rows; + // A = X + long long int LDA = N; + double VL = 0.0; + double VU = 0.0; + long long int IL = 0; + long long int IU = 0; + double ABSTOL = 0.0; + long long int M = N; + // W=eigval + // Z=eigvec + long long int LDZ = N; + arma::Col ISUPPZ(2*M); + // WORK length to be determined + double tmpWORK; + long long int LWORK = -1; // To be calculated + // IWORK length to be determined + long long int tmpIWORK = 0; + 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); + 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(), + &*eigvec.begin(),&LDZ,&*ISUPPZ.begin(),&*WORK.begin(),&LWORK,&*IWORK.begin(),&LIWORK,&INFO); + return INFO; // Return error code +} + +// Objective function for REML using the EMMA algorithm +Rcpp::List objREML(double param, Rcpp::List args){ + double df = args["df"]; + arma::vec eta = args["eta"]; + arma::vec lambda = args["lambda"]; + double value = df * log(sum(eta%eta/(lambda+param))); + value += sum(log(lambda+param)); + return Rcpp::List::create(Rcpp::Named("objective") = value, + Rcpp::Named("output") = 0); +} + +//' @title Read Matrix +//' +//' @description +//' Uses C++ to quickly read a matrix from a text +//' file. Requires knowledge of the number of rows +//' and columns in the file. +//' +//' @param fileName path to the file to 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 skipRows number of rows to skip +//' @param skipCols number of columns to skip +//' +//' @return a numeric matrix +//' +//' @export +// [[Rcpp::export]] +arma::mat readMat(std::string fileName, int rows, int cols, + char sep=' ', int skipRows=0, int skipCols=0){ + arma::mat output(rows,cols); + std::ifstream file(fileName.c_str()); + std::string line; + //Skip rows + for(arma::uword i=0; i0.0){ + numer = fabs(sum(VeNew.diag()-Ve.diag())); + if((numer/denom)=maxIter){ + Rf_warning("Reached maxIter without converging"); + break; + } + } + arma::mat HI = inv_sympd(kron(ZKZ, Vu)+ + kron(arma::eye(n,n), Ve)+ + tol*arma::eye(n*m,n*m)); + arma::mat E = Y.t() - B*X.t(); + arma::mat U = kron(K, Vu)*kron(Z.t(), + arma::eye(m,m))*(HI*vectorise(E)); //BLUPs + U.reshape(m,U.n_elem/m); + //Log Likelihood calculation + arma::mat ll = -0.5*arma::vectorise(E).t()*HI*vectorise(E); + ll -= double(n*m)/2.0*log(2*PI); + double value; + double sign; + log_det(value, sign, kron(ZKZ, Vu)+kron(arma::eye(n,n), Ve)); + ll -= 0.5*value*sign; + return Rcpp::List::create(Rcpp::Named("Vu")=Vu, + Rcpp::Named("Ve")=Ve, + Rcpp::Named("beta")=B.t(), + Rcpp::Named("u")=U.t(), + Rcpp::Named("LL")=arma::as_scalar(ll), + Rcpp::Named("iter")=iter); +} + +//' @title Solve Multivariate RR-BLUP +//' +//' @description +//' Solves a multivariate mixed model of form \eqn{Y=X\beta+Mu+e} +//' +//' @param Y a matrix with n rows and q columns +//' @param X a matrix with n rows and x columns +//' @param M a matrix with n rows and m columns +//' @param tol tolerance for convergence +//' @param maxIter maximum number of iteration +//' +//' @export +// [[Rcpp::export]] +Rcpp::List solveRRBLUPMV(const arma::mat& Y, const arma::mat& X, + const arma::mat& M, double tol=1e-6, + int maxIter=1000){ + int n = Y.n_rows; + int m = Y.n_cols; + arma::vec eigval(n); + arma::mat eigvec(n,n); + eigen2(eigval, eigvec, M*M.t()); + arma::mat Yt = Y.t()*eigvec; + arma::mat Xt = X.t()*eigvec; + arma::mat Vu = cov(Y)/2; + arma::mat Ve = Vu; + arma::mat W = Xt.t()*inv_sympd(Xt*Xt.t()); + arma::mat B = Yt*W; //BLUEs + arma::mat Gt(m,n), sigma(m,m), BNew, + VeNew(m,m), VuNew(m,m); + double denom, numer; + bool converging=true; + int iter=0; + while(converging){ + ++iter; + VeNew.fill(0.0); + VuNew.fill(0.0); + for(arma::uword i=0; i0.0){ + numer = fabs(sum(VeNew.diag()-Ve.diag())); + if((numer/denom)=maxIter){ + Rcpp::Rcerr<<"Warning: did not converge, reached maxIter\n"; + break; + } + } + arma::mat HI = inv_sympd(kron(M*M.t(), Vu)+ + kron(arma::eye(n,n), Ve)+ + tol*arma::eye(n*m,n*m)); + arma::mat E = Y.t() - B*X.t(); + arma::mat U = kron(arma::eye(M.n_cols,M.n_cols), Vu)*kron(M.t(), + arma::eye(m,m))*(HI*vectorise(E)); //BLUPs + U.reshape(m,U.n_elem/m); + //Log Likelihood calculation + arma::mat ll = -0.5*arma::vectorise(E).t()*HI*vectorise(E); + ll -= double(n*m)/2.0*log(2*PI); + double value; + double sign; + log_det(value, sign, kron(M*M.t(), Vu)+kron(arma::eye(n,n), Ve)); + ll -= 0.5*value*sign; + return Rcpp::List::create(Rcpp::Named("Vu")=Vu, + Rcpp::Named("Ve")=Ve, + Rcpp::Named("beta")=B.t(), + Rcpp::Named("u")=U.t(), + Rcpp::Named("LL")=arma::as_scalar(ll), + Rcpp::Named("iter")=iter); +} + +//' @title Solve Multikernel Model +//' +//' @description +//' Solves a univariate mixed model with multiple random effects. +//' +//' @param y a matrix with n rows and 1 column +//' @param X a matrix with n rows and x columns +//' @param Zlist a list of Z matrices +//' @param Klist a list of K matrices +//' @param maxIter maximum number of iteration +//' +//' @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 k = Klist.n_elem; + int n = y.n_rows; + int q = X.n_cols; + double df = double(n)-double(q); + arma::field V(k); + for(arma::uword i=0; i T(k); + sigma.fill(var(y.col(0))); + int iter = 0; + while(true){ + ++iter; + W0 = V(0)*sigma(0); + W0.diag() += sigma(k); + for(arma::uword i=1; i 1 & fabs(deltaLlik) < tol*10){ + break; + } + if(max(abs(qvec)) < tol){ + break; + } + if(iter >= maxIter){ + Rf_warning("Reached maxIter without converging"); + break; + } + } + while(sigma.min() < 0.0){ + sigma(sigma.index_min()) = 0.0; + } + arma::mat beta(q,1), ee(n,1); + arma::field u(k); + beta = solve(X.t()*W*X,X.t()*W*y); + ee = y - X*beta; + for(arma::uword i=0; i& Mlist, + int maxIter=40){ + double tol = 1e-4; + int k = Mlist.n_elem; + int n = y.n_rows; + int q = X.n_cols; + double df = double(n)-double(q); + arma::field V(k); + for(arma::uword i=0; i T(k); + sigma.fill(var(y.col(0))); + int iter = 0; + while(true){ + ++iter; + W0 = V(0)*sigma(0); + W0.diag() += sigma(k); + for(arma::uword i=1; i 1 & fabs(deltaLlik) < tol*10){ + break; + } + if(max(abs(qvec)) < tol){ + break; + } + if(iter >= maxIter){ + Rf_warning("Reached maxIter without converging"); + break; + } + } + while(sigma.min() < 0.0){ + sigma(sigma.index_min()) = 0.0; + } + arma::mat beta(q,1), ee(n,1); + arma::field u(k); + beta = solve(X.t()*W*X,X.t()*W*y); + ee = y - X*beta; + for(arma::uword i=0; i do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include +#include + +using namespace Rcpp; + +// readMat +arma::mat readMat(std::string fileName, int rows, int cols, char sep, int skipRows, int skipCols); +RcppExport SEXP _AlphaMME_readMat(SEXP fileNameSEXP, SEXP rowsSEXP, SEXP colsSEXP, SEXP sepSEXP, SEXP skipRowsSEXP, SEXP skipColsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< std::string >::type fileName(fileNameSEXP); + Rcpp::traits::input_parameter< int >::type rows(rowsSEXP); + Rcpp::traits::input_parameter< int >::type cols(colsSEXP); + Rcpp::traits::input_parameter< char >::type sep(sepSEXP); + Rcpp::traits::input_parameter< int >::type skipRows(skipRowsSEXP); + Rcpp::traits::input_parameter< int >::type skipCols(skipColsSEXP); + rcpp_result_gen = Rcpp::wrap(readMat(fileName, rows, cols, sep, skipRows, skipCols)); + return rcpp_result_gen; +END_RCPP +} +// solveUVM +Rcpp::List solveUVM(const arma::mat& y, const arma::mat& X, const arma::mat& Z, const arma::mat& K); +RcppExport SEXP _AlphaMME_solveUVM(SEXP ySEXP, SEXP XSEXP, SEXP ZSEXP, SEXP KSEXP) { +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 Z(ZSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type K(KSEXP); + rcpp_result_gen = Rcpp::wrap(solveUVM(y, X, Z, K)); + return rcpp_result_gen; +END_RCPP +} +// solveAniModel +Rcpp::List solveAniModel(const arma::mat& y, const arma::mat& X, const arma::mat& K); +RcppExport SEXP _AlphaMME_solveAniModel(SEXP ySEXP, SEXP XSEXP, SEXP KSEXP) { +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 K(KSEXP); + rcpp_result_gen = Rcpp::wrap(solveAniModel(y, X, K)); + return rcpp_result_gen; +END_RCPP +} +// solveRRBLUP +Rcpp::List solveRRBLUP(const arma::mat& y, const arma::mat& X, const arma::mat& M); +RcppExport SEXP _AlphaMME_solveRRBLUP(SEXP ySEXP, SEXP XSEXP, SEXP MSEXP) { +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_result_gen = Rcpp::wrap(solveRRBLUP(y, X, M)); + return rcpp_result_gen; +END_RCPP +} +// solveMVM +Rcpp::List solveMVM(const arma::mat& Y, const arma::mat& X, const arma::mat& Z, const arma::mat& K, double tol, int maxIter); +RcppExport SEXP _AlphaMME_solveMVM(SEXP YSEXP, SEXP XSEXP, SEXP ZSEXP, SEXP KSEXP, SEXP tolSEXP, SEXP maxIterSEXP) { +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 Z(ZSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type K(KSEXP); + Rcpp::traits::input_parameter< double >::type tol(tolSEXP); + Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP); + rcpp_result_gen = Rcpp::wrap(solveMVM(Y, X, Z, K, tol, maxIter)); + return rcpp_result_gen; +END_RCPP +} +// solveRRBLUPMV +Rcpp::List solveRRBLUPMV(const arma::mat& Y, const arma::mat& X, const arma::mat& M, double tol, int maxIter); +RcppExport SEXP _AlphaMME_solveRRBLUPMV(SEXP YSEXP, SEXP XSEXP, SEXP MSEXP, SEXP tolSEXP, SEXP maxIterSEXP) { +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 tol(tolSEXP); + Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP); + rcpp_result_gen = Rcpp::wrap(solveRRBLUPMV(Y, X, M, tol, maxIter)); + return rcpp_result_gen; +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) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat& >::type y(ySEXP); + Rcpp::traits::input_parameter< arma::mat& >::type X(XSEXP); + 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)); + return rcpp_result_gen; +END_RCPP +} +// solveRRBLUPMK +Rcpp::List solveRRBLUPMK(arma::mat& y, arma::mat& X, arma::field& Mlist, int maxIter); +RcppExport SEXP _AlphaMME_solveRRBLUPMK(SEXP ySEXP, SEXP XSEXP, SEXP MlistSEXP, SEXP maxIterSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat& >::type y(ySEXP); + Rcpp::traits::input_parameter< arma::mat& >::type X(XSEXP); + Rcpp::traits::input_parameter< arma::field& >::type Mlist(MlistSEXP); + Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP); + rcpp_result_gen = Rcpp::wrap(solveRRBLUPMK(y, X, Mlist, maxIter)); + return rcpp_result_gen; +END_RCPP +} +// calcG +arma::mat calcG(arma::mat X); +RcppExport SEXP _AlphaMME_calcG(SEXP XSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP); + rcpp_result_gen = Rcpp::wrap(calcG(X)); + return rcpp_result_gen; +END_RCPP +} +// calcD +arma::mat calcD(arma::mat X); +RcppExport SEXP _AlphaMME_calcD(SEXP XSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP); + rcpp_result_gen = Rcpp::wrap(calcD(X)); + return rcpp_result_gen; +END_RCPP +} +// calcGIbs +arma::mat calcGIbs(arma::mat X); +RcppExport SEXP _AlphaMME_calcGIbs(SEXP XSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP); + rcpp_result_gen = Rcpp::wrap(calcGIbs(X)); + return rcpp_result_gen; +END_RCPP +} +// fastDist +arma::mat fastDist(const arma::mat& X); +RcppExport SEXP _AlphaMME_fastDist(SEXP XSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); + rcpp_result_gen = Rcpp::wrap(fastDist(X)); + return rcpp_result_gen; +END_RCPP +} +// fastPairDist +arma::mat fastPairDist(const arma::mat& X, const arma::mat& Y); +RcppExport SEXP _AlphaMME_fastPairDist(SEXP XSEXP, SEXP YSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + rcpp_result_gen = Rcpp::wrap(fastPairDist(X, Y)); + return rcpp_result_gen; +END_RCPP +} +// gaussKernel +arma::mat gaussKernel(arma::mat& D, double theta); +RcppExport SEXP _AlphaMME_gaussKernel(SEXP DSEXP, SEXP thetaSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat& >::type D(DSEXP); + Rcpp::traits::input_parameter< double >::type theta(thetaSEXP); + rcpp_result_gen = Rcpp::wrap(gaussKernel(D, theta)); + return rcpp_result_gen; +END_RCPP +} diff --git a/src/alphamme.h b/src/alphamme.h new file mode 100644 index 0000000..fbb43b8 --- /dev/null +++ b/src/alphamme.h @@ -0,0 +1,7 @@ +#ifndef ALPHAMME_H +#define ALPHAMME_H + +#include +#include "optimize.h" + +#endif \ No newline at end of file diff --git a/src/init.c b/src/init.c new file mode 100644 index 0000000..901e43e --- /dev/null +++ b/src/init.c @@ -0,0 +1,48 @@ +#include +#include +#include // for NULL +#include + +/* FIXME: + Check these declarations against the C/Fortran source code. +*/ + +/* .Call calls */ +extern SEXP _AlphaMME_calcD(SEXP); +extern SEXP _AlphaMME_calcG(SEXP); +extern SEXP _AlphaMME_calcGIbs(SEXP); +extern SEXP _AlphaMME_fastDist(SEXP); +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_solveMVM(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _AlphaMME_solveRRBLUP(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}, + {NULL, NULL, 0} +}; + +void R_init_AlphaMME(DllInfo *dll) +{ + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/optimize.cpp b/src/optimize.cpp new file mode 100644 index 0000000..aafc0fd --- /dev/null +++ b/src/optimize.cpp @@ -0,0 +1,155 @@ +#include +using namespace Rcpp; +/* + * Modified version of Brent's method for numerical optimization + * objective is a pointer to a function returning a List + * objective takes a double and a List as arguments + * optimization is for the double + * the returned list includes "objective" and list of additional output + * args is a list which is passed to the objective + * l is the lower bound of optimization + * u is the upper bound of optimization + * maxIter is the maximum number of iterations + * the best obtained value is returned when maxIter is reached + * if firstLast is true, l and u are evaluated if convergence isn't reached + * 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, + bool evalU=false, bool evalL=false, double eps=1.0e-9){ + double MACHEPS_SQRT = sqrt(2.2204460492503131e-016); + double c = (3-sqrt(5))/2; + double x = l+c*(u-l); + double v = x; + double w = x; + double e = 0; + double lInt = l; + double uInt = u; + List fOut = objective(x, args); + List output = fOut["output"]; + double fx = fOut["objective"]; + 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); + q = (x-v)*(fx-fw); + p = (x-v)*q-(x-w)*r; + q = 2.0*(q-r); + if(q>0.0){ + p = -p; + }else{ + q = -q; + } + r = e; + e = d; + } + + if((std::abs(p)=tol){ + z = x+d; + }else if(d>0.0){ + z = x+tol; + }else{ + z = x-tol; + } + fOut = objective(z,args); + double funcu = fOut["objective"]; + if(maximize) funcu = -funcu; + + // Update + if(funcu<=fx){ + if(z