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

From Eigen
Jump to: navigation, search
(Some references)
(Q: Inclusive vs exclusive?)
Line 684: Line 684:
  
 
The symmetry advantage of the inclusive variant is probably strong enough to favor this choice.
 
The symmetry advantage of the inclusive variant is probably strong enough to favor this choice.
 +
 +
Now let's look at examples:
 +
 +
{| class="wikitable"
 +
|-
 +
! intents !! aseq (inclusive) !! aseqn !! aseqX (exclusive) !! python
 +
|-
 +
| first k elements || <tt>A(aseq(0,k+1))</tt> || <tt>A(aseqn(0,k)) </tt> || <tt>A(aseqX(0,k))</tt> || <tt>A[:k]</tt>
 +
|-
 +
| first k elements (descending)  || <tt>A(aseq(k-1,0,fix<-1>))</tt> || <tt>A(aseqn(k-1,k,fix<-1>)) </tt> || <tt>A(aseqX(k-1,-1,fix<-1>))</tt> || <tt>A[k-1::-1]</tt>
 +
|-
 +
| last k elements (ascending) || <tt>A(aseq(last-k+1,last))</tt> <br /> <tt>A(aseq(end-k,last))</tt> || <tt>A(aseqn(last-k+1,k)) </tt> <br/> <tt>A(aseqn(end-k,k)) </tt> || <tt>A(aseqX(last-k+1,last+1))</tt><br/><tt>A(aseqX(end-k,end))</tt> || <tt>A[-k:]</tt>
 +
|-
 +
|  last k elements (descending) || <tt>A(aseq(last,last-k+1,fix<-1>))</tt> <br/>  <tt>A(aseq(last,end-k,fix<-1>))</tt> || <tt>A(aseqn(last,k,fix<-1>)) </tt> || <tt>A(aseqX(last,last-k,fix<-1>))</tt> || <tt>A[:-(k+1):-1]</tt>
 +
|-
 +
|  elements || <tt>A(aseq(,))</tt> || <tt>A(aseqn(,)) </tt> || <tt>A(aseqX(,))</tt> || <tt>A[:]</tt>
 +
|-
 +
|  elements || <tt>A(aseq(,))</tt> || <tt>A(aseqn(,)) </tt> || <tt>A(aseqX(,))</tt> || <tt>A[:]</tt>
 +
|-
 +
|  elements || <tt>A(aseq(,))</tt> || <tt>A(aseqn(,)) </tt> || <tt>A(aseqX(,))</tt> || <tt>A[:]</tt>
 +
|-
 +
|  elements || <tt>A(aseq(,))</tt> || <tt>A(aseqn(,)) </tt> || <tt>A(aseqX(,))</tt> || <tt>A[:]</tt>
 +
|-
 +
|  elements || <tt>A(aseq(,))</tt> || <tt>A(aseqn(,)) </tt> || <tt>A(aseqX(,))</tt> || <tt>A[:]</tt>
 +
|-
 +
|  elements || <tt>A(aseq(,))</tt> || <tt>A(aseqn(,)) </tt> || <tt>A(aseqX(,))</tt> || <tt>A[:]</tt>
 +
|-
 +
|  elements || <tt>A(aseq(,))</tt> || <tt>A(aseqn(,)) </tt> || <tt>A(aseqX(,))</tt> || <tt>A[:]</tt>
 +
|-
 +
|}
  
 
==Function-based API==
 
==Function-based API==

Revision as of 09:27, 5 January 2017

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

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.

Arithmetic progressions (slicing)

This is the most tricky part as there exist many ways to define the same arithmetic progression (aka arithmetic sequence). 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)
  • size - the length of the slice (= exclusive(stop) - first = inclusive(stop) - first + 1)
  • incr - 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 incr 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 incr is given either explicitly or implicitly as ti defaults to 1. We can distinguish the following cases:

  1. Explicit-size:
    1. [start,size]
    2. [inclusive(stop),size]
    3. [exclusive(stop),size]
  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 sizes is crucial for performance reason, we thus must support at least one variant with explicit-size 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 size is known at compile-time then there is a high chance that the user better knows the size 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 SIZE and INCR.

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

  • range is 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 size here.

The word slice then refers to the generic notion.


