Closed paciorek closed 10 months ago
I ran the disease mapping example from our examples with HMC and results compared favorably (in terms of mixing) and were numerically similar to our default MCMC.
Not sure if we want an actual nimbleHMC test that assesses validity of the HMC results for a CAR model.
@paciorek I'm very happy to see this. Thank you.
I'm pushing some (very) minor changes to style/comments.
One minor question, should the line (# 293 of CAR.R
):
M <- rep(-1234, length = N)
instead be:
M <- rep(-1234L, length = N)
? It seems that without this, some elements of M (corresponding to zero neighbor islands) might persist as double value -1234, and have derivatives tracked. But perhaps this wouldn't cause any downstream problems, I'm not sure.
This is minor issue @paciorek but I noticed that in RCfunction_core.R, in fxnsNotAllowedInAD
, the dcar_normal
and dcar_proper
entries do not appear to be correctly entered. They have an extra "d". They were correctly entered in distsNotAllowedInAD
. This might mean that they were error-trapped in one place but not another. I am not tracking it down at the moment. And it might become moot from this branch.
Ok, I've fixed the issue Perry raised. I think it's ok in terms of M
. We can't change M
to be integer as it contains doubles and setting M
in ignore
causes a type mismatch.
@perrydv given you responded does means you looked over the changes and things seem fine to you such that I can merge this in?
Some tests were failing because of typo in numNumeric
.
@paciorek Thanks for this. My previous comment was actually something I just noticed while working on something else, so no I had not done a thorough code review. Now I've looked through it and here are my notes:
In CAR_calcNumIslands
nimInteger(..., init=FALSE)
if you are going to populate values in a for loop anyway. Someone somewhere will celebrate those microseconds ;-)nimInteger(N)
and rep(0L, N)
a few lines apart. They should be equivalent. It doesn't matter. If you want a more consistent style you could go with the first (with init=FALSE), or either.indToVisit <- nimInteger(nNeighbors, value = adj[adjStartInd:(adjStartInd+nNeighbors-1)])
, unless the AD processing gets confused about excluding all of that line from anything to do with AD.In CAR_calcCmatrix
init=FALSE
in nimMatrix
, although for Cmatrix
you might defensively really want the zeros because not every value will be filled by the following for loops.In the C++
Type(R_NaN)
we are using CppAD::numeric_limits<Type>::quiet_NaN()
if(Type(count) != Type(L))
should be if(count != L)
. It would take different syntax to have the conditional recorded on the AD tape, but if I follow these are fixed values that can be baked into the tape (and in this case would bake in a NaN and it would all be disaster anyway).The failure seems to be unrelated to this PR. I will look into it on devel.
This sets up templated versions of
dcar_normal
anddcar_proper
for AD and enablesbuildDerivs
for various auxiliary nimbleFunctions functions inCAR.R
. In the latter case, this requires various uses ofADbreak
,ignore
and setting up index variables as integers.Testing now includes tests that AD-enabled CAR models will build and basic (non-HMC) MCMC will run. Substantive HMC testing could be done in
nimbleHMC
.@perrydv a glance over the AD-enabled code would be helpful.
@danielturek just an FYI; I'm sure you'll be happy to see this.