0

I have a package on CRAN called 'metapack' and received a message from Professor Ripley that the error must be resolved for the package to be safely retained in CRAN. I've encountered this error on R-CMD-check on GitHub but cannot replicate it on any of my local machines to save my life. The reproducible code snippet is as follows (I use Rcpp/RcppArmadillo):

library(metapack)
data("cholesterol")
Outcome <- model.matrix(~ pldlc + phdlc + ptg, data = cholesterol)
SD <- model.matrix(~ sdldl + sdhdl + sdtg, data = cholesterol)
Trial <- cholesterol$Trial
Treat <- cholesterol$trt
Npt <- cholesterol$Npt
XCovariate <- model.matrix(~ 0 + bldlc + bhdlc + btg + age + durat + white + male + dm, data = cholesterol)
WCovariate <- model.matrix(~ trt, data = cholesterol)

fmodel <- 1
set.seed(2797542)
fit <- bayes.parobs(Outcome, SD, XCovariate, WCovariate, Treat, Trial,
  Npt, fmodel, mcmc = list(ndiscard = 1, nskip = 1, nkeep = 1),
  scale_x = TRUE, group = cholesterol$onstat, verbose = FALSE)

which gives me the following error:

warning: chol(): given matrix is not symmetric  
error: chol(): decomposition failed

This can't be true because (1) it runs without problem on other machines, (2) every matrix that is Cholesky-decomposed is ensured to be symmetric through Sig = 0.5 * (Sig + Sig.t());. Even with Sig = arma::symmatu(Sig) doesn't seem to guarantee symmetry.

Has anyone encountered something similar to this?

Daeyoung
  • 194
  • 7
  • 1
    Sorry to hear you are having issues. Reproducing things on _that machine_ can be a challenge. Did you try RHub's Solaris instance? (And please don't cross-post here and on r-package-devel where I just moderated your message in (as you were apparently a first-time poster needed moderation).) – Dirk Eddelbuettel Feb 18 '21 at 00:57
  • @DirkEddelbuettel oops. I won’t cross post. I wanted to get as much as possible. I have tried RHub’s Solaris instance and it ran without any warnings other than the size of the package. – Daeyoung Feb 18 '21 at 01:00
  • 1
    That's the worst case. Because now you can't reproduce. I have resorted to skipping tests on problematic machines. Not ideal, but better than getting the package thrown off CRAN. – Dirk Eddelbuettel Feb 18 '21 at 02:00
  • @DirkEddelbuettel If it's not rude to ask you this here (I don't want to bother you here), but that is actually an example code in the vignette, requested by CRAN maintainers. I don't know the formal route to let them know (or argue) that it should be skipped on Solaris. I'm sure they don't want me to email. – Daeyoung Feb 18 '21 at 02:34
  • 3
    Sorry I was between things. So basically in a vignette you can simple _precompute_ the chart or table and just referene it avoiding all the trouble. In unit tests I sometime just as 'am I currently running on operating sytem X' and then do things like skip a test or the remainder of the file. Dirty, but the best we can do. Look for example at my `anytime` package for how it is done (in the context of `tinytest` but I am sure you can adjust accordingly). – Dirk Eddelbuettel Feb 18 '21 at 02:58
  • 4
    You have a valgrind error in that exact spot: https://www.stats.ox.ac.uk/pub/bdr/memtests/valgrind/metapack/metapack-Ex.Rout. So it would be better to debug that instead of just skipping that on Solaris. Looking at the trace, it may be a bug in Armadillo that it's worth reporting. – Iñaki Úcar Feb 18 '21 at 09:37
  • 5
    @DaeyoungLim Following the valgrind log, the error is probably [here](https://github.com/cran/metapack/blob/c964d0bf04137d4ce1e9745cdf3a4b6aa4f8761d/src/random.cpp#L475-L481). The matrix is only partly initlalized. Try changing `arma::mat A(p,p)` to `arma::mat A(p,p,arma::fill:zeros)`. Also, why make your own `RNG::rwish()` and `RNG::riwish()`? Armadillo already has [wishrnd()](http://arma.sourceforge.net/docs.html#wishrnd) and [iwishrnd()](http://arma.sourceforge.net/docs.html#iwishrnd). – mtall Feb 18 '21 at 13:36
  • @mtall regarding my own samplers, as far as I know, Armadillo does not use the same seed as R. – Daeyoung Feb 18 '21 at 13:51
  • 1
    @DaeyoungLim RcppArmadillo appears to use the same seed as R. This gives strong hints: https://github.com/RcppCore/RcppArmadillo/blob/master/inst/include/RcppArmadillo/Alt_R_RNG.h – mtall Feb 22 '21 at 11:15

0 Answers0