Using functions with Nimblelists in model code

[philipcindamix]

Hi All,

I have a query about using nimbleList data structures as outputs for nimbleFunctions (which are included in model code). Based on the examples in Chapter 14 I wrote the following toy nimbleFunction:

mynf <- nimbleFunction(
  setup = function(){
    my_List_Def <- nimbleList(ouput1 = double(1), ouput2 = double(2))
  },
  run = function(input1 = double(1), input2 = double(2)){
   

    output_List <- my_List_Def$new()
    output_List$ouput1 <- input1
    output_List$ouput2 <- input2
   
    returnType(my_List_Def())
    return(output_List)
  }
)
Presumably because I have specified the "setup" argument I need to assign it to another object and then supply the arguments using $run(...) to get the required output?

mynf2 <- mynf()
mynf2$run(c(1,2), diag(2))

I'm not sure if this is the correct way to do it and furthermore how do I use this function in the model code?

Can I specify a nimbleList within the model code?

Many Thanks,

Philip  

[Perry de Valpine]

Dear Philip,

It is not currently possible to return a nimbleList in model code, or to use a nimbleFunction with setup code in model code.

One solution could be to use numeric objects (vectors, matrices, arrays) as data structures following a protocol that you use consistently.  E.g. if you need to keep input1 and input2 packed into a data structure together, you could make it a matrix with two columns.  If the sizes of relevant values need to change over the course of using the model, you would need the matrix that holds them to take the maximal possible size.  If you have other user-defined functions that need to use them and determine what the lengths of relevant values are, you could either hold the lengths in another model variable or you could encode the end of the relevant values in each column with -1 or some such value.  I hope that makes some sense.

If those ideas don't get you moving in a useful direction, feel free to explain more of what you need.

HTH

Perry

--
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/d11c4686-fcf5-4c39-b359-cf08992c3572n%40googlegroups.com.

[philipcindamix]

Great, Thanks for the prompt response.

Happily, the two outputs have similar dimensions so I can stick them together.

Thanks again,

Philip 

[philipcindamix]

Hi All, 

Unfortunately, I have encountered another problem. 

The objective of the function I mentioned in the previous post is to return the conditional mean vector and covariance matrix (stuck together as one matrix) based on the original mean vector, covariance matrix and conditional value(s) of mu (taking care to use inverse rather than solve):

conditional_MNV <- nimbleFunction(

   run = function(mu_vec = double(1), cov_mat = double(2), mu_vec_cond = double(1)){
     
     mat_size <- length(mu_vec)
     cond_vec_size <- length(mu_vec_cond)
     index_end <- mat_size-cond_vec_size
     
     #convert vectors to matrices
     mu_vec_mat <- matrix(0, nrow = mat_size, ncol = 1)
     mu_vec_mat[,1] <- mu_vec
     mu_vec_cond_mat <- matrix(0, nrow = cond_vec_size, ncol = 1)
     mu_vec_cond_mat[,1] <- mu_vec_cond
     
     resMatrix <- matrix(0, nrow = index_end, ncol = index_end+1)
     resMatrix[1:index_end,1:index_end] <- cov_mat[1:index_end,1:index_end] - cov_mat[1:index_end,(index_end+1):mat_size]%*%inverse(cov_mat[(index_end+1):mat_size,(index_end+1):mat_size])%*%cov_mat[(index_end+1):mat_size,1:index_end]
     resMatrix[1:index_end,index_end +1] <- mu_vec_mat[1:index_end,] + cov_mat[1:index_end,(index_end+1):mat_size]%*%inverse(cov_mat[(index_end+1):mat_size,(index_end+1):mat_size])%*%(mu_vec_cond_mat -      mu_vec_mat[(index_end+1):mat_size,])
   
    returnType(double(2))
    return(resMatrix)
  }
)

Defining the following values to test the function:

R = matrix(c(1, 0.7, 0.2,
             0.7, 1, -0.5,
             0.2, -0.5, 1), 3)

sd_vec <- c(.2,0.5,0.1)
cov_mat <-diag(sd_marg)%*%R%*%diag(sd_marg)
mu_vec <- c(5, 2,1)
mu_condit_val <- 3

However, when I want to include this in the nimbleCode I get a large number of compiler errors. 

mod_nimble2<- nimbleCode({
  condit_dist_1[1:2,1:3] <- conditional_MNV(mu_vec[1:3], cov_mat[1:3, 1:3], mu_condit_val[1:1])
})

