RTMB advector issue when assigning to Matrix object

[Sebastián Pardo]

Hi everyone,

I'm trying to implement a geostatistical model using RTMB, so in other words, I'm rewriting sdmTMB models in RTMB. I'm encountering an issue when incorporating anisotropy: creating the precision matrix Q needed for anisotropy causes issues with advectors, in that there are no methods yet for assigning a RHS advector value into a LHS sparse matrix via indexing:

> mat <- Matrix::Matrix(0, nrow = 2 , ncol = 2) > mat[1,1] <- advector(2) Error in .local(x, i, j, ..., value) : not-yet-implemented 'Matrix[<-' method

To highlight this issue in relation to geostatistical precision matrices, I've written a short-ish MWE that helps with understanding the barrier with implementing anisotropy in RTMB:

https://gist.github.com/sebpardo/9a7c63f2a9fd57b5e664c07dd14e7ec8

It seems to me there are a few  ways of potentially resolving this issue:

1. Implement 'Matrix[<-' method for advectors.

2. Compile R_inla.hpp as R objects and use those Q SPDE functions instead of R ones.

3. Move away from using sparse matrices for the precision matrix to a simpler "base" matrix (`advector`). I explored using this but encountered other errors along the way, and I'm not sure if it would be a fruitful avenue.

4. Calling `RTMB::getValues` within `Q_spde_aniso()` to remove RHS advector and avoid matrix assignment issue, for example replacing line 72 of the MWE with:

G1_aniso[TV[i, 2] , TV[i, 1] ] <- G1_aniso[TV[i, 2] , TV[i, 1] ] + RTMB:::getValues(Gtmp[i, 1, 2])
and so on until line 80. However, this is likely to mess up with the tape and is probably not advisable, and also throws an error ("Not compatible with requested type: [type=language; target=double]") when `Q_spde_aniso()`returns the final object.

I'm unsure how to proceed from here, so any thoughts and suggestions are most welcome!

Thanks,

Seb

[Kasper Kristensen]

The recently added RTMB::AD constructor can sort out such overload issues - example very close to what you want is on the help page ?AD.

In general, if you need to sub-assign to non-AD object from AD object, it's probably best to use e.g. :

G1_aniso <- AD(Matrix::Matrix(0, nrow = n_s, ncol = n_s))
Gtmp <- AD(array(0, dim = c(n_tri, 3, 3)))

In non-AD contexts the AD(.) is ignored, while in AD contexts an appropriate conversion takes place.

[Sebastián Pardo]

Thank you very much for this feature and for your reply, Kasper! Wrapping those objects with AD() solved the matrix assignment issue. I'm now encountering the following error at the end of Q_spde_aniso() when the precision matrix is being returned:

# Error: Not compatible with requested type: [type=language; target=double].
# > traceback()
# 7: stop(structure(list(message = "Not compatible with requested type: [type=language; target=double].", 
#        call = NULL, cppstack = NULL), class = c("Rcpp::not_compatible", 
#    "C++Error", "error", "condition")))
# 6: advec(x)
# 5: advector(y)
# 4: .MakeTape(mapfunc, obj$env$par)
# 3: MakeADFunObject(data, parameters, reportenv, ADreport = ADreport, 
#        DLL = DLL)
# 2: obj$retape()
# 1: MakeADFun(nll, par, random = c("rf1"), silent = FALSE)

I encountered this same error when wrapping matrixes with RTMB::getValues(), and have tried troubleshooting it,  albeit unsuccessfully. Any suggestions on how to address this?

Thanks again for all your help,

Seb

--
To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/tmb-users/15bf9487-8f75-4aa1-a0c2-0436839571f2n%40googlegroups.com.

[Kasper Kristensen]

rf1 ~ dgmrf(0, Q = tau_U^2 * Qa)

is a formula, think you meant this:

rf1 %~% dgmrf(0, Q = tau_U^2 * Qa)