Parallel Version/Thread Safe version

[mattfidler]

Hi @lbelzile,

Excellent package. For some reason I looked at cran for this package and got the wrong version:

https://github.com/cran/TruncatedNormal/blob/7364c5bc3f7c84d00eb4a767807b103b4232b648/R/ntail.R

This version is R only.

From there, I ported the package to a thread-safe/multi-threaded truncated normal distribution in RxODE using the threefry sitmo engine.

I noticed some deviations from the original version in this version, including truncated t-distribution.

Should we combine efforts? My main concern is having thread-safe routines for our threaded ode solving library RxODE. This would mean C level interfaces with all the routines.

[mattfidler]

See:

https://github.com/nlmixrdevelopment/RxODE/blob/pruneBranch/src/threefry.cpp

[lbelzile]

Thanks for reaching out, Matthew. The CRAN read-only versions are indeed not updated.

Zdravko Botev is the main author and I merely modified his original package (1.0) to be more R friendly and to include the truncated multivariate Student by porting Matlab code; I also maintain the R package for convenience.

I plan at some point to port the code to Rcpp (but C could be another option) because the current version is not fast enough to be competitive with existing packages that implement the separation-of-variables methods, so your proposal is interesting. I am overwhelmed at the moment, but can try to find some time to look into this after the academic semester ends.

[mattfidler]

The code I use is in RcppArmadillo; Here is the timing compared to tmvtnorm::rtmvnorm

In this example threading doesn't help too much:

mu <- c(0.5, 0.5)
sigma <- matrix(c(1, 0.8, 0.8, 2), 2, 2)
a <- c(-1, -Inf)
b <- c(0.5, 4)

mb <- microbenchmark::microbenchmark(tmvtnorm=tmvtnorm::rtmvnorm(n=10000, mean=mu, sigma=sigma, lower=a, upper=b),
                                                          rx1=RxODE::rxRmvn(n=10000,mu,sigma,a,b),
                                                          rx4=RxODE::rxRmvn(n=10000,mu,sigma,a,b,ncores=4))

print(mb)
#> Unit: milliseconds
#>      expr       min        lq      mean    median        uq       max neval
#>  tmvtnorm 27.107058 29.607081 31.881322 30.272680 30.950256 171.59679   100
#>       rx1  4.911288  5.036971 12.071839  5.122539  5.339175 688.88775   100
#>       rx4  5.024147  5.132115  5.375059  5.249840  5.445955   8.85388   100

plot(mb,log="y")

^{Created on 2019-12-06 by the reprex package (v0.3.0)}