search

nimble-dev / nimble

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

untrapped case of stochastic indexing producing a vector #1500

Closed paciorek closed 1 week ago

paciorek commented 1 month ago

Simplifying from a user post on 2024-10-22, this causes model building to fail because of AD. Note that in testing we have a latent dcat with a dependent dnorm and that is fine. So something about the two dcats may be the issue.

code <- nimbleCode({
    y ~ dcat(p[z,1:2])
    for(i in 1:2) {
      p[i,1]~dunif(0,1)
      p[i,2]<-1-p[i,1]
    }
    z ~ dcat(p[1,1:2])

})
m <- nimbleModel(code,buildDerivs=TRUE)
cm <- compileNimble(m)
using C++ compiler: ‘g++ (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
In file included from P_1_code_MID_1.h:6,
                 from P_1_code_MID_1.cpp:7:
/accounts/vis/paciorek/R/x86_64-pc-linux-gnu-library-ubuntu-22.04/4.4/nimble/include/nimble/NimArr.h: In instantiation of ‘NimArr<nDim, T>& VecNimArr<ndim, T>::operator[](unsigned int) [with int ndim = 1; T = double]’:
P_1_code_MID_1.cpp:28:10:   required from here
/accounts/vis/paciorek/R/x86_64-pc-linux-gnu-library-ubuntu-22.04/4.4/nimble/include/nimble/NimArr.h:1862:35: warning: format ‘%i’ expects argument of type ‘int’, but argument 3 has type ‘std::vector<NimArr<1, double>, std::allocator<NimArr<1, double> > >::size_type’ {aka ‘long unsigned int’} [-Wformat=]
 1862 |           "values.size() is only %i\n",
      |                                  ~^
      |                                   |
      |                                   int
      |                                  %li
 1863 |           i, values.size());
      |              ~~~~~~~~~~~~~         
      |                         |
      |                         std::vector<NimArr<1, double>, std::allocator<NimArr<1, double> > >::size_type {aka long unsigned int}
/accounts/vis/paciorek/R/x86_64-pc-linux-gnu-library-ubuntu-22.04/4.4/nimble/include/nimble/NimArr.h: In instantiation of ‘NimArr<nDim, T>& VecNimArr<ndim, T>::operator[](unsigned int) [with int ndim = 2; T = double]’:
P_1_code_MID_1.cpp:29:10:   required from here
/accounts/vis/paciorek/R/x86_64-pc-linux-gnu-library-ubuntu-22.04/4.4/nimble/include/nimble/NimArr.h:1862:35: warning: format ‘%i’ expects argument of type ‘int’, but argument 3 has type ‘std::vector<NimArr<2, double>, std::allocator<NimArr<2, double> > >::size_type’ {aka ‘long unsigned int’} [-Wformat=]
 1862 |           "values.size() is only %i\n",
      |                                  ~^
      |                                   |
      |                                   int
      |                                  %li
 1863 |           i, values.size());
      |              ~~~~~~~~~~~~~         
      |                         |
      |                         std::vector<NimArr<2, double>, std::allocator<NimArr<2, double> > >::size_type {aka long unsigned int}
/accounts/vis/paciorek/R/x86_64-pc-linux-gnu-library-ubuntu-22.04/4.4/nimble/include/nimble/NimArr.h: In instantiation of ‘NimArr<nDim, T>& VecNimArr<ndim, T>::operator[](unsigned int) [with int ndim = 1; T = CppAD::AD<double>]’:
P_1_code_MID_1.cpp:113:10:   required from here
/accounts/vis/paciorek/R/x86_64-pc-linux-gnu-library-ubuntu-22.04/4.4/nimble/include/nimble/NimArr.h:1862:35: warning: format ‘%i’ expects argument of type ‘int’, but argument 3 has type ‘std::vector<NimArr<1, CppAD::AD<double> >, std::allocator<NimArr<1, CppAD::AD<double> > > >::size_type’ {aka ‘long unsigned int’} [-Wformat=]
 1862 |           "values.size() is only %i\n",
      |                                  ~^
      |                                   |
      |                                   int
      |                                  %li
 1863 |           i, values.size());
      |              ~~~~~~~~~~~~~         
      |                         |
      |                         std::vector<NimArr<1, CppAD::AD<double> >, std::allocator<NimArr<1, CppAD::AD<double> > > >::size_type {aka long unsigned int}
