nimble function - issues with indexing

[Rowenna]

Hi all,

I've written a function to sample covariate values based on the nearest raster grid cell to a location. I'm getting an error related to indexing but can't figure it out. 

Any suggestions would be much appreciated!

Test data is attached:

y - x,y locations

S - x,y values from raster grid 

cov1-3 covariates from raster grids

#covariate function

covVals2 <- nimbleFunction(

run = function(y = double(2), 

S = double(2),

cov1= double(1),

cov2 = double(1),

cov3 = double(1)) {

dist <- numeric(nrowNim(S))

for(i in 1:nrowNim(S)) {

dist[i] <- sqrt(((y[,1] - S[i, 1]) ^2) + ((y[,2] - S[i, 2]) ^2)) #calculate nearest 

id <- which(dist == min(dist)) 

x1.val <- cov1[id] 

x2.val <- cov2[id]

x3.val <- cov3[id] }

cov.vals <- array(c(x1.val, x2.val, x3.val), dim = c(3, nrowNim(y)))

return(cov.vals)

returnType(double(2)) 

}

) 

#test in loop

bla <- data.frame(array(0, c(3, 6)))

for (t in 1:nrowNim(test2mu)) {

bla[t] <- covVals2(test2mu[t], testS, xTest1, xTest2, xTest3)

}

Error in y[, 1] : incorrect number of dimensions

Cheers,

Rowenna

[Daniel Eacker]

Hi Rowenna,

It appears that your distance calculation is producing a vector at each loop in your for loop, but you only have a scalar to hold the output. I didn't use your example data, but I created a simple example and your function seems to work fine:

require(nimble)
require(raster)

#Test data is attached:
#  y - x,y locations
#S - x,y values from raster grid
#cov1-3 covariates from raster grids

x1 = runif(100, 0, 1000)
y1 = runif(100, 0, 1000)
pts = SpatialPoints(cbind(x1, y1))

x2 = seq(0, 1000, 10)
y2 = seq(0, 1000, 10)
rast = raster(xmn = 0, xmx = 1000, ymn = 0, ymx = 1000, res = 100)
rast[] = rpois(ncell(rast), 20)
rast1 = rast
rast1[] = rpois(ncell(rast), 20)
rast2 = rast
rast2[] = rpois(ncell(rast), 20)

#covariate function
covVals2 <- nimbleFunction(
  run = function(y = double(2),
                 S = double(2),
                 cov1= double(1),
                 cov2 = double(1),
                 cov3 = double(1)) {

    dist <- matrix(0, nrow = length(S[,1]), ncol=length(y[,1]))
    for(i in 1:length(S[,1])) {
      dist[i,] <- sqrt(((y[,1] - S[i, 1]) ^2) + ((y[,2] - S[i, 2]) ^2)) #calculate nearest
    }
    # now loop over locations to get covariate ids
    id = numeric(length(y[,1]))
    for(i in 1:length(S[,1])) {
    id[i] <- which(dist[i,] == min(dist[i,]))

     }
    
      x1.val <- cov1[id]
      x2.val <- cov2[id]
      x3.val <- cov3[id]

    cov.vals <- array(c(x1.val, x2.val, x3.val), dim = c(3, length(y[,1])))
    return(cov.vals)
    returnType(double(2))
  }
)

covVals2(y = coordinates(pts), S = coordinates(rast), cov1 = rast@data@values, cov2 = rast1@data@values,cov3 = rast2@data@values)

Please let me know if this helps.

Thanks,

Dan

[Daniel Eacker]

Sorry had a bug in my code. Should be correct now:

require(nimble)
require(raster)

#Test data is attached:
#  y - x,y locations
#S - x,y values from raster grid
#cov1-3 covariates from raster grids

x1 = runif(200 0, 1000)
y1 = runif(200, 0, 1000)
pts = SpatialPoints(cbind(x1, y1))

rast = raster(xmn = 0, xmx = 1000, ymn = 0, ymx = 1000, res = 100)
rast[] = rpois(ncell(rast), 20)
rast1 = rast
rast1[] = rpois(ncell(rast), 20)
rast2 = rast
rast2[] = rpois(ncell(rast), 20)

#covariate function
covVals2 <- nimbleFunction(
  run = function(y = double(2),
                 S = double(2),
                 cov1= double(1),
                 cov2 = double(1),
                 cov3 = double(1)) {
    dist <- matrix(0, nrow = length(S[,1]), ncol=length(y[,1]))
    for(i in 1:length(S[,1])) {
      dist[i,] <- sqrt(((y[,1] - S[i, 1]) ^2) + ((y[,2] - S[i, 2]) ^2)) #calculate nearest
    }
    # now loop over locations to get covariate ids
    id = numeric(length(y[,1]))

    for(i in 1:length(y[,1])) {
    id[i] <- which(dist[,i] == min(dist[,i]))

     }
        
      x1.val <- cov1[id]
      x2.val <- cov2[id]
      x3.val <- cov3[id]
    cov.vals <- array(c(x1.val, x2.val, x3.val), dim = c(3, length(y[,1])))
    return(cov.vals)
    returnType(double(2))
  }
)