Q: 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.
    • Used in several languages/libraries: MatLab, Mathematica, Fortran, R, Pascal, Ada, Perl, boost::irange, etc.
  • Advantage of exclusive(stop)
    • closer to c++ begin/end paradigm.
    • Used in several langages/libraries: python, rust, GO, D, C++ range-v3, etc.

The symmetry advantage of the inclusive variant is probably strong enough to favor this choice.

Now let's look at examples:

intents aseq (inclusive) aseqn aseqX (exclusive) python
first k elements A(aseq(0,k+1)) A(aseqn(0,k)) A(aseqX(0,k)) A[:k]
first k elements (descending) A(aseq(k-1,0,fix<-1>)) A(aseqn(k-1,k,fix<-1>)) A(aseqX(k-1,-1,fix<-1>)) A[k-1::-1]
last k elements (ascending) A(aseq(last-k+1,last))
A(aseq(end-k,last))
A(aseqn(last-k+1,k))
A(aseqn(end-k,k))
A(aseqX(last-k+1,last+1))
A(aseqX(end-k,end))
A[-k:]
last k elements (descending) A(aseq(last,last-k+1,fix<-1>))
A(aseq(last,end-k,fix<-1>))
A(aseqn(last,k,fix<-1>)) A(aseqX(last,last-k,fix<-1>)) A[:-(k+1):-1]
elements A(aseq(,)) A(aseqn(,)) A(aseqX(,)) A[:]
elements A(aseq(,)) A(aseqn(,)) A(aseqX(,)) A[:]
elements A(aseq(,)) A(aseqn(,)) A(aseqX(,)) A[:]
elements A(aseq(,)) A(aseqn(,)) A(aseqX(,)) A[:]
elements A(aseq(,)) A(aseqn(,)) A(aseqX(,)) A[:]
elements A(aseq(,)) A(aseqn(,)) A(aseqX(,)) A[:]
elements A(aseq(,)) A(aseqn(,)) A(aseqX(,)) A[:]

Function-based API

In this section, we will discuss how to design an API of the form:

Matrix A;
auto subA = A(foo(...), foo(...));

with some proper names for the function foo and its argument.

Variant A: two function names

Here we'll start with the variant A that use two different names for the bound- and size- based versions.

Arithmetic progressions with explicit-size (span)

For explicit-size spans, a naive approach like:

