Working notes - Indexing++

From Eigen
Jump to: navigation, search

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 seq (inclusive) or seqn seqX (exclusive) python expression-based with comparison operator expr with predicate
(when appropriate)
first k elements A(seqn(0,k)) A(seqX(0,k)) A[:k] A(iota(k))
A(Iota<k)
first k elements (descending) A(seqn(k-1,k,fix<-1>))
A(seq(k-1,0,fix<-1>))
A(seqX(k-1,-1,fix<-1>)) A[k-1::-1] A(-Iota<k)
A(k-1-iota(k))
first elements to index i A(seq(0,i)) A(seqX(0,i+1)) A[:i+1] A(iota(i+1))
A(Iota<=i)
first elements to index i (descending) A(seq(i,0,fix<-1>))
A(seqn(i,i+1,fix<-1>))
A(seqX(i,-1,fix<-1>)) A[i::-1] A(i-iota(i+1))
A(-Iota<=i)
A(i>=-Iota)
last k elements (ascending) A(seq(end-k,last))
A(seqn(last-k+1,k))
A(seqn(end-k,k))
A(seqX(last-k+1,last+1))
A(seqX(end-k,end))
A[-k:] A(last-k<Iota)
A(reverse(last-iota(k)))
last k elements
(descending)
A(seqn(last,k,fix<-1>))
A(aseq(last,end-k,fix<-1>))
A(seqX(last,last-k,fix<-1>)) A[:-(k+1):-1] A(last-iota(k))
A(last-k<-Iota)
all but the last k A(seq(0,last-k))
A(seqn(0,size-k))
A(seqX(0,last-k+1)) A[:-k] A(iota(last+1-k))
A(Iota<=last-k)
A(Iota<end-k)
A(Iota|ubound(last-k))
A(Iota|_1<=last-k)
all but the last k
(descending)
A(seq(last-k,0,fix<-1>))
A(seqn(last-k,size-k,fix<-1>))
A(seqX(last-k,-1,fix<-1>)) A[-(k+1)::-1] A(-Iota<=last-k)
A(-Iota<end-k)
A(last-k-Iota(last+1-k))
subset of the first k
elements with step s
A(seq(0,k-1,s))
A(seqn(0,(k+s-1)/s,s))
A(seqX(0,k,s)) A[:k:s] A(Iota*s<k) A(Iota*s|ubound(k-1))
A(Iota*s|_1<k)
subset of the first k
elements with step s
(descending)
A(seq(k-1,fix<0>,-s))
A(seqn(k-1,(k+s-1)/s,-s))
A(seqX(k-1,-1,-s)) A(k-1-Iota*s)
A[k-1::-s]
A(-Iota*s<k) A(k-1-Iota*s|lbound<0>)
A(k-1-Iota*s|_1>=0)
pick last k elements with
step s (ascending)
A(seq(last-(k-1)*s,last,s))
A(seqn(last-(k-1)*s,k,s))
A(flip(seqn(last,k,-s)))
A(seqX(last-(k-1)*s,end,s)) A[-(k-1)*s-1::s] A(flip(last-iota(k)*s))
A(last-(k-1)*s+iota(k)*s)
pick last k elements with
step s (descending)
A(seq(last,end-(k-1)*s,-s))
A(seqn(last,k,-s))
A(seqX(last,last-(k-1)*s-1,-s)) A[:-(k-1)*s-2:-s] A(last-iota(k)*s)
elements A(seq(,)) A(seqn(,)) A(seqX(,)) A[:] A()
elements A(seq(,)) A(seqn(,)) A(seqX(,)) A[:] A()
elements A(seq(,)) A(seqn(,)) A(seqX(,)) A[:] A()
elements A(seq(,)) A(seqn(,)) A(seqX(,)) A[:] 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 (dismissed)

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.

Caveats

  • The syntax seq(a,b,incr) with inclusive bounds a and b looks symmetric regarding the bounds but they are not:
    • The sequence starts at a, but in general it is not guaranteed that the last element is b. The last element is: a+((b-a)/incr)*incr
    • For the same reason, the reverse of seq(a,b,incr) is not 'seq(b,a,-incr) but seq('a+((b-a)/incr)*incr,a,-incr)
    • As a consequence, we cannot call these arguments first and last. A more precise wording would be seq(from,up/down-to,by).

Generalization

