Difference between revisions of "Working notes - Indexing++"

From Eigen
Jump to: navigation, search
(Source code)
(Generated output)
Line 545: Line 545:
 
mat(all)                    -> all
 
mat(all)                    -> all
 
mat("")                      -> all
 
mat("")                      -> all
mat(eii)                    -> index-based: 3 1 6 5 [4]
 
mat(eii*2-1)                -> index-based: 5 1 11 9 [4]
 
 
mat(vali)                    -> index-based: 3 1 6 5 []
 
mat(vali)                    -> index-based: 3 1 6 5 []
 
mat(veci)                    -> index-based: 3 1 6 5 []
 
mat(veci)                    -> index-based: 3 1 6 5 []
 +
mat(eii)                    -> index-based: 3 1 6 5 [4]
 +
mat(eii*2-1)                -> index-based: 5 1 11 9 [4]
 
mat(eia>0)                  -> mask-based: 2 4 7 8 9 11 12
 
mat(eia>0)                  -> mask-based: 2 4 7 8 9 11 12
 
mat(vala>0.0)                -> mask-based: 2 4 7 8 9 11 12
 
mat(vala>0.0)                -> mask-based: 2 4 7 8 9 11 12
Line 565: Line 565:
 
mat(range(last-8,last-1,2))  -> strided-range: 4 6 8 10 []{}
 
mat(range(last-8,last-1,2))  -> strided-range: 4 6 8 10 []{}
 
mat(range(3,last-3,c<3>))    -> strided-range: 3 6 9 []{3}
 
mat(range(3,last-3,c<3>))    -> strided-range: 3 6 9 []{3}
 +
mat(range(3,last-3,c<3>()))  -> strided-range: 3 6 9 []{3}
 
mat(range(last-1,3,c<-2>))  -> strided-range: 11 9 7 5 3 []{-2}
 
mat(range(last-1,3,c<-2>))  -> strided-range: 11 9 7 5 3 []{-2}
 
mat(range(::end-1,3,c<-2>))  -> strided-range: 12 10 8 6 4 []{-2}
 
mat(range(::end-1,3,c<-2>))  -> strided-range: 12 10 8 6 4 []{-2}
Line 578: Line 579:
 
mat(span(9,3,-1))            -> strided-span: 9 8 7 []{}
 
mat(span(9,3,-1))            -> strided-span: 9 8 7 []{}
 
mat(span(9,3,-2))            -> strided-span: 9 7 5 []{}
 
mat(span(9,3,-2))            -> strided-span: 9 7 5 []{}
mat(span(9,c<3>,-2))   !14  -> strided-span: 9 7 5 [3]{}
+
mat(span(9,c<3>,-2))         -> strided-span: 9 7 5 [3]{}
mat(span(last,c<3>,-2)) !14  -> strided-span: 12 10 8 [3]{}
+
mat(span(last,c<3>,-2))     -> strided-span: 12 10 8 [3]{}
 
mat(span(last-1,3,-2))      -> strided-span: 11 9 7 []{}
 
mat(span(last-1,3,-2))      -> strided-span: 11 9 7 []{}
 
mat(span(last-1,3,c<-2>))    -> strided-span: 11 9 7 []{-2}
 
mat(span(last-1,3,c<-2>))    -> strided-span: 11 9 7 []{-2}
Line 600: Line 601:
 
mat(last-2*iota<3>())        -> strided-range: 12 10 8 []{}
 
mat(last-2*iota<3>())        -> strided-range: 12 10 8 []{}
 
mat(last-1-2*iota<3>())      -> strided-range: 11 9 7 []{}
 
mat(last-1-2*iota<3>())      -> strided-range: 11 9 7 []{}
mat(last-1-c<2>*iota<3>())  -> strided-range: 11 9 7 []{}
+
mat(9-2*iota(c<3>))         -> strided-range: 9 7 5 []{}
mat(9-2*iota(c<3>))     !14  -> strided-range: 9 7 5 []{}
+
 
mat(9-2*Iota<3>)        !14  -> strided-range: 9 7 5 []{}
 
mat(9-2*Iota<3>)        !14  -> strided-range: 9 7 5 []{}
 +
mat(last-1-c<2>*iota<3>()) !14  -> strided-range: 11 9 7 []{}
  
 
Some c++11 possibilities:
 
Some c++11 possibilities:
Line 609: Line 610:
 
mat({false,true,true,false})                -> mask-based: 1 2 []
 
mat({false,true,true,false})                -> mask-based: 1 2 []
 
mat(valarray<int>{5,2,5,6}*2-1)              -> index-based: 9 3 9 11 []
 
mat(valarray<int>{5,2,5,6}*2-1)              -> index-based: 9 3 9 11 []
mat(array<int,4>{5,2,5,6})                  -> index-based: 5 2 5 6 [4]
+
mat(array<int,4>{5,2,5,6})                  -> index-based: 5 2 5 6 []
 