span(start,size[,incr]);
span<SIZE>(start[,incr]);
span<INCR>(start,size);
span<SIZE,INCR>(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 anyways. This approach would thus require to introduce much more verbose function names to disambiguate the different possibilities.

In order to preserve the order of the arguments, we could pass compile-time values through the notion of integral_constant<N>, that is currently abbreviated as ic. The previous span(start,size[,incr]); example, can then be revisited as follows:

span(start, ic<SIZE>)
span(start, ic<SIZE>, incr)
span(start, size, ic<INCR>)
span(start, ic<SIZE>, ic<INCR>)

No ambiguity anymore! Of course, for concision this assumes the user has imported a short Eigen::ic in the current scope. And for indexing from the end:

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

Arithmetic progressions with lower-upper-bounds (range)

This case is simpler because there is fewer variants. Following the previous span variants, we coudl propose:

range(start, stop)
range(start, stop, incr)
range(start, stop, ic<INCR>)

Here, if we agree that only the increment can be specified at compile-time, a template overload like:

range<INCR>(start, stop)

would also be acceptable.

Variant B: single function name for arithmetic progressions

The previous variants needs two different names (currently, range and span) to refer to the same notion of arithmetic progressions. It is possible to unify them using named arguments. For now, let's call it aseq as a contraction for arithmetic sequence. Then we could have:

aseq(start,stop[,incr])              <-> range(start,stop[,incr])
aseq(start,stop,Incr<INCR>)          <-> range(start,stop,ic<incr>)
aseq(start, Size(size)[,incr])       <-> span(start,size[incr])
aseq(start, Size<SIZE>[,incr])       <-> span(start, ic<SIZE>[, incr])
aseq(start, Size(size),Incr<INCR>)   <-> span(start, size, ic<INCR>)
aseq(start, Size<SIZE>,Incr<INCR>)   <-> span(start, ic<SIZE>, ic<INCR>)

The main drawbacks is that it implies a weird API for the simple aseq(start, Size(size)) case, and it is less generalizable as it requires a new identifier in Eigen's namespace for each argument name, whereas ic<N> can be used everywhere.

Variant C: expression-based API

Since finding short and clear names is difficult, in this section we will discuss a completely different approach. The idea is to programmatically define an arithmetic progression using standard arithmetic/comparison operators, and a single iota function generating an arithmetic progression starting at 0 and with unit increment:

Iota      <-->  [0,1,2,...,last]
iota(n)   <-->  [0,1,2,...,n-1], n==size/length of the sequence
iota<N>   <-->   same but with compile-time size

By default, the sequence is thus bounded by the underlying matrix sizes, and equivalent to the notion of all. The default could also be spelled iota<> or iota(), or perhaps using a completely different name...

Then the idea consists in using arithmetic operators and comparison operators to set the start, stop and increment as detailed below.

With an explicit size

Here we start with the iota(n)/iota<N> versions that generates n integers from 0 to n-1. Then we can apply a scaling factor and an offset to define the incr and start respectively. The syntax is:

[start +] [incr*]iota(size)

For compile-time size:

[start +] [incr*]iota<SIZE>()    // c++98
[start +] [incr*]iota<SIZE>      // c++14

For a compile-time increment, we could use the notion of integral_constant<N>, that is, as before, currently abbreviated as ic:

[start +] ic<INCR>()*iota(size)    // c++98
[start +] ic<INCR>*iota(size)      // c++98 with shortcut

And you can of course combine compile-time incr and size 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 or reverse...

With lower-upper-bounds

Here we start with the default generator Iota that is unbounded, or rather implicitly bounded by the underlying matrix sizes once the expression is passed to the matrix. Then a scaling factor can be applied to define a non-default increment:

// All entries:
A(Iota)
// all even entries
A(2*Iota)
// All entries in reverse order (incr==-1 at compile-time)
A(-Iota)
// All even entries in reverse order (incr==-2 at compile-time)
A(2*-Iota)
A(ic<-2>*Iota)

Note that we can define a compile time -1 increment without messing with ic<-1>.

Then explicit bounds (either inclusive or exclusive) can be defined using comparison operators:

A(Iota<=6)            // the first 7 entries [0..6]  <--> A(range(ic<0>,6))
A(Iota<7)             // the first 7 entries         <--> A(range(ic<0>,6))
A(3<=Iota<=9)         // [3..9]                      <--> A(range(3,9))
A(3<=ic<2>*Iota<10)   // [3,5,7,9]                   <--> A(range(3,9,ic<2>))
A(3<=Iota)            // [3,...]                     <--> A(range(3,last))

Not too bad, except the 4th line A(3<=ic<2>*Iota<10) that is getting messy with all these < and > symbols.

Choice of identifier names

Now comes the difficult choice to pick good names for each novel identifiers. A good name is concise, clear, and not too polluting.

Integral constant For variants A and C
int_constant<N> Too verbose!
ic<N> A bit short and unclear the first time you see it, but sounds pretty good.
c<N> Even shorter! Clearly too polluting.
Fixed<N> Pretty good too.
Fix<N>
Cst<N>
IndexC<N> Reference to Eigen::Index
CIndex<N>
CstIndex<N>
Static<N> static is a C++ keyword
Stat<N> Shortcut for Static
(int[N]){} A name-free variant, just for fun, does not work with N<1
Lower-upper Lower-size For variant A only
range(i,j) span(i,n)
fromto(i,j) fromfor(i,n) fromto looks quite nice, but not fromfor.
fromto(i,j) fromlen(i,n) not really better than fromfor.
aseq(i,j) aseqn(i,n) suffix -n for the size-based version, sounds pretty good.
seq(i,j) seqn(i,n) sequence alone is a bit too general
fap(i,j) fapn(i,n)
range(i,j) rangen(i,n)
Arithmetic progression/sequence For variant B only.
aseq(...)
seq(...)
aprog(...)
fap(...) finite arithmetic progression
smseq(...) strictly monotonic sequence
range(...)
span(...)
Size For variant B only
Size<N> Very clear, but polluting
Length<N>
Len<N>
Count<N>
Card<N>
Increment For variant B only
Incr<N> Short and precise
Stride<N> Already used in Eigen (Stride<Outer,Inner> <-> Stride<LeadingDim,Incr>)
Step<N> step can mean a lot of different things
Full range (aka : in Matlab) For variants A and B
All
":" catch char* at compile-time, with runtime check in debug mode
Compile-time 0 For variants A, B, C
Zero
"0" catch char* at compile-time, with runtime check in debug mode

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 size 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.