nimble-dev / nimble

The base NIMBLE package for R
http://R-nimble.org
BSD 3-Clause "New" or "Revised" License
156 stars 23 forks source link

Failing LKJ test model when buildDerivs = TRUE #1403

Closed danielturek closed 8 months ago

danielturek commented 8 months ago

@perrydv @paciorek

There are two tests using a model with the dlkj_corr_cholesky distribution. Different versions of this model appear in the testing files of both the nimble package (in test-mcmc.R) and in the nimbleHMC package (in test-HMC.R).

In the nimble version of this test, the nimbleModel call does not build derivatives. In this case, the model builds and compiles fine.

In the nimbleHMC version of this test, the nimbleModel call includes buildDerivs = TRUE. Subsequently the model is built, and compilation of the model fails.

Reproducible code is provided below, which includes the buildDerivs = TRUE for the model, and fails at model compilation.

Noting also, this model involves a custom function (uppertri_mult_diag) which also has buildDerivs = TRUE, which could also be the source of a problem.

library(nimble)

R <- matrix(c(
    1, 0.9, .3, -.5, .1,
    0.9, 1, .15, -.3, .1,
    .3, .15, 1, .3, .1,
    -.5, -.3, .3, 1, .1,
    .1,.1,.1,.1, 1)
  , 5, 5)
U <- chol(R)
sds <- c(5, 4, 3, 2, 1)
set.seed(1)
Sigma <- diag(sds)%*%R%*%diag(sds)
n <- 100
p <- 5
y <- t(t(chol(Sigma))%*%matrix(rnorm(p*n),p,n))

uppertri_mult_diag <- nimbleFunction(
    run = function(mat = double(2), vec = double(1)) {
        returnType(double(2))
        p <- length(vec)
        out <- matrix(nrow = p, ncol = p, init = FALSE)
        for(i in 1:p)
            out[ , i] <- mat[ , i] * vec[i]
        return(out)
    },
    buildDerivs = list(run = list(ignore = 'i'))
)

code <- nimbleCode({
    for(i in 1:n)
        y[i, 1:p] ~ dmnorm(mu[1:p], cholesky = U[1:p, 1:p], prec_param = 0)
    U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
    Ustar[1:p,1:p] ~ dlkj_corr_cholesky(1.3, p)
})

m <- nimbleModel(code, constants = list(n = n, p = p, mu = rep(0, p)),
                 data = list(y = y), inits = list(sds = sds, Ustar = U))
                 ##buildDerivs = TRUE)

cm <- compileNimble(m, showCompilerOutput = TRUE)    ## FAILS

Error: Failed to create the shared library.
clang++ -mmacosx-version-min=10.13 -std=gnu++14 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -DR_NO_REMAP   -DEIGEN_MPL2_ONLY=1 -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/nimble/include" -Wno-misleading-indentation -Wno-ignored-attributes -Wno-deprecated-declarations  -I/usr/local/include   -fPIC  -Wall -g -O2  -c P_1_code_MID_1.cpp -o P_1_code_MID_1.o
clang++ -mmacosx-version-min=10.13 -std=gnu++14 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -DR_NO_REMAP   -DEIGEN_MPL2_ONLY=1 -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/nimble/include" -Wno-misleading-indentation -Wno-ignored-attributes -Wno-deprecated-declarations  -I/usr/local/include   -fPIC  -Wall -g -O2  -c P_1_code_MID_1_nfCode.cpp -o P_1_code_MID_1_nfCode.o
P_1_code_MID_1_nfCode.cpp:412:5: warning: variable 'Interm_25' set but not used [-Wunused-but-set-variable]
int Interm_25;
    ^
P_1_code_MID_1_nfCode.cpp:413:5: warning: variable 'Interm_26' set but not used [-Wunused-but-set-variable]
int Interm_26;
    ^
P_1_code_MID_1_nfCode.cpp:457:5: warning: variable 'Interm_25' set but not used [-Wunused-but-set-variable]
int Interm_25;
    ^
P_1_code_MID_1_nfCode.cpp:458:5: warning: variable 'Interm_26' set but not used [-Wunused-but-set-variable]
int Interm_26;
    ^
P_1_code_MID_1_nfCode.cpp:468:67: error: no matching function for call to 'stoch_ind_set'
 (eigenBlock_oPout_comma_1toInterm_26_comma_i_cP).block(0, i - 1, stoch_ind_set(out.dim(), 0), 1) = (Eig_ARG1_mat_Interm_27).block(0, i - 1, stoch_ind_get(ARG1_mat_.dim(), 0), 1) * stoch_ind_get(ARG2_vec_, (i) - 1);
                                                                  ^~~~~~~~~~~~~
P_1_code_MID_1_nfCode.cpp:505:13: note: in instantiation of function template specialization 'rcFun_R_GlobalEnv5<CppAD::AD<double>>' requested here
Interm_24 = rcFun_R_GlobalEnv5< CppAD::AD<double> >(Interm_22, Interm_23);
            ^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/nimble/include/nimble/nimDerivs_atomic_dyn_ind.h:265:22: note: candidate function template not viable: no known conversion from 'const int *' to 'NimArr<1, CppAD::AD<double>> &' for 1st argument
stoch_ind_set1_c<I_> stoch_ind_set(NimArr<1, CppAD::AD<double> > &x, const I_ &i) {
                     ^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/nimble/include/nimble/nimDerivs_atomic_dyn_ind.h:342:1: note: candidate function template not viable: requires 3 arguments, but 2 were provided
stoch_ind_set(NimArr<2, CppAD::AD<double> > &x,
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/nimble/include/nimble/nimDerivs_atomic_dyn_ind.h:431:1: note: candidate function template not viable: requires 4 arguments, but 2 were provided
stoch_ind_set(NimArr<3, CppAD::AD<double> > &x,
^
4 warnings and 1 error generated.
make: *** [P_1_code_MID_1_nfCode.o] Error 1
paciorek commented 8 months ago

@perrydv I'm not seeing what is going wrong here, but I'll note it seems to relate to lack of indexing. If I change to

out[1:p , i] <- mat[ 1:p, i] * vec[i]

it compiles.

I think this needs attention from you ASAP in case there is implication for release, which it feels like there will be.

perrydv commented 8 months ago

Thanks for the report. I just made a PR to try to fix this.