use of `inprod` with 1-length vectors fails

[paciorek]

We probably want this to work when n=1 so that a user's code could be agnostic to the value of n.

code <- nimbleCode({
x <- inprod(ug[1:n],   utraits[1:n])
})
m = nimbleModel(code, constants = list(n=1))
cm = compileNimble(m)

> printErrors()
P_29_code_MID_32_nfCode.cpp:18:64: error: invalid operands to binary expression ('CwiseNullaryOp<seqClass<Eigen::Matrix<double, -1, -1, 0>, int, int, int, useBy>, Eigen::Matrix<double, -1, -1, 0>>' and 'int')
(**model_x)[0] = eigenInprod((**model_ug)[(nimSeqByD(1,1,1,0)) - 1], (**model_utraits)[(nimSeqByD(1,1,1,0)) - 1]);
                                          ~~~~~~~~~~~~~~~~~~~~ ^ ~
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/v1/iterator:855:1: note: candidate template ignored: could not match 'reverse_iterator' against 'CwiseNullaryOp'
operator-(const reverse_iterator<_Iter1>& __x, const reverse_iterator<_Iter2>& __y)
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/v1/iterator:1311:1: note: candidate template ignored: could not match 'move_iterator' against 'CwiseNullaryOp'
operator-(const move_iterator<_Iter1>& __x, const move_iterator<_Iter2>& __y)
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/v1/iterator:1719:1: note: candidate template ignored: could not match '__wrap_iter' against 'CwiseNullaryOp'
operator-(const __wrap_iter<_Iter1>& __x, const __wrap_iter<_Iter2>& __y) _NOEXCEPT
^

[paciorek]

Use of inprod(x[1:n],y[1:n]) and setting n=1 is fine in a nimbleFunction.

The problem is in model code where n=1 causes use of nimSeqByD instead of the Eigen .block method. Need to understand the processing steps better...

[paciorek]

More specifically this happens because of the dropThisDim code in sizeIndexingBracket. At that point, we don't directly know that this is inside inprod, so not clear how to handle the inprod case differently from other cases where we drop (I'm not seeing at the moment what those other cases would be or why we want to drop in those cases).

I guess we would have to set some flag in sizeBinaryReduction that gets passed deeply into recurseSetSizes and exprClasses_setSizes so that sizeIndexingBracket can receive an argument that tells it not to drop.

Or (this feels hacky), modify 1:1 perhaps as 1:1.001 or the like to indicate that the dimension shouldn't be dropped and then change to 1:1 after the dimension dropping code. Or set the type to integer if dimension shouldn't be dropped. That feels fragile. Or something more explicit like 1:KEEP. One might have additional syntax wrapped around the args to inprod, like inprod(foo(x[1:3]),y[1:7]), so I think one would only modify 1:1 in arrays going directly into inprod. However, what about inprod(x[1:1,1:3],y[1:3]). Perhaps that would rightfully error out, just as any inner product of a matrix and vector would.

And of course one could have inprod(x[2:n],y[2:n]) or inprod(x[n:5],y[n:5]) and still face the same question.

[paciorek]

Ok, I'm inclined to do the following to handle this as a special case when we are in sizeBinaryReduction and we detect use of inprod directly:

• look for index positions in indexing used directly on variables passed to inprod where the begin and end indices are the same. Replace the end index with KEEP.
• then in sizeIndexingBracket avoid dropping the dimension and then replace KEEP with the first index.

That said, this would still leave other similar cases unhandled such as x[1:3,1:1] %*% y[1:1].

Would like some feedback from @danielturek @perrydv

[danielturek]

@paciorek I wish I had more useful feedback for you. The proposed solution would handle common cases of the inprod problem, which maybe is sufficient for the time being.

I'm not thinking of an (easy) solution that's much more general. One could (perhaps) be to introduce a drop = FALSE option into DSL array indexing, and use this in sizeIndexingBracket to control if dimensions are reduced... but this probably has other repercussions (and implementation challenges) that I'm not beginning to imagine.