DM-14740: Stop using ndarray::EigenView indirectly in C++ code #20
Conversation
r-owen
force-pushed
the
tickets/DM-14740
branch
4 times, most recently
from
2bf77d7
to
aae840e
Compare
| @@ -52,9 +52,9 @@ double buildModelsImpl( | |||
| ndarray::Array<T,2,-2> matrix = matrixT.transpose(); | |||
| double result = 0.0; | |||
| for (ndarray::Array<double const,2>::Iterator i = parameters.begin(); i != parameters.end(); ++i) { | |||
| ellipse.setParameterVector(i->asEigen()); | |||
| ellipse.setParameterVector(ndarray::asEigenMatrix(*i)); | |||
There was a problem hiding this comment.
You could avoid the awkward *i with a foreach loop.
There was a problem hiding this comment.
I would hold off on this for now, or at least be very careful testing that change. This just reminded me that https://github.com/ndarray/ndarray/blob/master/include/ndarray/detail/NestedIterator.h#L43 might need to be changed to the InputIterator concept equivalent, and not doing so might get in the way of higher-dimensionsal ndarray iterators working properly in range-based for.
Long version of the story is that proxy iterators in C++ aren't allowed to be RandomAccessIterators (or even ForwardIterators), due to a really annoying edge case in the standard. Once upon a time, Boost.Iterator and it's not-quite-the-same-as-the-standard's iterator categories were seen as the way towards a future version of the standard without those (and other) problems and that's why I used them here. AFAIK, that idea has gone by the wayside - I haven't seen anything about fixing proxy iterators in the latest standards (though I haven't paid super close attention) - and I think we just need to stop promising that these are RandomAccessIterators, even though they can do 99% of the things a RandomAccessIterator can do.
Sorry to dump that on you right before I disappear for a week, and I'm even more sorry I didn't fix this earlier - it is a serious problem, but I keep remembering it right when I can't fix it immediately. But for now, stay away from range-based for on 2-d or higher ndarray iterators.
There was a problem hiding this comment.
Sorry, I don't understand: as far as I can tell for-each loops should work just fine on input iterators, since they only need inequality, increment, and dereference capability. They certainly don't need random access.
There was a problem hiding this comment.
I'll revert the change.
There was a problem hiding this comment.
Yes, for-each will work on either kind of iterator, but the compiler is allowed to do different things under the hood in the two cases, and the code it generates for random-access iterators may yield incorrect behavior if the iterator really only models input iterator. I have never really nailed down what is going on here (doing that is part of why I never got around to fixing it, though I really should just do it and investigate later), but I have seen at least one case (in an early version of afw::geom::SpanSet) where having an incorrect category for a proxy iterator led to some weird behavior in a range-based for, and that weird behavior went away if we either fixed the category or switched to a regular for loop (obviously the former is the right fix, but it always bothered me that I couldn't figure out what the range-based for was doing to cause trouble).
src/MatrixBuilder.cc
Outdated
| @@ -479,7 +479,7 @@ class RemappedShapeletImpl : public ShapeletImpl<T> { | |||
| _ellipse.scale(_radius); | |||
| _lhs.setZero(); | |||
| ShapeletImpl<T>::buildMatrix(_lhs, _ellipse); | |||
| output.asEigen() += _lhs.matrix() * _remapMatrix.transpose(); // undo the transpose in the ctor | |||
| ndarray::asEigenMatrix(output) += _lhs.matrix() * _remapMatrix.transpose(); // undo the transpose in ctor | |||
There was a problem hiding this comment.
Consider moving the comment to above this line instead of abbreviating it.
| result.getComponents().back().getCoefficients().asEigen() | ||
| = i->getMatrix().asEigen() * coefficients.asEigen(); | ||
| ndarray::asEigenMatrix(result.getComponents().back().getCoefficients()) | ||
| = ndarray::asEigenMatrix(i->getMatrix()) * ndarray::asEigenMatrix(coefficients); |
There was a problem hiding this comment.
I realize this is outside the scope of the ticket, but do you even need Eigen here?
There was a problem hiding this comment.
I think the eigen matrix map is being used for the matrix multiplication
src/ShapeletFunction.cc
Outdated
| @@ -169,13 +169,14 @@ ShapeletFunctionEvaluator::ShapeletFunctionEvaluator( | |||
| ShapeletFunction ShapeletFunction::convolve(ShapeletFunction const & other) const { | |||
| GaussHermiteConvolution convolution(getOrder(), other); | |||
| afw::geom::ellipses::Ellipse newEllipse(_ellipse); | |||
| auto matrix = convolution.evaluate(newEllipse).asEigen(); | |||
| auto matrixNdArray = convolution.evaluate(newEllipse); | |||
| auto matrix = ndarray::asEigenMatrix(matrixNdArray); | |||
There was a problem hiding this comment.
Consider adding a comment saying that matrixNdArray must be a local variable, to prevent accidental refactoring. Though since matrix is only used at one point, you could also just do the evaluation and conversion there.
There was a problem hiding this comment.
Good point. I'll eliminate the matrix temporary.
Replace ndarray::Array::asEigen, which returns an ndarray::EigenView, with ndarray::asEigen, which returns an Eigen::Map
Replace ndarray::Array::asEigen, which returns an ndarray::EigenView,
with ndarray::asEigen, which returns an Eigen::Map