Skip to content

Commit

Permalink
Use AbcPMM everywhere, e.g. in function PMMLogLik use AbcPMMLogLik. T…
Browse files Browse the repository at this point in the history 
…his is necessary to facilitate creating a project using PMMUsingSPLITT as a template.
  • Loading branch information
Venelin Mitov committed Nov 18, 2018
1 parent 0a2ffbe commit 26ef35f
Show file tree
Hide file tree
Showing 14 changed files with 677 additions and 125 deletions.
22 changes: 11 additions & 11 deletions R/PMM.R → R/AbcPMM.R
View file
@@ -1,4 +1,4 @@
# PMM.R
# AbcPMM.R
# PMMUsingSPLITT
#
# Copyright 2018 Venelin Mitov
Expand Down Expand Up @@ -37,7 +37,7 @@
#' which is also set as a default value. Can be passed as argument to speed-up
#' the calculation.
#' @return the log-likelihood value.
PMMLogLik <- function(
AbcPMMLogLik <- function(
x, tree, x0, sigma2, sigmae2,
ord = reorder(tree, order = "postorder", index.only = TRUE)) {

Expand Down Expand Up @@ -85,14 +85,14 @@ PMMLogLik <- function(
a[N+1]*x0^2 + b[N+1]*x0 + c[N+1]
}

#' Calculate the PMM log-likelihood for a given tree, data and model parameters using the Rcpp module
#' @inheritParams PMMLogLik
#' @param cppObject a previously created object returned by \code{\link{NewPMMCppObject}}
#' Calculate the AbcPMM log-likelihood for a given tree, data and model parameters using the Rcpp module
#' @inheritParams AbcPMMLogLik
#' @param cppObject a previously created object returned by \code{\link{NewAbcPMMCppObject}}
#' @param mode an integer denoting the mode for traversing the tree, i.e. serial vs parallel.
#'
#' @return the log-likelihood value.
PMMLogLikCpp <- function(x, tree, x0, sigma2, sigmae2,
cppObject = NewPMMCppObject(x, tree),
AbcPMMLogLikCpp <- function(x, tree, x0, sigma2, sigmae2,
cppObject = NewAbcPMMCppObject(x, tree),
mode = getOption("SPLITT.postorder.mode", 0)) {
abc <- cppObject$TraverseTree(c(sigma2, sigmae2), mode)
abc[1]*x0^2 + abc[2]*x0 + abc[3]
Expand All @@ -101,9 +101,9 @@ PMMLogLikCpp <- function(x, tree, x0, sigma2, sigmae2,

#' Create an instance of the Rcpp module for a given tree and trait data
#'
#' @inheritParams PMMLogLik
#' @return an object to be passed as argument of the \link{PMMLogLikCpp} function.
#' @seealso \link{PMMLogLikCpp}
NewPMMCppObject <- function(x, tree) {
#' @inheritParams AbcPMMLogLik
#' @return an object to be passed as argument of the \link{AbcPMMLogLikCpp} function.
#' @seealso \link{AbcPMMLogLikCpp}
NewAbcPMMCppObject <- function(x, tree) {
PMMUsingSPLITT__TraversalTaskAbcPMM$new(tree, x[1:length(tree$tip.label)])
}
40 changes: 20 additions & 20 deletions R/MiniBenchmark.R
View file
Expand Up @@ -70,7 +70,7 @@ MiniBenchmark <- function(N = 10000, Ntests = 10) {

ord <- reorder(tree, order = "postorder", index.only = TRUE)

cppPMMObject <- NewPMMCppObject(x, tree)
cppPMMObject <- NewAbcPMMCppObject(x, tree)

if(cppPMMObject$algorithm$VersionOPENMP == 0) {
cat("It seems that OpenMP was disabled during C++ compilation. For parallel
Expand All @@ -87,63 +87,63 @@ MiniBenchmark <- function(N = 10000, Ntests = 10) {
cat("Measuring calculation times...\n")

# warm-up for mode AUTO
for(t in 1:1000) PMMLogLikCpp(x, tree, x0, sigma2, sigmae2, cppPMMObject, 0)
for(t in 1:1000) AbcPMMLogLikCpp(x, tree, x0, sigma2, sigmae2, cppPMMObject, 0)

tR <- system.time(for(t in seq_len(Ntests/10)) PMMLogLik(x, tree, x0, sigma2, sigmae2, ord = ord))[3] / (Ntests/10)
tR <- system.time(for(t in seq_len(Ntests/10)) AbcPMMLogLik(x, tree, x0, sigma2, sigmae2, ord = ord))[3] / (Ntests/10)
unname(tR)

measureTimePMMCpp <- function(mode) {
measureTimeAbcPMMCpp <- function(mode) {
unname(
system.time(
for(t in seq_len(Ntests))
PMMLogLikCpp(x, tree, x0, sigma2, sigmae2, cppPMMObject, mode)
AbcPMMLogLikCpp(x, tree, x0, sigma2, sigmae2, cppPMMObject, mode)
)[3] / Ntests*1000
)
}

resultsPMM <- rbind(
resultsAbcPMM <- rbind(
data.frame(model = "PMM",
mode = "R (serial)",
time.ms = tR*1000),
data.frame(model = "PMM",
mode = "C++ (AUTO)",
time.ms = measureTimePMMCpp(0)),
time.ms = measureTimeAbcPMMCpp(0)),
data.frame(model = "PMM",
mode = "C++ (SINGLE_THREAD_LOOP_POSTORDER)",
time.ms = measureTimePMMCpp(10)),
time.ms = measureTimeAbcPMMCpp(10)),
data.frame(model = "PMM",
mode = "C++ (SINGLE_THREAD_LOOP_PRUNES)",
time.ms = measureTimePMMCpp(11)),
time.ms = measureTimeAbcPMMCpp(11)),
data.frame(model = "PMM",
mode = "C++ (SINGLE_THREAD_LOOP_VISITS)",
time.ms = measureTimePMMCpp(12)),
time.ms = measureTimeAbcPMMCpp(12)),
data.frame(model = "PMM",
mode = "C++ (MULTI_THREAD_LOOP_PRUNES)",
time.ms = measureTimePMMCpp(21)),
time.ms = measureTimeAbcPMMCpp(21)),
data.frame(model = "PMM",
mode = "C++ (MULTI_THREAD_LOOP_VISITS)",
time.ms = measureTimePMMCpp(22)),
time.ms = measureTimeAbcPMMCpp(22)),
data.frame(model = "PMM",
mode = "C++ (MULTI_THREAD_LOOP_VISITS_THEN_LOOP_PRUNES)",
time.ms = measureTimePMMCpp(23)),
time.ms = measureTimeAbcPMMCpp(23)),
data.frame(model = "PMM",
mode = "C++ (MULTI_THREAD_VISIT_QUEUE)",
time.ms = measureTimePMMCpp(24)),
time.ms = measureTimeAbcPMMCpp(24)),
data.frame(model = "PMM",
mode = "C++ (MULTI_THREAD_LOOP_PRUNES_NO_EXCEPTION)",
time.ms = measureTimePMMCpp(25)),
time.ms = measureTimeAbcPMMCpp(25)),
data.frame(model = "PMM",
mode = "C++ (HYBRID_LOOP_PRUNES)",
time.ms = measureTimePMMCpp(31)),
time.ms = measureTimeAbcPMMCpp(31)),
data.frame(model = "PMM",
mode = "C++ (HYBRID_LOOP_VISITS)",
time.ms = measureTimePMMCpp(32)),
time.ms = measureTimeAbcPMMCpp(32)),
data.frame(model = "PMM",
mode = "C++ (HYBRID_LOOP_VISITS_THEN_LOOP_PRUNES)",
time.ms = measureTimePMMCpp(33))
time.ms = measureTimeAbcPMMCpp(33))
)