mat( vector<bool>{false,true,true,false})    -> mask-based: 1 2
 
mat( vector<bool>{false,true,true,false})    -> mask-based: 1 2
 
mat(!valarray<bool>{false,true,true,false})  -> mask-based: 0 3
 
mat(!valarray<bool>{false,true,true,false})  -> mask-based: 0 3

Revision as of 12:54, 29 December 2016

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), but sometimes it is also convenient to specify the start+length similar to what Eigen::Block offers.
  • slicing: as in A(3:2:9). Same as for range, we sometime also want to specify start+length+stride
  • negative indices: as in A(end-3) or A(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 and maybe strides.

Some references

  • Some overviews what other programming languages do: [1] [2]
  • This post about the Range v3 library.
  • This comment describes an interesting API for range and slices inspired by APL.
  • libigl extends Eigen by a colon and a slice function
  • Several questions/feature request on StackOverflow: [3] [4] [5] [6]

API playground

To get some support for the discussions, let's start with a demo of a work-in-progress demo.

Source code

Click expand to see the implementation details.

#if __cplusplus >= 201103L || defined(_MSC_VER)
#define USE_CXX11
#endif
 
#include <iostream>
#ifdef USE_CXX11
#include <array>
#endif
#include <valarray>
#include <vector>
#include <cassert>
 
#define USE_EIGEN
 
#ifdef USE_EIGEN
#include <Eigen/Dense>
using namespace Eigen;
using internal::is_same;
#else
typedef std::ptrdiff_t Index;
static const int Dynamic = -1;
namespace internal {
template<typename T, int Value> class variable_if_dynamic
{
  public:
    explicit variable_if_dynamic(T v) { assert(v == T(Value)); }
    static T value() { return T(Value); }
    void setValue(T) {}
};
 
template<typename T> class variable_if_dynamic<T, Dynamic>
{
    T m_value;
    variable_if_dynamic() { assert(false); }
  public:
    explicit variable_if_dynamic(T value) : m_value(value) {}
    T value() const { return m_value; }
    void setValue(T value) { m_value = value; }
};
using std::is_same;
using std::enable_if;
using std::is_integral;
template<typename T> struct remove_all { typedef T type; };
template<typename T> struct remove_all<const T>   { typedef typename remove_all<T>::type type; };
template<typename T> struct remove_all<T const&>  { typedef typename remove_all<T>::type type; };
template<typename T> struct remove_all<T&>        { typedef typename remove_all<T>::type type; };
template<typename T> struct remove_all<T const*>  { typedef typename remove_all<T>::type type; };
template<typename T> struct remove_all<T*>        { typedef typename remove_all<T>::type type; };
}
#endif
 
using namespace std;
 
struct all_t { all_t() {} };
static const all_t all;
 
struct shifted_last {
  explicit shifted_last(int o) : offset(o) {}
  int offset;
  shifted_last operator+ (int x) const { return shifted_last(offset+x); }
  shifted_last operator- (int x) const { return shifted_last(offset-x); }
  int operator- (shifted_last x) const { return offset-x.offset; }
};
 
struct last_t {
  shifted_last operator- (int offset) const { return shifted_last(-offset); }
  shifted_last operator+ (int offset) const { return shifted_last(+offset); }
  int operator- (last_t) const { return 0; }
  int operator- (shifted_last x) const { return -x.offset; }
};
static const last_t last;
 
 
struct shifted_end {
  explicit shifted_end(int o) : offset(o) {}
  int offset;
  shifted_end operator+ (int x) const { return shifted_end(offset+x); }
  shifted_end operator- (int x) const { return shifted_end(offset-x); }
  int operator- (shifted_end x) const { return offset-x.offset; }
};
 
struct end_t {
  shifted_end operator- (int offset) const { return shifted_end (-offset); }
  shifted_end operator+ (int offset) const { return shifted_end ( offset); }
  int operator- (end_t) const { return 0; }
  int operator- (shifted_end x) const { return -x.offset; }
};
static const end_t end;
 
template<int N> struct Index_c {
  static const int value = N;
  operator int() const { return value; }
  Index_c (Index_c<N> (*)() ) {}
  Index_c() {}
  // Needed in C++14 to allow c<XX>():
  Index_c operator() () const { return *this; }
};
 
//--------------------------------------------------------------------------------
// Range(first,last) and Slice(first,step,last)
//--------------------------------------------------------------------------------
 
template<typename FirstType=Index,typename LastType=Index,typename StepType=Index_c<1> >
struct Range_t {
  Range_t(FirstType f, LastType l) : m_first(f), m_last(l) {}
  Range_t(FirstType f, LastType l, StepType s) : m_first(f), m_last(l), m_step(s) {}
 
  FirstType m_first;
  LastType  m_last;
  StepType  m_step;
 
  enum { SizeAtCompileTime = -1 };
 
  Index size() const { return (m_last-m_first+m_step)/m_step; }
  Index operator[] (Index k) const { return m_first + k*m_step; }
};
 
template<typename T> struct cleanup_slice_type { typedef Index type; };
template<> struct cleanup_slice_type<last_t> { typedef last_t type; };
template<> struct cleanup_slice_type<shifted_last> { typedef shifted_last type; };
template<> struct cleanup_slice_type<end_t> { typedef end_t type; };
template<> struct cleanup_slice_type<shifted_end> { typedef shifted_end type; };
template<int N> struct cleanup_slice_type<Index_c<N> > { typedef Index_c<N> type; };
template<int N> struct cleanup_slice_type<Index_c<N> (*)() > { typedef Index_c<N> type; };
 
template<typename FirstType,typename LastType>
Range_t<typename cleanup_slice_type<FirstType>::type,typename cleanup_slice_type<LastType>::type >
range(FirstType f, LastType l)  {
  return Range_t<typename cleanup_slice_type<FirstType>::type,typename cleanup_slice_type<LastType>::type>(f,l);
}
 
template<typename FirstType,typename LastType,typename StepType>
Range_t<typename cleanup_slice_type<FirstType>::type,typename cleanup_slice_type<LastType>::type,typename cleanup_slice_type<StepType>::type >
range(FirstType f, LastType l, StepType s)  {
  return Range_t<typename cleanup_slice_type<FirstType>::type,typename cleanup_slice_type<LastType>::type,typename cleanup_slice_type<StepType>::type>(f,l,typename cleanup_slice_type<StepType>::type(s));
}
 
 
template<typename T, int Default=-1> struct get_compile_time {
  enum { value = Default };
};
 
template<int N,int Default> struct get_compile_time<Index_c<N>,Default> {
  enum { value = N };
};
 
template<typename T> struct is_compile_time         { enum { value = false }; };
template<int N> struct is_compile_time<Index_c<N> > { enum { value = true }; };
 
template<typename FirstType=Index,typename SizeType=Index,typename StepType=Index_c<1> >
struct Span_t {
  Span_t(FirstType first, SizeType size) : m_first(first), m_size(size) {}
  Span_t(FirstType first, SizeType size, StepType step) : m_first(first), m_size(size), m_step(step) {}
 
  FirstType m_first;
  SizeType  m_size;
  StepType  m_step;
 
  enum { SizeAtCompileTime = get_compile_time<SizeType>::value };
 
  Index size() const { return m_size; }
  Index operator[] (Index k) const { return m_first + k*m_step; }
};
 
template<typename FirstType,typename SizeType,typename StepType>
Span_t<typename cleanup_slice_type<FirstType>::type,typename cleanup_slice_type<SizeType>::type,typename cleanup_slice_type<StepType>::type >
span(FirstType first, SizeType size, StepType step)  {
  return Span_t<typename cleanup_slice_type<FirstType>::type,typename cleanup_slice_type<SizeType>::type,typename cleanup_slice_type<StepType>::type>(first,size,step);
}
 
template<typename FirstType,typename SizeType>
Span_t<typename cleanup_slice_type<FirstType>::type,typename cleanup_slice_type<SizeType>::type >
span(FirstType first, SizeType size)  {
  return Span_t<typename cleanup_slice_type<FirstType>::type,typename cleanup_slice_type<SizeType>::type>(first,size);
}
 
//--------------------------------------------------------------------------------
// APL-like implementation based on iota(n) <=> {0,1,2,3,...n-1}
//--------------------------------------------------------------------------------
 
template<typename StartType> struct iota2slice {
  typedef Range_t<Index,Index,Index> type;
  typedef Index last_type;
};
template<> struct iota2slice<shifted_last> {
  typedef Range_t<shifted_last,shifted_last,Index> type;
  typedef shifted_last last_type;
};
 
template<int Size=Dynamic,int Step=1,typename StartType=Index_c<0> >
struct iota_t {
  iota_t(Index size=Size) : m_size(size), m_step(Step) {}
  iota_t(Index size,Index step,StartType start) : m_size(size), m_step(step), m_start(start) {}
  internal::variable_if_dynamic<Index,Size> m_size;
  internal::variable_if_dynamic<Index,Step> m_step;
  StartType m_start;
 
  typedef typename iota2slice<StartType>::type SliceType;
  typedef typename iota2slice<StartType>::last_type last_type;
 
  SliceType range() const {
    return SliceType(m_start,
                     m_start+(m_size.value()-1)*m_step.value(),
                     m_step.value());
  }
 
  operator SliceType() const { return range(); }
 
  // set a runtime step
  friend iota_t<Size,Dynamic,StartType> operator*(int stride,const iota_t &x) {
    return iota_t<Size,Dynamic,StartType>(x.m_size.value(),x.m_step.value()*stride,x.m_start);
  }
 
  // negate the step
  iota_t<Size,Step==Dynamic?Dynamic:-Step,StartType> operator-() {
    return iota_t<Size,Step==Dynamic?Dynamic:-Step,StartType>(m_size.value(),-m_step.value(),m_start);
  }
 
  // set first element
  iota_t<Size,Step,Index> operator+(int offset) const {
    return iota_t<Size,Step,Index>(m_size.value(),m_step.value(),m_start+offset);
  }
 
  // set first element
  friend iota_t<Size,Step,Index> operator+(int offset,const iota_t &x) {
    return iota_t<Size,Step,Index>(x.m_size.value(),x.m_step.value(),x.m_start+offset);
  }
 
  // set first index and negative step
  friend iota_t<Size,Step==Dynamic?Dynamic:-Step,Index> operator-(int offset,const iota_t &x) {
    return iota_t<Size,Step==Dynamic?Dynamic:-Step,Index>(x.m_size.value(),-x.m_step.value(),x.m_start+offset);
  }
 
  // set first index to last element and negate the step
  friend iota_t<Size,Step==Dynamic?Dynamic:-Step,shifted_last> operator-(last_t,const iota_t &x) {
    return iota_t<Size,Step==Dynamic?Dynamic:-Step,shifted_last>(x.m_size.value(),-x.m_step.value(),shifted_last(x.m_start));
  }
 
  // set first index to last element and negate the step
  friend iota_t<Size,Step==Dynamic?Dynamic:-Step,shifted_last> operator-(shifted_last l,const iota_t &x) {
    return iota_t<Size,Step==Dynamic?Dynamic:-Step,shifted_last>(x.m_size.value(),-x.m_step.value(),shifted_last(l.offset+x.m_start));
  }
};
 
iota_t<> iota(Index size) { return iota_t<>(size); }
template<int Size>
iota_t<Size> iota() { return iota_t<Size>(); }
 
#if __cplusplus > 201103L
template<int N>
static const Index_c<N> c{};
#else
template<int N>
inline Index_c<N> c() { return Index_c<N>(); }
#endif
 
#if __cplusplus > 201103L
template<int Size>
static const iota_t<Size> Iota{};
#endif
 
template<int Size>
iota_t<Size> iota(Index_c<Size>) { return iota_t<Size>(); }
 
//--------------------------------------------------------------------------------
// Matrix implementation
//--------------------------------------------------------------------------------
 
// NOTE that declaring value_type here is needed to make SFINAE works properly in case T does not exhibit value_type
// Another solution would be to specialize has_boolean_value_type and has_index_value_type for Range, Slice, etc.
template<typename T,typename value_type=typename T::value_type>
struct has_boolean_value_type {
  enum {
    value =  bool(internal::is_same<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 = internal::is_integral<T>::value }; };
template<> struct is_index<bool> { enum { value = false }; };
 
template<typename T,typename value_type=typename T::value_type>
struct has_index_value_type {
  enum {
    value = is_index<value_type>::value
  };
};
 
 
template<typename T, typename EnableIf = void> struct get_size {
  enum { value = -1 };
};
 
template<typename T> struct get_size<T,typename internal::enable_if<((T::SizeAtCompileTime&0)==0)>::type> {
  enum { value = T::SizeAtCompileTime };
};
 
#ifdef USE_CXX11
template<typename T,int N> struct get_size<std::array<T,N> > {
  enum { value = N };
};
#endif
 
struct MyMat {
 
  MyMat(Index n) : m_size(n) {}
 
  void operator()(Index k) { cout << "singleton: "; cout << k << '\n'; }
  void operator()(const char*) { cout << "all\n"; }
  void operator()(all_t) { cout << "all\n"; }
 
  template<typename FirstType,typename LastType,typename StepType>
  void operator()(Range_t<FirstType,LastType,StepType> ids) {
    cout << "strided-range: ";
    print_indices(Range_t<Index,Index,StepType>(getIndex(ids.m_first),getIndex(ids.m_last),ids.m_step));
    if(is_compile_time<StepType>::value)
      cout << "{" << get_compile_time<StepType,0>::value << "}\n";
    else
      cout << "{}\n";
  }
 
  template<typename FirstType,typename SizeType,typename StepType>
  void operator()(Span_t<FirstType,SizeType,StepType> ids) {
    cout << "strided-span: ";
    print_indices(Span_t<Index,SizeType,StepType>(getIndex(ids.m_first),ids.m_size,ids.m_step));
    if(is_compile_time<StepType>::value)
      cout << "{" << get_compile_time<StepType,0>::value << "}\n";
    else
      cout << "{}\n";
  }
 
  template<int Size,int Step,typename StartType>
  void operator()(iota_t<Size,Step,StartType> x) {
    return operator()(x.range());
  }
 
  template<typename Indices>
  typename internal::enable_if<has_index_value_type<Indices>::value>::type
  operator()(const Indices &ids) { cout << "index-based: "; print_indices(ids); cout << "\n"; }
 
  template<typename Mask>
  typename internal::enable_if<has_boolean_value_type<Mask>::value>::type
  operator()(const Mask &mask) {
    cout << "mask-based: ";
    for(Index k=0; k<Index(mask.size()); ++k)
      if(mask[k]) cout << k << " ";
    cout << "\n";
  }
 
  template<typename T,Index N>
  void operator()(const T (&ids)[N] ) {
    cout << "index-based: ";
    for(Index k=0; k<N; ++k)
      cout << ids[k] << " ";
    cout << " [" << N << "]\n";
  }
 
  template<Index N>
  void operator()(const bool (&mask)[N] ) {
    cout << "mask-based: ";
    for(Index k=0; k<N; ++k)
      if(mask[k]) cout << k << " ";
    cout << "[]\n";
  }
 
protected:
  Index getIndex(Index x) const       { return x; }
  Index getIndex(last_t) const        { return m_size-1; }
  Index getIndex(shifted_last x) const  { return m_size+x.offset-1; }
  Index getIndex(end_t) const        { return m_size; }
  Index getIndex(shifted_end x) const  { return m_size+x.offset; }
 
  template<typename T>
  void print_indices(const T& ids) {
    for(Index k=0; k<Index(ids.size()); ++k)
      cout << ids[k] << " ";
 
    int size = get_size<T>::value;
    if(size!=-1)  cout << "[" << size << "]";
    else          cout << "[]";
  }
 
  Index m_size;
};
 
#define PRINT(X,P) cout << #X << P << "  -> "; X
 
int main()
{
  int n = 13;
#ifdef USE_EIGEN
  ArrayXd eia(n); eia.setRandom();
  Array4i eii(4); eii << 3, 1, 6, 5;
  valarray<double> vala(n); Map<ArrayXd>(&vala[0],n) = eia;
  valarray<int> vali(4); Map<ArrayXi>(&vali[0],4) = eii;
  vector<int> veci(4); Map<ArrayXi>(veci.data(),4) = eii;
#else
  valarray<double> vala{3, 1, 6, 5};
  valarray<int> vali{3, 1, 6, 5};
  vector<int> veci{3, 1, 6, 5};
#endif
 
 
  MyMat mat(n);
  cout << "\nAll, indices, masking:\n"
       <<   "----------------------\n";
  PRINT( mat(5),"                     " );  // singleton
  PRINT( mat(all),"                   " );  // all
  PRINT( mat(""),"                    " );  // all
 
  PRINT( mat(vali),"                  " );  // index-based
  PRINT( mat(veci),"                  " );  // index-based
#ifdef USE_EIGEN
  PRINT( mat(eii),"                   " );  // index-based
  PRINT( mat(eii*2-1),"               " );  // index-based
  PRINT( mat(eia>0),"                 " );  // mask
#endif
  PRINT( mat(vala>0.0),"              " );  // mask
 
 
  cout << "\nLower-upper-bound API using range:\n"
       <<   "----------------------------------\n";
  PRINT( mat(range(3,9)),"            " );
  PRINT( mat(range(3,last)),"         " );
  PRINT( mat(range(3,last-2)),"       " );
  PRINT( mat(range(9,3)),"            " );  // empty range (as in MatLab), need to use a slice:
  PRINT( mat(range(9,3,-1)),"         " );
  PRINT( mat(range(9,1,-2)),"         " );
  PRINT( mat(range(last,3,-2)),"      " );
  PRINT( mat(range(last-1,3,-2)),"    " );
  PRINT( mat(range(3,last-3,3)),"     " );
  PRINT( mat(range(last-8,last-1,2)),"" );
  PRINT( mat(range(3,last-3,c<3>)),"  " );
  PRINT( mat(range(3,last-3,c<3>())),"" );
  PRINT( mat(range(last-1,3,c<-2>))," " );
  PRINT( mat(range(::end-1,3,c<-2>)),"" );
  cout << "The even elements among the last six in ascending order:\n";
  PRINT( mat(range(last-6,last,2)),"    " );
  PRINT( mat(range(::end-7,::end-1,2)),"" );
 
 
  cout << "\nLength-based API using span:\n"
       <<   "----------------------------\n";
  PRINT( mat(span(0,3)),"             " );
  PRINT( mat(span(2,3)),"             " );
  PRINT( mat(span(3,3,2)),"           " );
  PRINT( mat(span(9,3,-1)),"          " );
  PRINT( mat(span(9,3,-2)),"          " );
  PRINT( mat(span(9,c<3>,-2)),"       " );
  PRINT( mat(span(last,c<3>,-2)),"    " );
  PRINT( mat(span(last-1,3,-2)),"     " );
  PRINT( mat(span(last-1,3,c<-2>)),"  " );
  PRINT( mat(span(1,c<3>,c<2>)),"     " );
  cout << "The last m=4 even (or odd) elements in ascending order:\n";
  int m=4;
  PRINT( mat(span(last-2*(m-1),m,2)),"" );
  PRINT( mat(span(::end+1-2*m,m,2))," " );
  cout << "The last m=4 elements with stride s=3 in ascending order:\n";
  int s=3;
  PRINT( mat(span(last-s*(m-1),m,s))," " );
  PRINT( mat(span(::end-1+s-s*m,m,s)),"" );
 
 
  cout << "\nLength-based API using iota:\n"
       <<   "----------------------------\n";
  PRINT( mat(iota(3)),"               " );
  PRINT( mat(2+iota(3)),"             " );
  PRINT( mat(3+2*iota(3)),"           " );
  PRINT( mat(9-iota(3)),"             " );
  PRINT( mat(9-2*iota(3)),"           " );
  PRINT( mat(9-2*iota<3>()),"         " );
  PRINT( mat(last-2*iota<3>()),"      " );
  PRINT( mat(last-1-2*iota<3>()),"    " );
#if __cplusplus > 201103L
  PRINT( mat(9-2*iota(c<3>)),"        " );
  // int c; // check naming collision (ok with gcc, ko with clang)
  PRINT( mat(9-2*Iota<3>),"        !14" );
  PRINT( mat(last-1-c<2>*iota<3>())," !14" );
#endif
 
#ifdef USE_CXX11
  cout << "\nSome c++11 possibilities:\n"
       << "---------------------------\n";
  PRINT( mat({5,2,5,6}),"                             " );  // index-based
  PRINT( mat({false,true,true,false}),"               " );  // mask
  PRINT( mat(valarray<int>{5,2,5,6}*2-1),"            " );  // index-based
  PRINT( mat(array<int,4>{5,2,5,6}),"                 " );  // index-based
  PRINT( mat( vector<bool>{false,true,true,false}),"  " );  // mask
  PRINT( mat(!valarray<bool>{false,true,true,false}),"" );  // mask
#endif
}

Generated output

The values on the right indicate the extracted compile time size and step as: [SIZE]{STEP}. Here, [] means a runtime size, and {} a runtime step.

All, indices, masking:
----------------------
mat(5)                       -> singleton: 5
mat(all)                     -> all
mat("")                      -> all
mat(vali)                    -> index-based: 3 1 6 5 []
mat(veci)                    -> index-based: 3 1 6 5 []
mat(eii)                     -> index-based: 3 1 6 5 [4]
mat(eii*2-1)                 -> index-based: 5 1 11 9 [4]
mat(eia>0)                   -> mask-based: 2 4 7 8 9 11 12
mat(vala>0.0)                -> mask-based: 2 4 7 8 9 11 12

Lower-upper-bound API using range:
----------------------------------
mat(range(3,9))              -> strided-range: 3 4 5 6 7 8 9 []{1}
mat(range(3,last))           -> strided-range: 3 4 5 6 7 8 9 10 11 12 []{1}
mat(range(3,last-2))         -> strided-range: 3 4 5 6 7 8 9 10 []{1}
mat(range(9,3))              -> strided-range: []{1}
mat(range(9,3,-1))           -> strided-range: 9 8 7 6 5 4 3 []{}
mat(range(9,1,-2))           -> strided-range: 9 7 5 3 1 []{}
mat(range(last,3,-2))        -> strided-range: 12 10 8 6 4 []{}
mat(range(last-1,3,-2))      -> strided-range: 11 9 7 5 3 []{}
mat(range(3,last-3,3))       -> strided-range: 3 6 9 []{}
mat(range(last-8,last-1,2))  -> strided-range: 4 6 8 10 []{}
mat(range(3,last-3,c<3>))    -> strided-range: 3 6 9 []{3}
mat(range(3,last-3,c<3>()))  -> strided-range: 3 6 9 []{3}
mat(range(last-1,3,c<-2>))   -> strided-range: 11 9 7 5 3 []{-2}
mat(range(::end-1,3,c<-2>))  -> strided-range: 12 10 8 6 4 []{-2}
The even elements among the last six in ascending order:
mat(range(last-6,last,2))      -> strided-range: 6 8 10 12 []{}
mat(range(::end-7,::end-1,2))  -> strided-range: 6 8 10 12 []{}

Length-based API using span:
----------------------------
mat(span(0,3))               -> strided-span: 0 1 2 []{1}
mat(span(2,3))               -> strided-span: 2 3 4 []{1}
mat(span(3,3,2))             -> strided-span: 3 5 7 []{}
mat(span(9,3,-1))            -> strided-span: 9 8 7 []{}
mat(span(9,3,-2))            -> strided-span: 9 7 5 []{}
mat(span(9,c<3>,-2))         -> strided-span: 9 7 5 [3]{}
mat(span(last,c<3>,-2))      -> strided-span: 12 10 8 [3]{}
mat(span(last-1,3,-2))       -> strided-span: 11 9 7 []{}
mat(span(last-1,3,c<-2>))    -> strided-span: 11 9 7 []{-2}
mat(span(1,c<3>,c<2>))       -> strided-span: 1 3 5 [3]{2}
The last m=4 even (or odd) elements in ascending order:
mat(span(last-2*(m-1),m,2))  -> strided-span: 6 8 10 12 []{}
mat(span(::end+1-2*m,m,2))   -> strided-span: 6 8 10 12 []{}
The last m=4 elements with stride s=3 in ascending order:
mat(span(last-s*(m-1),m,s))   -> strided-span: 3 6 9 12 []{}
mat(span(::end-1+s-s*m,m,s))  -> strided-span: 3 6 9 12 []{}

Length-based API using iota:
----------------------------
mat(iota(3))                 -> strided-range: 0 1 2 []{}
mat(2+iota(3))               -> strided-range: 2 3 4 []{}
mat(3+2*iota(3))             -> strided-range: 3 5 7 []{}
mat(9-iota(3))               -> strided-range: 9 8 7 []{}
mat(9-2*iota(3))             -> strided-range: 9 7 5 []{}
mat(9-2*iota<3>())           -> strided-range: 9 7 5 []{}
mat(last-2*iota<3>())        -> strided-range: 12 10 8 []{}
mat(last-1-2*iota<3>())      -> strided-range: 11 9 7 []{}
mat(9-2*iota(c<3>))          -> strided-range: 9 7 5 []{}
mat(9-2*Iota<3>)        !14  -> strided-range: 9 7 5 []{}
mat(last-1-c<2>*iota<3>()) !14  -> strided-range: 11 9 7 []{}

Some c++11 possibilities:
---------------------------
mat({5,2,5,6})                               -> index-based: 5 2 5 6  [4]
mat({false,true,true,false})                 -> mask-based: 1 2 []
mat(valarray<int>{5,2,5,6}*2-1)              -> index-based: 9 3 9 11 []
mat(array<int,4>{5,2,5,6})                   -> index-based: 5 2 5 6 []
mat( vector<bool>{false,true,true,false})    -> mask-based: 1 2
mat(!valarray<bool>{false,true,true,false})  -> mask-based: 0 3

Discussions

nD-Arrays

For concision, this demo only considers the 1D case, but of course the idea is to implement it in the nD settings. For instance, with a 3D tensor, we could imagine:

mat(span(0,2,last), "", vec>0)

which is why an API like mat.slice(...) is a no go. This is also a crucial limitation of the current Eigen::Block-based API.

Index-based

This case is relatively easy to handle. We mostly have to check for the value_type of the input, and extract the compile-time size if available. The input must provide a size() method and operator[](int). This way we can support Eigen::Matrix, std::vector, std::valarray, std::array, etc.

Mask-based

Masking is similar to the previous Index-based case except that we have to check whether value_type==bool and workaround some clang issues to deal with the return type of the comparison operators of valarray. Here we probably do not want to support multi-dimentional masking as in matlab: A=rand(n,n); A(A<0)=0;. There is DenseBase::select for that purpose. Select could be extended to be a lvalue too.

Slicing

This is the most tricky part as there exist many ways to express the same thing. Let's start with some definitions:

  • start - the index of the first element in the slice
  • stop - either the index of the last element in the slice (inclusive), or one more than the index of last element in the slice (exclusive)
  • len - the length of the slice (= exclusive(stop) - first = inclusive(stop) - first + 1)
  • step - the increment to go from one element to the next one (default 1)

Among those 5 quantities (because stop has two variants), only three are necessary to define a slice. The cases where step can be inferred from the others is already handled by LinSpace and does not appear to be relevant for indexing, so let's focus on the other cases where step is given either explicitly or implicitly as ti defaults to 1. We can distinguish the following cases:

  1. Explicit-length:
    1. [start,len]
    2. [inclusive(stop),len]
    3. [exclusive(stop),len]
  2. Lower-upper-bounds:
    1. Inclusive upper-bound: [start,inclusive(stop)]
    2. Exclusive upper-bound: [start,exclusive(stop)]

Those five cases are equivalent in theory but depending on the context one might be more convenient than the other, and since being able to define compile-time lengths is crucial for performance reason, we thus must support at least one variant with explicit-length and one with lower-upper-bounds. On the other hand, there is probably no need to strive to support compile-time-size with 'lower-upper-bounds'. The rationale is that if the length is known at compile-time then there is a high chance that the user better knows the length than the index bounds.

In order to define relative indices from the end, we also need to introduce pseudo keywords referring to some notion of the end of the underlying array. Again, we have to consider two possible variants:

  1. last refers to the index of the last element (inclusive)
  2. end refers to the index of the last element + 1 (exclusive)

In the following, we'll use upper case to denote compile-time quantities such as LEN and STEP.

So now the challenge is to provide an unambiguous and concise API to specify at least one of the explicit-length (with optional compile-time length and step) and one lower-upper-bounds variant, with, in both cases, an optional step. A first step consist in defining concise names for these two paradigms. I propose to use span to refer to explicit-length case, and range to refer to the lower-upper-bounds variant. Here is the rationale:

  • range is (for me) more related to the notions of interval, limits, gamut, etc. that are naturally defined by their bounds.
  • span is related to the notion of period of time, distance, width, extent, etc. and thus the notion of length here.

The word slice then refers to the generic notion.


Inclusive vs exclusive

At this point it is pretty clear that in order to reduce the ambiguities we will have to choose among the inclusive/exclusive variants of stop. At a first glance, we can make the following remarks:

  • Advantage of inclusive(stop)
    • closer to matlab
    • symmetry with first
    • permits to introduce a keyword Eigen::last to refer to the last element of the underlying array whereas the exclusive approach would suggest introducing a Eigen::end that conflicts with std::end.
  • Advantage of exclusive(stop)
    • closer to c++ begin/end paradigm.

More pros/cons will come later.

Span (explicit-length)

For explicit-length spans, a naive approach like:

span(start,len[,step]);
span<LEN>(start[,step]);
span<STEP>(start,len);
span<LEN,STEP>(start);

is probably not the way to go because the argument orders are not consistent across all compile-time/optional variants, and the second and third cases are ambiguous. This approach would thus require to introduce much more verbose function names to disambiguate the different possibilities.

In the previous demo we proposed a completely different approach that incrementally builds a complex slice starting from a uniform list of indices via iota(n) that generates n integers from 0 to n-1. Then we can apply a scaling factor and an offset to define the step and start respectively. The syntax is:

[start +] [step*]iota(len)

For compile-time length:

[start +] [step*]iota<LEN>()    // c++98
[start +] [step*]iota<LEN>      // c++14, perhaps possible with c++11

For compile-time step (not implemented in the demo):

[start +] Index_c<STEP>()*iota(len)    // c++98
[start +] c<STEP>()*iota(len)       // c++98 with shortcut
[start +] c<STEP>*iota(len)      // c++14, perhaps possible with c++11

And you can of course combine compile-time step and length without ambiguity. This API also allows to seamlessly index from the end of the array:

// the last n elements in descending order
last-iota(n);         
end-1-iota(n);
// the last n elements in ascending order
last-(n-1)+iota(n);
end-n+iota(n); 
// the last n even (or odd) elements in ascending order
last-2*(n-1)+2*iota(n);   
last-2*(n-1-iota(n));
end+1-2*n+2*iota(n);

Here we see that using an inclusive last keyword is more convenient that an exclusive end keyword. Perhaps we could imagine adding some manipulation functions, like flip...


Taking inspiration from the Index_c<STEP>() or c<STEP> approach to specify a quantity at compile time, we could revisit the previous span(start,len[,step]); with:

span(start, c<LEN>)
span(start, c<LEN>, step)
span(start, len, c<STEP>)
span(start, c<LEN>, c<STEP>)

no ambiguity anymore. Of course, for concision this assumes a c++14 compiler, and the user has to import Eigen::c in the current scope. The long version could be Eigen::Index_c<len>(). And for indexing from the end:

// the last n elements in descending order
span(last, len, -1)
span(end-1, len, -1)
// the last n elements in ascending order
span(last-(len-1), len)
span(end-len, len)
// the last n even (or odd) elements in ascending order
span(last-2*(len-1), len, 2)
span(end+1-2*n, len, 2)

Range (lower-upper-bounds)

TODO


Others

  1. Using "" for all might looks odd, but it would avoid introducing a polluting identifier all (nobody wants to write Eigen::all). We have the same issue with last though.
  2. APL-like API:
    • Pros:
      • easy to specify compile-time info
      • no argument-order confusion
      • can be built on top of the range/slice API (hence, can be done outside Eigen)
    • Cons:
      • cumbersome?
      • needs to specify length instead of first-last

Internal implementation details

Performance

One could be tempted to leverage AVX2 gather/scatter instructions, but those are horribly slow. Better emulate them for now.