The expressivity power of seq/seqN could be greatly improved by the following additions:

  1. Just like Eigen's ArrayXi::LinSpaced inherits ArrayBase, it would be nice to have seq/seqN inherits ArrayBase (at least conceptually because here we have to deal with symbolic computations). This way the returned object of seq/seqN can be manipulated just like an ArrayXi.
  2. We could add one more overload to seqN with implicit start=0: seqN(n)=[0,1,2,...,n-1]. Yes, we are defaulting the first parameter of the function seqN(start,size) but I don't see how this could be miss-interpreted.

Combining these two minor additions we can reproduce the APL/iota based syntax, but without having to introduce new concepts and with a gradual transition from the pure seq/seqN syntax to the APL one and benefits from the available ArrayXi features. Some random stupid examples just to demonstrate what could be written:

A( a + d*seqN<n> )
A( (a + d*seqN<n>).reverse()+1 )
A( (seq(a,b,d).reverse()*3+7).head(5) )
A( (seqN(a,fix<n>,d).replicate<4>() )     // ! not n arithmetic sequence anymore because of the modulo

Actually, we can even go further by recursively using the current API we are designing:

A( seq(a,b,d)(seqN(last,n,-1)) ) == A( seq(a,b,d).tail(n).reverse() )

Note that here the keyword last refers to 'seq(a,b,d)', not 'A'.

Also, when indexing 1D arrays or vector, this is similar to calling operator() and other methods to the vector itself:

A( seq(a,b,d).reverse().head(5) )  <-->  A( seq(a,b,d) ).reverse().head(5)

This approach also converges towards some aspects of the filtering approach but with a more classic syntax and more consistency/reuse of the current Eigen's API.

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 two placeholders representing either :

  • any integer enumerated in increasing order: anyZ
  • any natural number (0 included) enumerated in increasing order: anyN

In some sense anyZ and anyN represent some unbounded (or half bounded) arithmetic progressions. When passed to a vector or matrix using operator(), then the unbounded extremities are implicitly bounded by the size of the underlying vector or matrix. Some first examples with implicit bounds:

A(anyZ)     ==  A(anyN)      == A(seq(0,last)         == A(all) 
A(-anyZ)    ==  A(last-anyN) == A(seq(last,0,fix<-1>) == A.reverse()
A(2*anyZ)   ==  A(2*anyN)    == A(seq(0,last,2))      // all even entries
A(2*anyZ+1) ==  A(2*anyN+1)  == A(seq(1,last,2)       // all odd entries

Then, one can add explicit inclusive or exclusive bounds using standard comparison operators. The general syntax:

a <= phase+stride*anyZ <= b

represents all integers including the integer 'phase', with increment stride that are in the range [a,b]. This syntax is especially useful for cases for which the bounds a and b and phase are not correlated. A typical example is taking all odd numbers in a given range:

A(a<= 2*anyZ+1<=b)  == A(seq(a+(1-a%2),b, 2)        // all odd entries in [a,b]
A(a<=-2*anyZ+1<=b)  == A(seq(b+(b%2-1),a,-2)        // all odd entries in [a,b] in reverse order

If you want to start at 'a' or finish at 'b' exactly then you can simply pick one of the extremities as the phase:

1:   a <= a+stride*anyZ <= b   // [a, a+stride, a+2stride, ...]
2:   a <= b+stride*anyZ <= b   // [..., b-2stride, b-stride, b]

The first case corresponds to the more classic seq(a,b,stride) that you can more simply mimic with:

(a+stride*anyN)<=b

It might be surprising that with this syntax the bounds a and b does not plays the same role, whereas the syntax seq(a,b,stride) seems to imply symmetry. Actually, this is because, in general, in seq(a,b,stride) the last element of this sequence is not the value b, but a+((b-a)/stride)*stride (unless the upper-bound 'b' has already been carefully computed this way). To reverse this sequence using seq, one might be tempted to write seq(b,a,-stride), which is wrong in general for the same reason. To be safe, one as to write: seq('a+((b-a)/stride)*stride,a,-stride). With the anyN syntax, this is more consistant:

reverse((a+stride*anyN)<=b)  == seq('a+((b-a)/stride)*stride,a,-stride)
     ==((a-stride*anyN)<=b)

Finally, another common use case is taking the last elements from a given index a:

A(a+stride*anyN)        == A(seq(a,last,stride));
A(a-stride*anyZ>=a)     == A(seq(a+((last-a)/stride)*stride,a,stride));   // reverse order

Because the anyZ syntax permits to seamlessly fix the first or last element of the sequence, the reverse order case is much simpler to write.

Now, what if we want the elements from the last one down-to index a:

A(last-stride*anyZ>=a)  == A(seq(last,a,-stride));
A(last+stride*anyZ>=a)  == A(seq(last-((last-a)/stride)*stride,a,-stride));   // reverse order

A bounded arithmetic sequence can still be manipulated using arithmetic operators, just like the half-bounded anyN:

A( (2*(5<anyZ)+1) )       // a random, not well thought example with one bound
A( (2*(5<anyZ<14)+1) )    // a random, not well thought example with two bounds

At this stage, one could simply define anyN as the following alias: auto anyN = (anyZ>=0);. If so, this means that one could put the bounding comparison operators at any stage of the expression:

3:  A( (2*(anyZ>5)+1))<=23 )    // a random, not well thought example!

Of course, it would never be allowed to add a bound to an extremity that is already bounded:

A(anyN>2)                 // compilation error
A(anyN<2)                 // OK
A((13-((2*anyZ+1)<13))>4) // compilation error

Nonetheless, the example #3 above does not read nicely, but after all, it's up to the user to use it wisely?

Finally, if a bounded extremity is too large regarding the size of the underlying matrix/vector, then an out-of-range assertion will be triggered:

VectorXd v(10);
v(anyN<13);      // assertion, the upper-bound is explicitly 12 
v(anyN-1);       // assertion, the lower-bound is explicitly -1
v(-anyN-1);      // assertion, the upper-bound is explicitly -1

Variant D: first attempt to design an expression-based API (kind of broken)

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

Using comparison operator to indicate 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>.

YZ: Regarding expression A(2*-Iota), This is fishy, Iota means (0, 1, ...) but -Iota means (last, last-1, last-2, ...)?
GG: Not at all, Iota means (... -2 -1 0 1 2 ...) and -Iota means (... 2 1 0 -1 -2 ...) --> we should really rename Iota to something meaning any i in N, (N=set of natural integers)
YZ: Say we call it AnyInt. But then A(3 + AnyInt * 4) becomes A(1), A(4), A(7), .... That's a bit unintuitive consider A(3 + Iota(N) *4) starts with A(3). Moreover, to properly implement seq(a, b, incr) (incr >0), you need to write a <= a + AnyInt * incr <= b?

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.

The problem with this notation is that if we view the result of Iota and iota(n) as a pseudo vector, then the meaning of the expression iota < 10 should not be different from v < 10 where v is a regular Eigen vector. If we want to allow boolean mask indexing of the form v(v >= 3 && v < 5), then the type of the expression v < 10 is clearly a boolean mask, not another vector formed by elements of v that are less than 10.

Using filter and predicate to indicate lower-upper-bounds

Here is an alternative syntax for specifying bounds. Here iota(n) follows the regular definition of 0, 1, ... n-1. iota without parameters is the infinite sequence 0, 1, 2, ...

A predicate P is a boolean expression that can be evaluated on a number.

An ordered sequence S, such as a regular Eigen vector or an iota sequence, can be filtered by a predicate P, written as S | P = (s | s in S and P(s)). (...) denotes an ordered sequence here. |, which we call the "such that" operator, is chosen since it is used in set notation { x | predicate on x }.

To form simple predicate, one can use a placeholder like _1 (alias from std::placeholder, or introduce a new one, or call it _x maybe). The expression _1 < 10 returns a predicate that returns true for anything number less than 10. One can also give names to commonly used predicates such as

ubound(u) <=> leq(u) <=> _1 < u
ubound<U> <=> leq<U> <=> _1 < fix<U>
bound(l, u) <=> l <= _1 && _1 <= u
bound<L, U> <=> fix<L> <= _1 && _1 <= fix<U>

Then the bounded arithmetic progression seq(a, b, s) can simply be written as

// Assume s > 0:
seq(a, b, s) = a + iota * s | _1 <= b  
              = a + iota * s | leq(b)
              = a + iota * s | ubound(b)

// When s can be both positive or negative
seq(a, b, s) = a + iota * s | sgn(s) * _1 <= sgn(s) * b 

Moreover the such-that operator | (Sequence v, Predicate p) is properly defined for any sequence and predicate. In the most general case, p is an actual function, so v | p needs to evaluate p on every element of v to create the new sequence. However, with proper overloads, we can recognize expressions such as "10 <= _1 && _1 < n" as bounding predicate. Any arithmetic progression filtered by a bounding predicate will simply be another arithmetic progression.

In addition to conceptual clarity, this syntax allows composition of additional predicates, such that

   A(2 + iota * s | _1 <= upper | _1 % 3 == 0)

Of course, when the predicate is no longer a simple bounding predicate, an actual new sequence needs to be created.

Variant E: Iota, Pipe, and Sequence-wide function

Part of the reason that we are considering implementing bounds either using either relationship operators (<, >, <=, >=) or predicates is that we want to easily represent a bounded arithmetic progression seq(first, last, incr). However, at the end of the day, neither allows a completely generic seq to be written trivially, since the bounding expression needs to take into account the sign of incr:

seq(a, b, d) = a <= a + anyZ * d <= b  if d > 0
seq(a, b, d) = a >= a + anyZ * d >= b  if d < 0

seq(a, b, d) = a + iota * d | ubound(b) if d > 0
seq(a, b, d) = a + iota * d | lbound(b) if d < 0

Another possible convention (e.g., range v3) is that S | F means invoking a sequence-wide function F on the entire sequence of S, relying on the unix pipe connotation.

As before let iota(n) be the finite arithmetic progression 0, 1, 2, ..., n-1. Let iota be the infinite arithmetic progression 0, 1, 2, ... . Then a + iota * s defines an infinite AP. To apply a general predicate like the predicate based expression, we write

a + iota * s | filter(P)

So it's not as simple as a + iota * s | P, But the benefit is that more general purpose sequence-wide function can be used, not just predicates. Moreover, the definition of previously named predicate like ubound and lbound can be changed to work on the whole sequence, so in practice one can still write

a + iota * s | ubound(b)

Now for any finite or infinite AP, we can define a sequence wide function until(bound) that truncate a finite or infinite AP at bound, so

a + iota * s | until(b) = seq(a, b, s)

It automatically takes into consideration the direction of the progression (increasing, or decreasing), so no special handling based on the sign of incr is needed, in contrast to the case of relationship operator and predicate based notations. Also one can have an exclusive version with a different name, e.g., untilX(b).

One can think of additional sequence wide operators, and also form longer expressions, e.g.

(1 + iota * 3 | until(10) | every(2)) + 4 | reverse
= (1 + [0, 3, 6, ...] | until(10) | every(2)) + 4 | reverse 
= ([1, 4, 7, 10] | every(2)) + 4 | reverse
= [1, 7] + 4 | reverse
= [5, 12] | reverse
= [12, 5]

This use of the pipe operator is completely in line with range v3, which might one day becomes STL2.

Interestingly, we are almost revisiting the "named parameter" API here:

  • Iota : Count
  • "A+" : Start, default to 0 if omitted
  • "*S" : Incr, default to 1 if omitted
  • "until(b)" : Last

A few key advantages:

  • Optional argument with default values can just be omitted.
  • No confusion on ordering of parameter.
  • At least for "Count" and "Last", we can omit fixed<N>, and just write iota<Count> and until<Last>.
  • The whole expression has a well-defined meaning. It can be seen as a special case of indexing with a concrete vector. It can be composed with other operations on vectors.

We have some composition rules that preserve AP property:

AP | until(N) = AP   // truncate 
AP | every(k) = AP   // takes very k element
AP | reverse = AP
AP * int = AP
AP + int = AP
AP | ubound(U) = AP  // upper bound
AP | lbound(L) = AP  // lower bound
AP | first(N) = AP   // take first N elements

More generic sequence operations, such as filter, transform, ..., can only be applied to finite AP, which gets turned into into a concrete vector, i.e.,

finite AP | transform(int -> T) = vec<T>
finite AP | filter(int -> boolean) = vec<int>

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
fix<N> best choice so far
ic<N> A bit short and unclear the first time you see it, but sounds pretty good.
c<N> Even shorter! Clearly too polluting.
int_constant<N> Too verbose!
int_<N> Used by boost::MPL
integral_c<N> Used by boost::MPL
int_c<N> Used by boost::hana
Fixed<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
Arithmetic progression/sequence For variant A only
Lower-upper Lower-size
aseq(i,j) aseqn(i,n)
aseqN(i,n)
suffix -n for the size-based version, sounds pretty good.
seq(i,j) seqn(i,n)
seqN(i,n)
without the 'a' prefix, pretty good too, though sequence alone is a bit too general
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.
fap(i,j) fapN(i,n) finite arithmetic progression
range(i,j) rangeN(i,n)
span(i,j) spanN(i,n)
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
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

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.