covVals2(y = coordinates(pts), S = coordinates(rast), cov1 = rast@data@values, cov2 = rast1@data@values,cov3 = rast2@data@values)

[Rowenna Gryba]

Many thanks!!!

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/c7c5d6ef-1685-4698-b253-11c98e0f52a2n%40googlegroups.com.

[Rowenna Gryba]

Hi again Daniel,

Thanks again for the above. One issue I'm still having is trying to get the function to work when looping over locations of y but outside of the function. The reason I want to do this is to be able to use the previous time step of a latent state.

So ideally I would be able to run the function in a loop like:

for (t in 2:length(y[,1])) {
bla[t] <- covVals2(y[t-1,1:2], testS, xTest1, xTest2, xTest3)
}

I just can't work around in my head how to modify it without getting the Error in y[, 1] : incorrect number of dimensions. Sorry that I didn't include this in the initial post!

Cheers,

Rowenna

[Rowenna Gryba]

Solved! 

Needed to index y with [[]] instead of [,1].

covVals2 <- nimbleFunction(
run = function(y = double(2),
S = double(2),
cov1= double(1),
cov2 = double(1),
cov3 = double(1)) {

dist <- matrix(0, nrow = length(S[,1]), ncol=length(y[[1]]))

for(i in 1:length(S[,1])) {

dist[i,] <- sqrt(((y[[1]] - S[i, 1]) ^2) + ((y[[2]] - S[i, 2]) ^2)) #calculate nearest

}
id <- which(dist == min(dist))

x1.val <- cov1[id]
x2.val <- cov2[id]
x3.val <- cov3[id]

cov.vals <- array(c(x1.val, x2.val, x3.val), dim = c(3, length(y[[1]])))
return(cov.vals)
returnType(double(2))
}
)

bla <- data.frame(array(0, c(3, 5)))
for (t in 2:length(test2mu[,1])) {
bla[t] <- covVals2(test2mu[t-1,1:2], testS, xTest1, xTest2, xTest3)
}

[Perry de Valpine]

Hi Rowenna,

You shouldn't need (or want) to use double brackets in a nimbleFunction.  Does that compile?

I think the possible confusion is that when you provide y[t-1, 1:2] as an argument for covVals2 in your model code, that creates a vector (not a matrix) because the first index is a single value.  So the argument type to covVals2 should be "y = double(1)", not "y = double(2)".  Does that work?

Also it is important to test this using compiled code to ensure the types work as intended.  I couldn't quite tell if you are testing these uncompiled or compiled.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CADKu3jodfK1UhM%3DVEw475feukTos7Pw%3D%2BR3_BO%2BtQa_X5xUwSg%40mail.gmail.com.

[Rowenna Gryba]

Hi Perry,

I had forgotten to try and compile this version! It works/compiles when I make the edits "y=double(1)" and change the index on my y variables. 

Thank you!

[Rowenna Gryba]

Apologies - I spoke too soon - I thought it was on it's way to compiling but there is an error. 

Here is my modified code:

covVals2 <- nimbleFunction(
run = function(y = double(1),

S = double(2),
cov1= double(1),
cov2 = double(1),
cov3 = double(1)){
dist <- matrix(0, nrow = length(S[,1]), ncol=length(y[,1]))
for(i in 1:length(S[,1])) {
dist[i,] <- sqrt(((y[,1] - S[i, 1]) ^2) + ((y[,2] - S[i, 2]) ^2)) #calculate nearest
}

id <- which(dist == min(dist))

x1.val <- cov1[id]
x2.val <- cov2[id]
x3.val <- cov3[id]

cov.vals <- array(c(x1.val, x2.val, x3.val), dim = c(3, length(y[,1])))
return(cov.vals)
returnType(double(2))
}
)

> CcovVals2 <- compileNimble(covVals2, showCompilerOutput = TRUE)
compiling... this may take a minute. On some systems there may be some compiler warnings that can be safely ignored.
Error: Error, wrong number of indices provided for ARG1_y_[, 1].
 This occurred for: ARG1_y_[, 1]
 This was part of the call:  ARG1_y_[, 1] - ARG2_S_[i, 1]

Thank you again for any suggestions and for bearing with me!

Rowenna

[Perry de Valpine]

It looks like something is still not right about dimensions (number of indices) of y.

I don't know if I misread your code in my previous reply.  Somewhere I thought I saw you trying to provide y as an argument to covVals2 as a vector.  But the way it is used inside of covVals2, it does look like you want it to be a matrix.  I may have been hasty.

I am a little lost at this point on exactly what you need, but maybe we can still see a solution.  Daniel's code provided above can be compiled with the following modification:

id[i] <- which(dist[,i] == min(dist[,i]))[1]

The reason this didn't compile as it was written is that `which` returns a vector, but id[i] is a scalar.  nimbleFunctions become statically typed when they compile, so types need to match, and this can be confusing in some cases.  Adding [1] at the end ensures that the type of the result will be a scalar.

(Also there is a missing comma in

x1 = runif(200, 0, 1000)

)

So now you are back to y being a matrix, which I think is what you wanted.

Thanks Daniel for putting together that code.