/accounts/vis/paciorek/R/x86_64-pc-linux-gnu-library-ubuntu-22.04/4.4/nimble/include/nimble/NimArr.h: In instantiation of ‘NimArr<nDim, T>& VecNimArr<ndim, T>::operator[](unsigned int) [with int ndim = 2; T = CppAD::AD<double>]’:
P_1_code_MID_1.cpp:114:10:   required from here
/accounts/vis/paciorek/R/x86_64-pc-linux-gnu-library-ubuntu-22.04/4.4/nimble/include/nimble/NimArr.h:1862:35: warning: format ‘%i’ expects argument of type ‘int’, but argument 3 has type ‘std::vector<NimArr<2, CppAD::AD<double> >, std::allocator<NimArr<2, CppAD::AD<double> > > >::size_type’ {aka ‘long unsigned int’} [-Wformat=]
 1862 |           "values.size() is only %i\n",
      |                                  ~^
      |                                   |
      |                                   int
      |                                  %li
 1863 |           i, values.size());
      |              ~~~~~~~~~~~~~         
      |                         |
      |                         std::vector<NimArr<2, CppAD::AD<double> >, std::allocator<NimArr<2, CppAD::AD<double> > > >::size_type {aka long unsigned int}
P_1_code_MID_1_nfCode.cpp: In member function ‘virtual CppAD::AD<double> y_L1_UID_4::calculate_ADproxyModel(const indexedNodeInfo&) const’:
P_1_code_MID_1_nfCode.cpp:120:37: error: invalid ‘static_cast’ from type ‘CppAD::AD<double>’ to type ‘int’
  120 | Interm_8.setMap((**ADproxyModel_p), static_cast<int>(((**ADproxyModel_z)[0] - 1) * (**ADproxyModel_p).strides()[0]), (**ADproxyModel_p).strides()[1], 2);
      |                                     ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
perrydv commented 3 weeks ago

I believe this is simply not supported and should be error-trapped and noted in documentation.

The issue is not the nestedness of dcats but rather the non-scalar stochastic indexing of p[z, 1:2]. Note the following (adapted from the last example in test-ADmodels) also does not work:

set.seed(1)
code <- nimbleCode({
    for(i in 1:5) {
        y[i, 1:2] ~ dmnorm(mu[z[i], 1:2], cov = cov[1:2, 1:2])
        z[i] ~ dcat(p[1:3])
    }
    for(i in 1:3)
      for(j in 1:2)
        mu[i,j] ~ dnorm(0,1)
    p[1:3] ~ ddirch(alpha[1:3])
    cov[1,1:2] <- c(1, 0)
    cov[2, 1:2] <- c(0,1)
})

model <- nimbleModel(code, 
  data = list(y = matrix(rnorm(10), ncol=2)),
  inits = list(z = c(1,2,2,1,3), 
  mu = matrix(rnorm(6), ncol=2), 
  alpha = 1:3, p = c(.2,.3,.5)),
                     buildDerivs=TRUE)

The problem is non-scalar stochastic indexing: mu[z[i], 1:2].

Making AD work through stochastic indexing is non-trivial, so we currently only have it working when the result is a scalar.

paciorek commented 3 weeks ago

So that suggests we could suggest a work-around along these lines:

    ptmp[1:2] <- c(p[z,1],p[z,2])

@perrydv is that correct?

perrydv commented 3 weeks ago

Hmm, I guess so. Typically there will be another index (through time) and more categories than two, so it could get a bit cumbersome and potentially inefficient but that does look like a work-around for now. Since this looks like a discrete Markov model, another potential way forward is the HMM distributions in nimbleEcology.