Working notes - Indexing++
The goal of this page is to summarize the different ideas and working plan to (finally!) provide support for flexible row/column indexing in Eigen. See this bug report.
Needs
We aim to support various indexing mechanisms. For each dimension, we would like to be able to have any of:
- all: as in
A(:)
- singleton: as in
A(3)
(already supported!) - index-array: as in
A([3 1 5])
- range-based: as in
A(3:9)
(partly supported via Block) - slicing: as in
A(3:2:9)
- negative indices: as in
A(end-3)
orA(3:2:end-1)
- boolean masking: as in
A(A>0)
All of this with a concise and unambiguous API! We should also be able to provide compile-time information such as lengths.
C++11 API
Since achieving a lean API in C++03 seems to be rather impossible, let's first see what could be done in C++11 and provide more verbose fallbacks for C++03.
Better than a long story, here is a showcase demo:
#include <iostream> #include <array> #include <valarray> #include <vector> #include <Eigen/Dense> using namespace Eigen; using namespace std; struct Range { Range(Index f, Index l) : m_first(f), m_last(l) {} Index m_first, m_last; Index size() const { return m_last-m_first+1; } Index operator[] (Index k) { return m_first + k; } }; struct Slice { Slice(Index f, Index s, Index l) : m_first(f), m_step(s), m_last(l) {} Index m_first, m_step, m_last; Index size() const { return (m_last-m_first+m_step)/m_step; } Index operator[] (Index k) { return m_first + k*m_step; } }; template<typename T> struct has_boolean_value_type { enum { value = bool(internal::is_same<typename T::value_type,bool>::value) #if EIGEN_COMP_CLANG // Workaround clang's issue: // (valarray::operator<)::value_type is int instead of bool, // so we check the return type of (valarray::operator<).operator[], which is bool! || bool(internal::is_same<typename internal::remove_all<decltype(declval<T>()[0])>::type,bool>::value) #endif }; }; template<typename T> struct is_index { enum { value = std::is_integral<T>::value }; }; template<> struct is_index<bool> { enum { value = false }; }; template<typename T> struct has_index_value_type { enum { value = is_index<typename T::value_type>::value }; }; struct MyMat { void operator()(Range ids) { cout << "range-based: "; print_indices(ids); } void operator()(Slice ids) { cout << "slice-based: "; print_indices(ids); } template<typename Indices> typename internal::enable_if<has_index_value_type<Indices>::value>::type operator()(const Indices &ids) { cout << "index-based: "; print_indices(ids); } template<typename Mask> typename internal::enable_if<has_boolean_value_type<Mask>::value>::type operator()(const Mask &mask) { cout << "mask-based: "; for(int k=0; k<mask.size(); ++k) if(mask[k]) cout << k << " "; cout << "\n"; } protected: template<typename T> void print_indices(T& ids) { for(int k=0; k<ids.size(); ++k) cout << ids[k] << " "; cout << "\n"; } }; int main() { ArrayXd eia(10); eia.setRandom(); valarray<double> vala(10); Map<ArrayXd>(&vala[0],10) = eia; ArrayXi eii(4); eii << 3, 1, 7, 5; MyMat mat; mat({3,9}); // range mat({9,3}); // empty range (as in MatLab), need to use a slice: mat({9,-1,3}); // slice mat({9,-2,1}); // slice mat({3,2,9}); // slice mat({3,2,10}); // slice mat(valarray<bool>{true,false,false,true}); // mask mat(vector<bool>{true,false,false,true}); // mask mat(eia>0); // mask mat(vala>0.0); // mask mat(eii); // index-based mat(valarray<int>{5,2,5,6}); // index-based mat(array<int,4>{5,2,5,6}); // index-based }
that returns:
range-based: 3 4 5 6 7 8 9 range-based: slice-based: 9 8 7 6 5 4 3 slice-based: 9 7 5 3 1 slice-based: 3 5 7 9 slice-based: 3 5 7 9 mask-based: 0 3 mask-based: 0 3 mask-based: 2 4 7 8 9 mask-based: 2 4 7 8 9 index-based: 3 1 7 5 index-based: 5 2 5 6 index-based: 5 2 5 6
Still need to work on compile-time sizes, all, and end, and triple check for possible conflicts.
What does not work:
mat({5,2,6,7,1}); // needs to explicitly use a valarray, array, or ... mat({true,false,false,true,true,false,true}); // needs to explicitly use a valarray, array, or ...
C++03 fallback API
Internal implementation details
Performance
One could be tempted to leverage AVX2 gather/scatter instructions, but those are horribly slow. Better emulate them for now.