How to support `model$getParam(nodes[i], 'mean')`

[perrydv]

Here is a question from starting work on model$getParam. In nimble a while back we encountered the need to use getParam for any element from a vector of nodes, so we built support for model$getParam(nodes[i], 'mean'), where nodes is a setup output and nodes[i] at run-time uses the i-th node. This allowed an algorithm to iterate over a vector of nodes, using the same parameter for each. IIRC in compilation we essentially use model$getParam(nodes, i, 'mean'), i.e. we rearrange so that the index is an argument we manage, and we create a nodeFxnVector for all of nodes, which allows in C++ to pick out the i-th one dynamically. Also note that for this purpose, nodes could be ordered arbitrarily and there is no need for topological sorting, and moreover the ordering of nodes shouldn't be changed in any way the user doesn't realize, because they may have logic to how it is ordered in relation to how they are ordering related data structures. This all suggests we need to support arbitrary ordering even when it is somewhat unnatural with the new node abstractions. I would imagine often a user's (or algorithm's) ordering will be "natural" and that will match our compact representations (e.g. x[1:10000]), but we need to support anything. Does this make sense @paciorek?

[paciorek]

I'm feeling rather fuzzy here, but a few offhand comments that may not be helpful:

• we use getParam with indexing of a node vector in Laplace (via normTooling.R in nimble) when we want to get the known analytic gradient rather than using AD
• nodeRanges don't have any particular ordering. I think ordering only comes in when we get to the calcRule/calcRange stage of things.
• do we need something like a simplified instrList without the ordering?

[perrydv]

Had a good discussion with @paciorek about this today. Here is some further thinking:

• We agree for backward compatibility we need to support "nodes[i]" where nodes is not the old-style character vector of node names but may be a nodeRange (list of node Ranges?) or nList(instr_nClass) of compact representations. This means we need to maintain an ordering that we can somewhat easily access for an overall index.
• The plot seems to thicken if we are to handle all of the indexing machinery. If nodes is a nodeRange from a declaration with 3 indices and non-simple arrangement (i.e., one of the kinds of iteration cases under way for calculate), and say it represents a large number of nodes, then it may be non-trivial to find an arbitrary iteration within it. For example, the 50,000th node in nodes might be y[561, 253, 2] by some complicated arrangement. Getting to that from nodes[50000] might not be trivial.
• We do need to handle this for backward compatibility, even if we accept that it may be inefficient. But then for cases where one wants to iterate through param values of a large set of values and do it in a newer way, what can we provide as a good way to do this?
• One idea would be to allow getParam to return a container of all the param values, either an nList or a padded dimension of whatever numeric object would otherwise be returned (a matrix for a set of vectors, etc). A problem is that we don't really want to force allocation and population of huge memory objects if instead one could iterate through one at a time.
• Another idea would be to expose iteration or iterators in the DSL in some way. So that a DSL programmer could say they want to step through getParam results node by node from a nodeRange and be able to explicitly call for the next one. This would take care for supporting parallelization too. Right now this idea seems a bit complicated, but it would allow efficient programming.
• We don't really want to create new features right now, and iterating node by node is not very widespread, but it does happen I think in BNP and also conjugate samplers, both important cases.
• This is all relevant for model$calculate(nodes[i]) and similar methods.

[paciorek]

I am very likely missing something obvious, but regarding how to get to the 50000th node:

If we have a single instr_nClass object, when we do calculate, we iterate over the length of the instr (len) and can readily get the idx value we need at each iteration: idx = values[1]+i or idx = values[i], which is then passed into calc_one.

So for getParam, let's assume we will also always be working with an nList(instr_nClass) (but not necessarily in sortID order) since we have those in compiled code and don't have nodeRanges. So say we have getParam(depNodes[i], 'mean') and this becomes getParam(one_instr, i, 'mean') in compiled code. Can't we pass the value of i into something akin to calc_1_seq, calc_1_mat, etc. but where instead of iterating over all the elements in the instr, we are just picking out the ith one via simple calculations with the values of the instr? I.e., we still have idx = one_instr->values[1]+i or idx=one_instr->values[i]. In multi-index cases, it's a bit more complicated as we have to work out idx[0],idx[1], etc., based on determining the index value for each of the index slots, but we know the lens of each indexRange, so it's simple arithmetic. So {561, 253, 2} comes directly from the machinery I have just built.

(I am glossing over that depNodes might involve multiple instrs. But we know the number of elements of each, so again it seems straightforward.)

Granted for the multi-index case there is overhead from the fact that the for loop is outside of all this and for each iteration we have to do the lookup arithmetic, unlike standard calculate where the for loop is "internal", but the lookup is just a bit of arithmetic. (And perhaps there is a way to cache something useful.)

What am I missing?

[perrydv]

Ok I see what you're saying. That part of the problem is not hard. I think this "random access" need to index anywhere into a set of instructions might still introduce some changes or extensions to what we've got. Let's discuss in person.