Any help would be greatly appreciated.

Philip 

[Perry de Valpine]

Hi Philip,

You can test the function and its compilation outside of a model.  I am guessing sd_diag was supposed to be sd_vec, so I've made that change.

sd_vec <- c(.2,0.5,0.1)
cov_mat <-diag(sd_vec)%*%R%*%diag(sd_vec)

mu_vec <- c(5, 2,1)
mu_condit_val <- 3

conditional_MNV(mu_vec, cov_mat, mu_condit_val) # uncompiled: works

C_conditional_MNV <- compileNimble(conditional_MNV) # compile it
C_conditional_MNV(mu_vec, cov_mat, mu_condit_val) # compiled: works

So the problem is how you are using it in a model.

Unfortunately a friction point in nimble is that compilation supports true scalars, but R does not have true scalars (only length 1 vectors).  As a result, the single value mu_condit_val[1:1] is not passed as a vector successfully.  If you will always have a scalar here, declare it as double(0) instead of double(1).  Alternatively, if sometimes you. need a vector, you can pad it with an extra 0, which you can then not use within the conditional_MVN function.

Also I would suggest a different strategy within conditional_MVN if computational efficiency is a concern.  Instead of repeating the step

cov_mat[1:index_end,(index_end+1):mat_size]%*%inverse(cov_mat[(index_end+1):mat_size,(index_end+1):mat_size])

you could do that once and re-use the result in both places where it's needed.

Another strategy would be to do this calculation within the model (or as a separate user-defined function) and then have separate functions to obtain the conditional mean and conditional covariance, using this shared part of the calculation as input.

HTH

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/e2bc249d-35a0-4d31-97e3-92167836c968n%40googlegroups.com.

[philipcindamix]

Many thanks for explaining the issue and pointing out the computational inefficiencies.

I modified the code and it did work, however, I ended up having to pad the matrix object in for both rows and columns. It was giving me an error when I returned a 1x2 matrix for example. 

conditional_MNV <- nimbleFunction(
 
  run = function(mu_vec = double(1), cov_mat = double(2), mu_vec_cond_pad = double(1)){
   
   
    mu_vec_cond <- numeric(length(mu_vec_cond_pad)-1, value = -999)
    mu_vec_cond[1:(length(mu_vec_cond_pad)-1)] <- mu_vec_cond_pad[1:(length(mu_vec_cond_pad)-1)]

    mat_size <- length(mu_vec)
    cond_vec_size <- length(mu_vec_cond)
    index_end <- mat_size-cond_vec_size
   
    #convert vectors to matrices
    mu_vec_mat <- matrix(0, nrow = mat_size, ncol = 1)
    mu_vec_mat[,1] <- mu_vec
    mu_vec_cond_mat <- matrix(0, nrow = cond_vec_size, ncol = 1)
    mu_vec_cond_mat[,1] <- mu_vec_cond
   

    #Inverse VCOV
    inv_cov_mat <- matrix(0, nrow = cond_vec_size, ncol = cond_vec_size)
    inv_cov_mat[1:cond_vec_size,1:cond_vec_size] <- inverse(cov_mat[(index_end+1):mat_size,(index_end+1):mat_size])
   
    resMatrix <- matrix(0, nrow = index_end+1, ncol = index_end+1) #padding row too
    resMatrix[1:index_end,1:index_end] <- cov_mat[1:index_end,1:index_end] - cov_mat[1:index_end,(index_end+1):mat_size]%*%inv_cov_mat[1:cond_vec_size,1:cond_vec_size]%*%cov_mat[(index_end+1):mat_size,1:index_end]
    resMatrix[1:index_end,index_end +1] <- mu_vec_mat[1:index_end,] + cov_mat[1:index_end,(index_end+1):mat_size]%*%inv_cov_mat[1:cond_vec_size,1:cond_vec_size]%*%(mu_vec_cond_mat-mu_vec_mat[(index_end+1):mat_size,])
   
    returnType(double(2))
    return(resMatrix)
  }
)

Thanks again,

Philip 

[Perry de Valpine]

Thanks.  I'm guessing this issue is due to code on the nimbleModel side.  In any case, you've made it work.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/41a3ff8c-f0db-44d1-ba4a-435ca5f55a78n%40googlegroups.com.