rownames(resultsPMM) <- NULL
rownames(resultsAbcPMM) <- NULL

resultsPMM
resultsAbcPMM
}
4 changes: 2 additions & 2 deletions README.Rmd
View file
Expand Up @@ -55,8 +55,8 @@ g <- rTraitCont(tree, model = "OU", root.value = x0,
x <- g + rnorm(n = N, mean = 0, sd = sqrt(sigmae2))
cat("logLikelihood using R:", PMMLogLik(x, tree, x0, sigma2, sigmae2), "\n")
cat("logLikelihood using R:", PMMLogLikCpp(x, tree, x0, sigma2, sigmae2), "\n")
cat("logLikelihood using R:", AbcPMMLogLik(x, tree, x0, sigma2, sigmae2), "\n")
cat("logLikelihood using R:", AbcPMMLogLikCpp(x, tree, x0, sigma2, sigmae2), "\n")
```

* Performing a benchmark to measure the likelihood calculation times using different parallelization strategies:
Expand Down
54 changes: 27 additions & 27 deletions README.md
View file
Expand Up @@ -43,9 +43,9 @@ g <- rTraitCont(tree, model = "OU", root.value = x0,

x <- g + rnorm(n = N, mean = 0, sd = sqrt(sigmae2))

cat("logLikelihood using R:", PMMLogLik(x, tree, x0, sigma2, sigmae2), "\n")
cat("logLikelihood using R:", AbcPMMLogLik(x, tree, x0, sigma2, sigmae2), "\n")
#> logLikelihood using R: -1551.016
cat("logLikelihood using R:", PMMLogLikCpp(x, tree, x0, sigma2, sigmae2), "\n")
cat("logLikelihood using R:", AbcPMMLogLikCpp(x, tree, x0, sigma2, sigmae2), "\n")
#> logLikelihood using R: -1551.016
```

Expand All @@ -64,19 +64,19 @@ MiniBenchmark(N = 100, Ntests = 100)
#> Number of threads: 8
#> Measuring calculation times...
#> model mode time.ms
#> 1 PMM R (serial) 2.00
#> 1 PMM R (serial) 2.10
#> 2 PMM C++ (AUTO) 0.03
#> 3 PMM C++ (SINGLE_THREAD_LOOP_POSTORDER) 0.01
#> 4 PMM C++ (SINGLE_THREAD_LOOP_PRUNES) 0.02
#> 3 PMM C++ (SINGLE_THREAD_LOOP_POSTORDER) 0.02
#> 4 PMM C++ (SINGLE_THREAD_LOOP_PRUNES) 0.01
#> 5 PMM C++ (SINGLE_THREAD_LOOP_VISITS) 0.02
#> 6 PMM C++ (MULTI_THREAD_LOOP_PRUNES) 0.05
#> 7 PMM C++ (MULTI_THREAD_LOOP_VISITS) 0.08
#> 6 PMM C++ (MULTI_THREAD_LOOP_PRUNES) 0.08
#> 7 PMM C++ (MULTI_THREAD_LOOP_VISITS) 0.06
#> 8 PMM C++ (MULTI_THREAD_LOOP_VISITS_THEN_LOOP_PRUNES) 0.10
#> 9 PMM C++ (MULTI_THREAD_VISIT_QUEUE) 1.69
#> 10 PMM C++ (MULTI_THREAD_LOOP_PRUNES_NO_EXCEPTION) 0.25
#> 11 PMM C++ (HYBRID_LOOP_PRUNES) 0.14
#> 12 PMM C++ (HYBRID_LOOP_VISITS) 0.24
#> 13 PMM C++ (HYBRID_LOOP_VISITS_THEN_LOOP_PRUNES) 0.38
#> 9 PMM C++ (MULTI_THREAD_VISIT_QUEUE) 1.64
#> 10 PMM C++ (MULTI_THREAD_LOOP_PRUNES_NO_EXCEPTION) 0.05
#> 11 PMM C++ (HYBRID_LOOP_PRUNES) 0.08
#> 12 PMM C++ (HYBRID_LOOP_VISITS) 0.12
#> 13 PMM C++ (HYBRID_LOOP_VISITS_THEN_LOOP_PRUNES) 0.14
MiniBenchmark(N = 1000, Ntests = 10)
#> Performing a mini-benchmark of the PMM log-likelihood calculation with
#> a tree of size N= 1000 ;
Expand All @@ -86,19 +86,19 @@ MiniBenchmark(N = 1000, Ntests = 10)
#> Number of threads: 8
#> Measuring calculation times...
#> model mode time.ms
#> 1 PMM R (serial) 20.0
#> 1 PMM R (serial) 19.0
#> 2 PMM C++ (AUTO) 0.1
#> 3 PMM C++ (SINGLE_THREAD_LOOP_POSTORDER) 0.1
#> 3 PMM C++ (SINGLE_THREAD_LOOP_POSTORDER) 0.2
#> 4 PMM C++ (SINGLE_THREAD_LOOP_PRUNES) 0.1
#> 5 PMM C++ (SINGLE_THREAD_LOOP_VISITS) 0.1
#> 6 PMM C++ (MULTI_THREAD_LOOP_PRUNES) 0.2
#> 7 PMM C++ (MULTI_THREAD_LOOP_VISITS) 0.1
#> 5 PMM C++ (SINGLE_THREAD_LOOP_VISITS) 0.2
#> 6 PMM C++ (MULTI_THREAD_LOOP_PRUNES) 0.1
#> 7 PMM C++ (MULTI_THREAD_LOOP_VISITS) 0.5
#> 8 PMM C++ (MULTI_THREAD_LOOP_VISITS_THEN_LOOP_PRUNES) 0.2
#> 9 PMM C++ (MULTI_THREAD_VISIT_QUEUE) 16.3
#> 9 PMM C++ (MULTI_THREAD_VISIT_QUEUE) 17.8
#> 10 PMM C++ (MULTI_THREAD_LOOP_PRUNES_NO_EXCEPTION) 0.6
#> 11 PMM C++ (HYBRID_LOOP_PRUNES) 0.6
#> 12 PMM C++ (HYBRID_LOOP_VISITS) 0.6
#> 13 PMM C++ (HYBRID_LOOP_VISITS_THEN_LOOP_PRUNES) 0.9
#> 11 PMM C++ (HYBRID_LOOP_PRUNES) 0.1
#> 12 PMM C++ (HYBRID_LOOP_VISITS) 0.5
#> 13 PMM C++ (HYBRID_LOOP_VISITS_THEN_LOOP_PRUNES) 0.1
MiniBenchmark(N = 10000, Ntests = 10)
#> Performing a mini-benchmark of the PMM log-likelihood calculation with
#> a tree of size N= 10000 ;
Expand All @@ -113,12 +113,12 @@ MiniBenchmark(N = 10000, Ntests = 10)
#> 3 PMM C++ (SINGLE_THREAD_LOOP_POSTORDER) 1.1
#> 4 PMM C++ (SINGLE_THREAD_LOOP_PRUNES) 1.1
#> 5 PMM C++ (SINGLE_THREAD_LOOP_VISITS) 1.2
#> 6 PMM C++ (MULTI_THREAD_LOOP_PRUNES) 0.3
#> 6 PMM C++ (MULTI_THREAD_LOOP_PRUNES) 1.4
#> 7 PMM C++ (MULTI_THREAD_LOOP_VISITS) 0.3
#> 8 PMM C++ (MULTI_THREAD_LOOP_VISITS_THEN_LOOP_PRUNES) 0.7
#> 9 PMM C++ (MULTI_THREAD_VISIT_QUEUE) 154.0
#> 10 PMM C++ (MULTI_THREAD_LOOP_PRUNES_NO_EXCEPTION) 0.2
#> 11 PMM C++ (HYBRID_LOOP_PRUNES) 0.3
#> 8 PMM C++ (MULTI_THREAD_LOOP_VISITS_THEN_LOOP_PRUNES) 0.4
#> 9 PMM C++ (MULTI_THREAD_VISIT_QUEUE) 151.0
#> 10 PMM C++ (MULTI_THREAD_LOOP_PRUNES_NO_EXCEPTION) 0.3
#> 11 PMM C++ (HYBRID_LOOP_PRUNES) 0.7
#> 12 PMM C++ (HYBRID_LOOP_VISITS) 0.2
#> 13 PMM C++ (HYBRID_LOOP_VISITS_THEN_LOOP_PRUNES) 0.3
#> 13 PMM C++ (HYBRID_LOOP_VISITS_THEN_LOOP_PRUNES) 0.4
```

0 comments on commit 26ef35f

Please sign in to comment.