Getting started with the PCMBaseCpp R-package

2020-03-22

How to use the package?

There are two ways to use PCMBaseCpp:

Passing the function PCMInfoCpp as a metaI argument of PCMLik and/or PCMCreateLikelihood

library(PCMBase)
library(PCMBaseCpp)

## Loading required package: Rcpp

system.time(llR <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab))

##    user  system elapsed 
##   0.088   0.001   0.090

system.time(llCpp <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = PCMInfoCpp))

##    user  system elapsed 
##   0.006   0.000   0.005

print(llR)

## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1

print(llCpp)

## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1

logLikFunR <- PCMCreateLikelihood(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab)

logLikFunCpp <- PCMCreateLikelihood(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = PCMInfoCpp)

metaICpp <- PCMInfoCpp(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab)

logLikFunCpp2 <- PCMCreateLikelihood(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaICpp)

set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
randParam <- PCMParamRandomVecParams(PCMBaseTestObjects$model_MixedGaussian_ab)

system.time(llR <- logLikFunR(randParam))

##    user  system elapsed 
##   0.069   0.000   0.070

system.time(llCpp <- logLikFunCpp(randParam))

##    user  system elapsed 
##   0.002   0.000   0.002

system.time(llCpp2 <- logLikFunCpp2(randParam))

##    user  system elapsed 
##   0.002   0.000   0.002

print(llR)

## [1] -598.092
## attr(,"X0")
## [1] -4.689827 -2.557522  1.457067

print(llCpp)

## [1] -598.092
## attr(,"X0")
## [1] -4.689827 -2.557522  1.457067

print(llCpp2)

## [1] -598.092
## attr(,"X0")
## [1] -4.689827 -2.557522  1.457067

Passing the meta-information object returned by PCMInfoCpp as a metaI argument of PCMLik and PCMCreateLikelihood

This is the recommended usage in the case of multiple likelihood evaluations, e.g. during model inference:

metaIR <- PCMInfo(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab)

metaICpp <- PCMInfoCpp(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab)

system.time(llR <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaIR))

##    user  system elapsed 
##    0.07    0.00    0.07

system.time(llCpp <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaICpp))

##    user  system elapsed 
##   0.001   0.000   0.001

print(llR)

## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1

print(llCpp)

## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1

logLikFunR <- PCMCreateLikelihood(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaIR)

logLikFunCpp <- PCMCreateLikelihood(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaICpp)

system.time(llR <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaIR))

##    user  system elapsed 
##   0.068   0.000   0.068

system.time(llCpp <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaICpp))

##    user  system elapsed 
##   0.001   0.000   0.001

print(llR)

## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1

print(llCpp)

## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1