need to handle varying blocks of indices into a variable

[paciorek]

Consider this case:

   code <- quote({
        for(i in 1:3)
            y[i, n1[i]:n2[i]] ~ dmnorm(mu[n1[i]:n2[i]], pr[n1[i]:n2[i],n1[i]:n2[i]])
    })
    model <- nimbleModel(code, constants = list(n1 = c(3,1,2), n2 = c(6,3,3)))

This is allowed in nimble so long as n1 and n2 are constants.

Currently it is error-trapped in modelDecl.R even when provided as constants.

Noting also that if I simply turn the error trap off, we end up with an error later on:

Error in names(boundExprs)[names(boundExprs) %in% c("lower", "upper")] <- paste0(names(boundExprs)[names(boundExprs) %in%  : 
  attempt to set an attribute on NULL

in reparameterizeDists.

Need to come back to this and fix up. I am expecting it to be reasonably straightforward.

[paciorek]

This is looking to be a big headache.

The basic problem is illustrated by considering the example above. We have nodes y[1,3:6], y[2,1:3], y[3,1:2]. In the nimbleModel system, we want those as three nodes, with external indexing 1:3 and internal indexing that depends on the external indexing. But that violates how we set things up with separable index sets. If we try to use current nimbleModel functionality we end up unrolling things and end up with 9 "nodes", y[1,3], y[1,4], etc. The graph processing based on that seems correct (we get correct dependencies) but then we feed the wrong info into calculation because we have 9 nodes and not 3.

I'll note that the declaration functions we generate are what we want, namely they look up the indexing needed based on the values of n1 and n2 and we can easily ensure that n1 and n2 are variables in the model.

One path forward (but not clear I should work on it now or even before first release) is to extend the notion of a matrix index range to allow for block indices, as well as extending the functionality of the indexRuleArbitraryClass to allow for this.

I think it would be feasible for the indexRangeMatrix class to have the values be a matrix of lists rather than a matrix of numbers. Then in this case we would have these elements in a 3x2 matrix, with the nodeRange containing a single three-element indexRange (and indexRangeMatrix):

1, 3:6
2, 1:3
3, 1:2

Similar things would happen in the elements of setupResults in the indexRuleArbitraryClass and we'd have to do some work on indexRuleArbitrary_applyMatrix.

Another possibility is to live with having "nodes" that are not really nodes as the result of getNodes() and then when we go to make calcRanges, fix things up so that we have the correct calculation indexing. E.g. in the example above if we calculate on all of 'y', we end up with a calcRange with this indexing, so if we take the unique values, it seems we would do the correct calculation so long as as the decl fxn idx value goes over 1,2,3:

Browse[1]>  ranges[[1]]$indexingRange$indexRanges[[1]]$values
      [,1]
 [1,]    1
 [2,]    1
 [3,]    1
 [4,]    1
 [5,]    2
 [6,]    2
 [7,]    2
 [8,]    3
 [9,]    3

and if we calculate on dependencies of mu[2] (which are y[2,1:3], y[3,1:2]), we get:

Browse[1]>  ranges[[1]]$indexingRange$indexRanges[[1]]$values
     [,1]
[1,]    2
[2,]    2
[3,]    2
[4,]    3
[5,]    3

so using unique would also get us the correct result.

If we did take this second path, we'd want to make sure that having two sets of non-constant blocks, e.g.,
y[i, n1[i]:n2[i], n3[i]:n4[i]] ~ dwish(....) also works.

And if we live with incorrect nodes, we also are accepting that default MCMC sampler assignment will not know which elements constitute a node. In that case we'd want to warn users they may want to manually assign samplers. Perhaps we can live with that.

@perrydv let's discuss before you are off, if possible.

[paciorek]

Also noting that in contradiction of my first post here, if I turn off the error trapping, the model builds and we can run calculate, but we get incorrect calculation results because we loop over 9 calls to calc_one, not 3.

[paciorek]

We decided for now to punt and just get calculation working correctly and live with "nodes" being node elements in this case.

So it turns out on 2023-06-10, I put code into originalIndexRules$apply to avoid removing duplicates. From a somewhat deep dive now, that seems to have been incorrect and motivated by wanting calculate() to return as many elements as were computed in any given deterministic calculation, e.g., y[1:2] <- foo() returning 2 elements, but we actually don't return deterministic results anyway.

So we now remove duplicate indices when applying the indexingRules, which addresses this issue.

I also need to fix a bug in indexRuleAll_setup in that it couldn't handle a case of a fixed RHS node and a varying block LHS, such as:
y[i,n1[i]:n2[i]] ~ dmnorm(mu[1:3],pr[1:3,1:3]). We need to do a bit of gymnastics to cross each value of i with the values n1[i]:n2[i] to pull out all the elements of the y node as dependencies of mu or pr.

I've also added a warning (which will probably need wordsmithing) upon model-building that node queries will return node elements and not actual nodes in such cases.