Also it is better to do dim(X)[1] instead of length(x[,1]) because the latter will make a copy of x[,1] just for the purpose of determining its length.

-Perry

[Rowenna Gryba]

Hi Perry - thanks for circling back to this again!

What I really want to do is use the covVals function in a nimble model with a loop over the locations. The reason I want to do this is because I'm estimating the locations (they are the latent state) and would like to use the previous location (y_t-1) to get the covariate value and use that in an equation to estimate the next location. 

e.g.

for (t in 2:length(y[,1])) {

covar[t] <- covVals2(y[t-1,1:2], testS, xTest1, xTest2, xTest3)

}

Daniel's code (thank you again) works but I believe only if you are supplying all of the y values at once as a matrix and doesn't work in a loop like above.

You did see in my initial code that did try to define y as double(2) but even when I edit it to double(1) I'm still getting the indexing error when I compile. Am I right in thinking that stochastic indexing isn't supported by nimble? 

Thank you again for your help on this!

Rowenna

[Daniel Eacker]

Hi Rowenna,

I wrote some more simulation code and slicked out the function just a bit after you described exactly what you wanted to do.

Perry, thanks for the help and catching my coding errors. I need to slow down a bit here :)

I've included an example that compiles and runs the MCMC:

Hopefully this solves your problem:

x1 = runif(200, 0, 1000)

y1 = runif(200, 0, 1000)
pts = SpatialPoints(cbind(x1, y1))

rast = raster(xmn = 0, xmx = 1000, ymn = 0, ymx = 1000, res = 100)
rast[] = rpois(ncell(rast), 20)
rast1 = rast
rast1[] = rpois(ncell(rast), 20)
rast2 = rast
rast2[] = rpois(ncell(rast), 20)

#covariate function
covVals2 <- nimbleFunction(

  run = function(y = double(1),
                 S = double(2),
                 cov1= double(1),
                 cov2 = double(1),
                 cov3 = double(1)) {

    dist <- sqrt(((y[1] - S[1:dim(S)[1], 1]) ^2) + ((y[2] - S[1:dim(S)[1], 2]) ^2)) #calculate nearest

    # now loop over locations to get covariate ids

    id = which(dist == min(dist))

    x1.val <- cov1[id]
    x2.val <- cov2[id]
    x3.val <- cov3[id]

    cov.vals <- c(x1.val, x2.val, x3.val)
    return(cov.vals)
    returnType(double(1))
  }
)

# test the function with first point
covVals2(y = coordinates(pts)[1,], S = coordinates(rast), cov1 = rast@data@values, cov2 = rast1@data@values,cov3 = rast2@data@values)

# to simulate the problem in compilied code
nimcode <- nimbleCode({
 
  for(i in 1:nLocs){
  y[i,1] ~ dunif(0, 1000)
  y[i,2] ~ dunif(0, 1000)
  }
  for(i in 2:nLocs){
  covs[1:3,i] <- covVals2(y=y[i-1,],S=rast.coords[1:nRast,1:2], cov1 = v1[1:nRast], cov2 = v2[1:nRast], cov3 = v3[1:nRast])
  }
})

nimData = list(rast.coords = coordinates(rast),v1 = rast@data@values, v2 = rast1@data@values,v3 = rast2@data@values, covs = matrix(0, nrow=3, ncol=200))
constants = list(nLocs = 200, nRast = prod(dim(rast)[1:2]))

testR <- nimbleModel(nimcode,data = nimData, constants = constants, inits = list(y = matrix(runif(constants$nLocs*2),ncol=2)))

# complete initialization of model
testR$calculate()
testR$initializeInfo()

# compile model
testC <- compileNimble(testR, showCompilerOutput = FALSE)
# MCMC sampler configurations
mcmcspec<-configureMCMC(testR, monitors=c("y","covs"))
# build the MCMC specifications
testMCMC <- buildMCMC(mcmcspec)
# complile the code
CtestMCMC <- compileNimble(testMCMC, project = testR,resetFunctions = TRUE)
# run MCMC
tic()
samples <- runMCMC(CtestMCMC, niter = 10000, nburnin=1000,thin=1,nchains=1)
toc()


colNames <- grep("covs\\[",colnames(samples))
(covs <- samples[1,colNames]) # look at the first iteration

[Daniel Eacker]

Also, the code runs fine without including an empty matrix "covs" in the data (i.e,. covs = matrix(0, nrow=3, ncol=200) ).

Hope it works out for you!

Sincerely,

Daniel

[Daniel Eacker]

Also, here is a small test to show that the function works the same as the extract() function in the raster package:

covVals2(y = coordinates(pts)[1,], S = coordinates(rast), cov1 = rast@data@values, cov2 = rast1@data@values,cov3 = rast2@data@values)

c(raster::extract(rast, pts[1,]),raster::extract(rast1, pts[1,]),raster::extract(rast2, pts[1,]))

[Rowenna Gryba]

Hi Daniel,

Thank you so much and thank you for testing in MCMC!!!!

Have a great weekend,

Rowenna

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/4e8601ab-99b5-453c-a7c8-b0011d4f7220n%40googlegroups.com.