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

From Eigen
Jump to: navigation, search
(Discussions: removed duplicate APL-inspired API)
(Generalization)
 
(110 intermediate revisions by 3 users not shown)
Line 1: Line 1:
 
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 [http://eigen.tuxfamily.org/bz/show_bug.cgi?id=329 bug report].
 
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 [http://eigen.tuxfamily.org/bz/show_bug.cgi?id=329 bug report].
  
== Needs ==
+
= Needs =
 
We aim to support various indexing mechanisms. For each dimension, we would like to be able to have any of:
 
We aim to support various indexing mechanisms. For each dimension, we would like to be able to have any of:
 
* '''all''': as in <code>A(:)</code>
 
* '''all''': as in <code>A(:)</code>
Line 14: Line 14:
 
We should also be able to provide compile-time information such as ''lengths'' and maybe ''strides''.
 
We should also be able to provide compile-time information such as ''lengths'' and maybe ''strides''.
  
==Some references==
+
=Some references=
  
 +
* Some overviews what other programming languages do: [https://en.wikipedia.org/wiki/Array_slicing] [https://en.wikipedia.org/wiki/Comparison_of_programming_languages_(array)#Slicing]
 
* This [http://ericniebler.com/2014/12/07/a-slice-of-python-in-c post] about the Range v3 library.
 
* This [http://ericniebler.com/2014/12/07/a-slice-of-python-in-c post] about the Range v3 library.
 
* This [http://ericniebler.com/2014/12/07/a-slice-of-python-in-c/#comment-130118 comment] describes an interesting API for range and slices inspired by [https://en.wikipedia.org/wiki/APL_(programming_language) APL].
 
* This [http://ericniebler.com/2014/12/07/a-slice-of-python-in-c/#comment-130118 comment] describes an interesting API for range and slices inspired by [https://en.wikipedia.org/wiki/APL_(programming_language) APL].
 
* [https://libigl.github.io/libigl/matlab-to-eigen.html libigl] extends Eigen by a colon and a slice function
 
* [https://libigl.github.io/libigl/matlab-to-eigen.html libigl] extends Eigen by a colon and a slice function
* Several questions/feature request on StackOverflow: [http://stackoverflow.com/questions/13540147/submatrices-and-indices-using-eigen] [http://stackoverflow.com/questions/16253415/eigen-boolean-array-slicing] [http://stackoverflow.com/questions/40074738/eigen-extracting-submatrix-from-vector-of-indices]
+
* Several questions/feature request on StackOverflow: [http://stackoverflow.com/questions/13540147/submatrices-and-indices-using-eigen] [http://stackoverflow.com/questions/16253415/eigen-boolean-array-slicing] [http://stackoverflow.com/questions/26267194/how-to-extract-a-subvector-of-a-eigenvector-from-a-vector-of-indices-in-eige] [http://stackoverflow.com/questions/40074738/eigen-extracting-submatrix-from-vector-of-indices]
 +
* [http://stackoverflow.com/questions/11364533/why-are-slice-and-range-upper-bound-exclusive Why are slice and range upper bound exclusive in Python?]
  
==C++11 API==
+
=API playground=
Since achieving a lean API in C++03 seems to be rather impossible, let's first see what could be done in C++11 and provide more verbose fallbacks for C++03.
+
 
+
Better than a long story, here is a showcase demo.
+
  
 +
To get some support for the discussions, let's start with a demo of a work-in-progress demo.
  
 
===Source code===
 
===Source code===
Line 33: Line 33:
 
<div class="mw-collapsible-content">
 
<div class="mw-collapsible-content">
 
<source lang="cpp">
 
<source lang="cpp">
 +
 +
#if __cplusplus >= 201103L || defined(_MSC_VER)
 +
#define USE_CXX11
 +
#endif
  
 
#include <iostream>
 
#include <iostream>
 +
#ifdef USE_CXX11
 
#include <array>
 
#include <array>
 +
#endif
 
#include <valarray>
 
#include <valarray>
 
#include <vector>
 
#include <vector>
 +
#include <cassert>
 +
 +
#define USE_EIGEN
 +
 +
#ifdef USE_EIGEN
 
#include <Eigen/Dense>
 
#include <Eigen/Dense>
 
using namespace Eigen;
 
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;
 
using namespace std;
  
 
struct all_t { all_t() {} };
 
struct all_t { all_t() {} };
static const all_t all{};
+
static const all_t all;
  
 
struct shifted_last {
 
struct shifted_last {
 +
  explicit shifted_last(int o) : offset(o) {}
 
   int offset;
 
   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 {
 
struct last_t {
   shifted_last operator- (int offset) const { return shifted_last{offset}; }
+
   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{};
+
static const last_t last;
  
struct Range {
+
 
   Range(Index f, Index l) : m_first(f), m_last(l) {}
+
struct shifted_end {
   Index m_first, m_last;
+
   explicit shifted_end(int o) : offset(o) {}
   Index size() const { return m_last-m_first+1; }
+
   int offset;
   Index operator[] (Index k) { return m_first + k; }
+
   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 RangeLast {
+
struct end_t {
   RangeLast(Index f, last_t) : m_first(f) {}
+
   shifted_end operator- (int offset) const { return shifted_end (-offset); }
   Index m_first;
+
  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;
  
struct RangeShiftedLast {
+
template<int N> struct Index_c {
   RangeShiftedLast(Index f, shifted_last l) : m_first(f), m_last(l) {}
+
   static const int value = N;
   Index m_first;
+
  operator int() const { return value; }
   shifted_last m_last;
+
  Index_c (Index_c<N> (*)() ) {}
 +
   Index_c() {}
 +
   // Needed in C++14 to allow c<XX>():
 +
  Index_c operator() () const { return *this; }
 
};
 
};
  
struct Slice {
+
//--------------------------------------------------------------------------------
  Slice(Index f, Index s, Index l) : m_first(f), m_step(s), m_last(l) {}
+
// Range(first,last) and Slice(first,step,last)
   Index m_first, m_step, m_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 size() const { return (m_last-m_first+m_step)/m_step; }
   Index operator[] (Index k) { return m_first + k*m_step; }
+
   Index operator[] (Index k) const { return m_first + k*m_step; }
 
};
 
};
struct SliceLast {
+
 
  SliceLast(Index f, Index s, last_t) : m_first(f), m_step(s) {}
+
template<typename T> struct cleanup_slice_type { typedef Index type; };
   Index m_first, m_step;
+
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 };
 
};
 
};
  
struct SliceShiftedLast {
+
template<int N,int Default> struct get_compile_time<Index_c<N>,Default> {
   SliceShiftedLast(Index f, Index s, shifted_last l) : m_first(f), m_step(s), m_last(l) {}
+
   enum { value = N };
  Index m_first, m_step;
+
  shifted_last m_last;
+
 
};
 
};
  
Range            span(Index f, Index l)                  { return Range(f,l); }
+
template<typename T> struct is_compile_time        { enum { value = false }; };
RangeLast        span(Index f, last_t l)                { return RangeLast(f,l); }
+
template<int N> struct is_compile_time<Index_c<N> > { enum { value = true }; };
RangeShiftedLast  span(Index f, shifted_last l)          { return RangeShiftedLast(f,l); }
+
Slice            span(Index f, Index s, Index l)        { return Slice(f,s,l); }
+
SliceLast        span(Index f, Index s, last_t l)        { return SliceLast(f,s,l); }
+
SliceShiftedLast  span(Index f, Index s, shifted_last l)  { return SliceShiftedLast(f,s,l); }
+
  
template<int Size=Dynamic,int Step=1,int Start=0>
+
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 {
 
struct iota_t {
   iota_t(Index size=Size) : m_size(size), m_step(Step), m_start(Start) {}
+
   iota_t(Index size=Size) : m_size(size), m_step(Step) {}
   iota_t(Index size,Index step,Index start) : m_size(size), m_step(step), m_start(start) {}
+
   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,Size> m_size;
 
   internal::variable_if_dynamic<Index,Step> m_step;
 
   internal::variable_if_dynamic<Index,Step> m_step;
   internal::variable_if_dynamic<Index,Start> m_start;
+
   StartType m_start;
  
   operator Slice() const {
+
   typedef typename iota2slice<StartType>::type SliceType;
    return Slice(m_start.value(), m_step.value(),
+
  typedef typename iota2slice<StartType>::last_type last_type;
                m_start.value()+(m_size.value()-1)*m_step.value()); }
+
  
   friend iota_t<Size,Dynamic,Start> operator*(int stride,const iota_t &x) {
+
   SliceType range() const {
     return iota_t<Size,Dynamic,Start>{x.m_size.value(),x.m_step.value()*stride,x.m_start.value()};
+
     return SliceType(m_start,
 +
                    m_start+(m_size.value()-1)*m_step.value(),
 +
                    m_step.value());
 
   }
 
   }
  
   iota_t<Size,Step==Dynamic?Dynamic:-Step,Start> operator-() {
+
   operator SliceType() const { return range(); }
     return iota_t<Size,Step==Dynamic?Dynamic:-Step,Start>{m_size.value(),-m_step.value(),m_start.value()};
+
 
 +
  // 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);
 
   }
 
   }
  
   iota_t<Size,Step,Dynamic> operator+(int offset) const {
+
  // negate the step
     return iota_t<Size,Step,Dynamic>{m_size.value(),m_step.value(),m_start.value()+offset};
+
   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);
 
   }
 
   }
  
   friend iota_t<Size,Step,Dynamic> operator+(int offset,const iota_t &x) {
+
   // set first element
     return iota_t<Size,Step,Dynamic>{x.m_size.value(),x.m_step.value(),x.m_start.value()+offset};
+
  iota_t<Size,Step,Index> operator+(int offset) const {
 +
     return iota_t<Size,Step,Index>(m_size.value(),m_step.value(),m_start+offset);
 
   }
 
   }
  
   friend iota_t<Size,Step==Dynamic?Dynamic:-Step,Dynamic> operator-(int offset,const iota_t &x) {
+
  // set first element
     return iota_t<Size,Step==Dynamic?Dynamic:-Step,Dynamic>{x.m_size.value(),-x.m_step.value(),x.m_start.value()+offset};
+
  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));
 
   }
 
   }
 
};
 
};
Line 132: Line 281:
 
template<int Size>
 
template<int Size>
 
iota_t<Size> iota() { return iota_t<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
 
// NOTE that declaring value_type here is needed to make SFINAE works properly in case T does not exhibit value_type
Line 148: Line 317:
 
};
 
};
  
template<typename T> struct is_index { enum { value = std::is_integral<T>::value }; };
+
template<typename T> struct is_index { enum { value = internal::is_integral<T>::value }; };
 
template<> struct is_index<bool> { enum { value = false }; };
 
template<> struct is_index<bool> { enum { value = false }; };
  
Line 157: Line 326:
 
   };
 
   };
 
};
 
};
 +
 +
 +
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 {
 
struct MyMat {
  
   void operator()(Index k) { cout << "singleton: "; cout << k << '\n'; }
+
   MyMat(Index n) : m_size(n) {}
  
 +
  void operator()(Index k) { cout << "singleton: "; cout << k << '\n'; }
 
   void operator()(const char*) { cout << "all\n"; }
 
   void operator()(const char*) { cout << "all\n"; }
 
 
   void operator()(all_t) { cout << "all\n"; }
 
   void operator()(all_t) { cout << "all\n"; }
  
   void operator()(Range ids) { cout << "range-based: "; print_indices(ids);  }
+
  template<typename FirstType,typename LastType,typename StepType>
  void operator()(RangeLast ids) { return operator()({ids.m_first,m_size-1}); }
+
   void operator()(Range_t<FirstType,LastType,StepType> ids) {
  void operator()(RangeShiftedLast ids) { return operator()({ids.m_first,m_size-ids.m_last.offset-1}); }
+
    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";
 +
  }
  
   void operator()(Slice ids) { cout << "slice-based: "; print_indices(ids);  }
+
  template<typename FirstType,typename SizeType,typename StepType>
  void operator()(SliceLast ids) { return operator()(Slice{ids.m_first,ids.m_step,m_size-1}); }
+
   void operator()(Span_t<FirstType,SizeType,StepType> ids) {
   void operator()(SliceShiftedLast ids) { return operator()(Slice{ids.m_first,ids.m_step,m_size-ids.m_last.offset-1}); }
+
    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>
 
   template<typename Indices>
 
   typename internal::enable_if<has_index_value_type<Indices>::value>::type
 
   typename internal::enable_if<has_index_value_type<Indices>::value>::type
   operator()(const Indices &ids) { cout << "index-based: "; print_indices(ids); }
+
   operator()(const Indices &ids) { cout << "index-based: "; print_indices(ids); cout << "\n"; }
  
 
   template<typename Mask>
 
   template<typename Mask>
Line 182: Line 383:
 
   operator()(const Mask &mask) {
 
   operator()(const Mask &mask) {
 
     cout << "mask-based: ";
 
     cout << "mask-based: ";
     for(int k=0; k<mask.size(); ++k)
+
     for(Index k=0; k<Index(mask.size()); ++k)
 
       if(mask[k]) cout << k << " ";
 
       if(mask[k]) cout << k << " ";
 
     cout << "\n";
 
     cout << "\n";
 
   }
 
   }
  
   template<typename T>
+
   template<typename T,Index N>
   void F(T) {}
+
   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:
 
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>
 
   template<typename T>
   void print_indices(T& ids) {
+
   void print_indices(const T& ids) {
     for(int k=0; k<ids.size(); ++k)
+
     for(Index k=0; k<Index(ids.size()); ++k)
 
       cout << ids[k] << " ";
 
       cout << ids[k] << " ";
     cout << "\n";
+
 
 +
     int size = get_size<T>::value;
 +
    if(size!=-1)  cout << "[" << size << "]";
 +
    else          cout << "[]";
 
   }
 
   }
  
   Index m_size = 11;
+
   Index m_size;
 
};
 
};
  
Line 205: Line 428:
 
int main()
 
int main()
 
{
 
{
   ArrayXd eia(10); eia.setRandom();
+
  int n = 13;
   ArrayXi eii(4); eii << 3, 1, 7, 5;
+
#ifdef USE_EIGEN
   valarray<double> vala(10); Map<ArrayXd>(&vala[0],10) = eia;
+
   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;
 
   valarray<int> vali(4); Map<ArrayXi>(&vali[0],4) = eii;
 
   vector<int> veci(4); Map<ArrayXi>(veci.data(),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;
+
   MyMat mat(n);
   cout << "\nc++98 API:\n"
+
   cout << "\nAll, indices, masking:\n"
       <<  "-------------\n";
+
       <<  "----------------------\n";
 
   PRINT( mat(5),"                    " );  // singleton
 
   PRINT( mat(5),"                    " );  // singleton
 
   PRINT( mat(all),"                  " );  // all
 
   PRINT( mat(all),"                  " );  // all
 
   PRINT( mat(""),"                    " );  // all
 
   PRINT( mat(""),"                    " );  // all
  PRINT( mat(span(3,9)),"            " );  // range
+
 
  PRINT( mat(span(3,last-2)),"        " );  // range
+
  PRINT( mat(span(9,3)),"            " );  // empty range (as in MatLab), need to use a slice:
+
  PRINT( mat(span(9,-1,3)),"          " );  // slice
+
  PRINT( mat(span(3,3,last-3)),"      " );  // slice
+
  PRINT( mat(eii),"                  " );  // index-based
+
 
   PRINT( mat(vali),"                  " );  // index-based
 
   PRINT( mat(vali),"                  " );  // index-based
 
   PRINT( mat(veci),"                  " );  // 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
 
   PRINT( mat(eia>0),"                " );  // mask
 +
#endif
 
   PRINT( mat(vala>0.0),"              " );  // mask
 
   PRINT( mat(vala>0.0),"              " );  // mask
  
  cout << "\nc++11 shortcuts:\n"
 
      <<  "-------------\n";
 
  PRINT( mat({3,9}),"                " );  // range
 
  PRINT( mat({3,last}),"              " );  // range
 
  PRINT( mat({3,last-2}),"            " );  // range
 
  PRINT( mat({9,3}),"                " );  // empty range (as in MatLab), need to use a slice:
 
  PRINT( mat({9,-1,3}),"              " );  // slice
 
  PRINT( mat({9,-2,1}),"              " );  // slice
 
  PRINT( mat({3,2,9}),"              " );  // slice
 
  PRINT( mat({3,2,last}),"            " );  // slice
 
  PRINT( mat({3,3,last-3}),"          " );  // slice
 
  PRINT( mat(valarray<int>{5,2,5,6}),"" );  // 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
 
  
   cout << "\nAPL-like API (c++98 compatible):\n"
+
   cout << "\nLower-upper-bound API using range:\n"
       <<  "--------------------------------\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(iota(3)),"              " );
 
   PRINT( mat(2+iota(3)),"            " );
 
   PRINT( mat(2+iota(3)),"            " );
Line 252: Line 510:
 
   PRINT( mat(9-2*iota(3)),"          " );
 
   PRINT( mat(9-2*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
 
}
 
}
 +
  
 
</source>
 
</source>
Line 260: Line 537:
  
 
===Generated output===
 
===Generated output===
 +
 +
The values on the right indicate the extracted compile time size and step as: <tt>[SIZE]{STEP}</tt>. Here, <tt>[]</tt> means a runtime size, and <tt>{}</tt> a runtime step.
  
 
<pre>
 
<pre>
c++98 API:
+
All, indices, masking:
----------
+
----------------------
 
mat(5)                      -> singleton: 5
 
mat(5)                      -> singleton: 5
 
mat(all)                    -> all
 
mat(all)                    -> all
 
mat("")                      -> all
 
mat("")                      -> all
mat(span(3,9))              -> range-based: 3 4 5 6 7 8 9
+
mat(vali)                   -> index-based: 3 1 6 5 []
mat(span(3,last-2))          -> range-based: 3 4 5 6 7 8
+
mat(veci)                   -> index-based: 3 1 6 5 []
mat(span(9,3))              -> range-based:
+
mat(eii)                    -> index-based: 3 1 6 5 [4]
mat(span(9,-1,3))            -> slice-based: 9 8 7 6 5 4 3
+
mat(eii*2-1)                 -> index-based: 5 1 11 9 [4]
mat(span(3,3,last-3))        -> slice-based: 3 6
+
mat(eia>0)                  -> mask-based: 2 4 7 8 9 11 12
mat(eii)                    -> index-based: 3 1 7 5
+
mat(vala>0.0)                -> mask-based: 2 4 7 8 9 11 12
mat(vali)                    -> index-based: 3 1 7 5
+
mat(veci)                   -> index-based: 3 1 7 5
+
mat(eia>0)                  -> mask-based: 2 4 7 8 9
+
mat(vala>0.0)                -> mask-based: 2 4 7 8 9
+
  
c++11 shortcuts:
+
Lower-upper-bound API using range:
----------------
+
----------------------------------
mat({3,9})                   -> range-based: 3 4 5 6 7 8 9
+
mat(range(3,9))              -> strided-range: 3 4 5 6 7 8 9 []{1}
mat({3,last})               -> range-based: 3 4 5 6 7 8 9 10
+
mat(range(3,last))          -> strided-range: 3 4 5 6 7 8 9 10 11 12 []{1}
mat({3,last-2})             -> range-based: 3 4 5 6 7 8
+
mat(range(3,last-2))        -> strided-range: 3 4 5 6 7 8 9 10 []{1}
mat({9,3})                   -> range-based:
+
mat(range(9,3))              -> strided-range: []{1}
mat({9,-1,3})               -> slice-based: 9 8 7 6 5 4 3
+
mat(range(9,3,-1))          -> strided-range: 9 8 7 6 5 4 3 []{}
mat({9,-2,1})               -> slice-based: 9 7 5 3 1
+
mat(range(9,1,-2))          -> strided-range: 9 7 5 3 1 []{}
mat({3,2,9})                 -> slice-based: 3 5 7 9
+
mat(range(last,3,-2))        -> strided-range: 12 10 8 6 4 []{}
mat({3,2,last})             -> slice-based: 3 5 7 9
+
mat(range(last-1,3,-2))      -> strided-range: 11 9 7 5 3 []{}
mat({3,3,last-3})           -> slice-based: 3 6
+
mat(range(3,last-3,3))      -> strided-range: 3 6 9 []{}
mat(valarray<int>{5,2,5,6}) -> index-based: 5 2 5 6
+
mat(range(last-8,last-1,2))  -> strided-range: 4 6 8 10 []{}
mat(array<int,4>{5,2,5,6})   -> index-based: 5 2 5 6
+
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( vector<bool>{false,true,true,false})    -> mask-based: 1 2
 
mat(!valarray<bool>{false,true,true,false})  -> mask-based: 0 3
 
mat(!valarray<bool>{false,true,true,false})  -> mask-based: 0 3
 +
</pre>
  
APL-like API (c++98 compatible):
+
=Discussions=
--------------------------------
+
 
mat(iota(3))                -> slice-based: 0 1 2
+
==nD-Arrays==
mat(2+iota(3))              -> slice-based: 2 3 4
+
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(3+2*iota(3))            -> slice-based: 3 5 7
+
<pre>mat(span(0,2,last), "", vec>0)</pre>
mat(9-iota(3))               -> slice-based: 9 8 7
+
which is why an API like <code>mat.slice(...)</code> is a no go. This is also a crucial limitation of the current <code>Eigen::Block</code>-based API.
mat(9-2*iota(3))            -> slice-based: 9 7 5
+
 
mat(9-2*iota<3>())           -> slice-based: 9 7 5
+
==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: <code>A=rand(n,n); A(A<0)=0;</code>. 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:
 +
*<tt>start</tt> - the index of the first element in the slice
 +
*<tt>stop</tt> - either the index of the last element in the slice (inclusive), or one more than the index of last element in the slice (exclusive)
 +
*<tt>size</tt> - the length of the slice (= exclusive(stop) - first = inclusive(stop) - first + 1)
 +
*<tt>incr</tt> - the increment to go from one element to the next one (default 1)
 +
 
 +
Among those 5 quantities (because <tt>stop</tt> has two variants), only three are necessary to define a slice. The cases where <tt>incr</tt> 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 <tt>incr</tt> is given either explicitly or implicitly as ti defaults to 1. We can distinguish the following cases:
 +
# Explicit-size:
 +
## [<tt>start</tt>,<tt>size</tt>]
 +
## [<tt>inclusive(stop)</tt>,<tt>size</tt>]
 +
## [<tt>exclusive(stop)</tt>,<tt>size</tt>]
 +
# Lower-upper-bounds:
 +
## Inclusive upper-bound: [<tt>start</tt>,<tt>inclusive(stop)</tt>]
 +
## Exclusive upper-bound: [<tt>start</tt>,<tt>exclusive(stop)</tt>]
 +
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:
 +
# <tt>last</tt> refers to the index of the last element (inclusive)
 +
# <tt>end</tt> 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 <tt>SIZE</tt> and <tt>INCR</tt>.
 +
 
 +
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 <tt>incr</tt>. 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.
 +
 
 +
<!--
 +
Moreover, there is no need to follow the same logic for the two main cases. Indeed, if we go to something like:
 +
<pre>
 +
foo(start,size[,incr]);
 +
bar(start,last[,incr]);
 
</pre>
 
</pre>
 +
then we would have no pick function names that explicitly convey the expected arguments, e.g.:
 +
<pre>
 +
sliceStartLen(start,size[,incr]);
 +
sliceBounds(start,last[,incr]);
 +
</pre>
 +
is not very concise and it becomes cumbersome when <tt>size</tt> is specified at compile-time:<tt>sliceStartLen<size>(start[,incr]);</tt>. If on top of that we want to be able to specify the <tt>incr</tt> at compile-time, then we're definitely screwed.
 +
-->
 +
 +
==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 <tt>stop</tt>.
 +
At a first glance, we can make the following remarks:
 +
* Advantage of <tt>inclusive(stop)</tt>
 +
** closer to matlab
 +
** '''symmetry''' with <tt>first</tt>
 +
** permits to introduce a ''keyword'' <tt>Eigen::last</tt> to refer to the last element of the underlying array whereas the exclusive approach would suggest introducing a <tt>Eigen::end</tt> that conflicts with <tt>std::end</tt>.
 +
** Used in several languages/libraries: MatLab, Mathematica, Fortran, R, Pascal, Ada, Perl, boost::irange, etc.
 +
* Advantage of <tt>exclusive(stop)</tt>
 +
** 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:
 +
 +
{| class="wikitable"
 +
|-
 +
! intents !! seq (inclusive) or seqn !! seqX (exclusive) !! python !! expression-based with comparison operator !! expr with predicate<br/>(when appropriate)
 +
|-
 +
| first k elements || <tt>A(seqn(0,k)) </tt> || <tt>A(seqX(0,k))</tt> || <tt>A[:k]</tt>  || <tt>A(iota(k))</tt> <br/>  <tt>A(Iota<k)</tt> ||
 +
|-
 +
| first k elements (descending)  || <tt>A(seqn(k-1,k,fix<-1>)) </tt> <br/> <tt>A(seq(k-1,0,fix<-1>))</tt> || <tt>A(seqX(k-1,-1,fix<-1>))</tt> || <tt>A[k-1::-1]</tt> ||<tt>A(-Iota<k)</tt><br/><tt>A(k-1-iota(k))</tt>
 +
|-
 +
| first elements to index i || <tt>A(seq(0,i))</tt> || <tt>A(seqX(0,i+1))</tt> || <tt>A[:i+1]</tt>  || <tt>A(iota(i+1))</tt><br/><tt>A(Iota<=i)</tt> ||
 +
|-
 +
| first elements to index i (descending)  || <tt>A(seq(i,0,fix<-1>))</tt>  <br/> <tt> A(seqn(i,i+1,fix<-1>)) </tt> || <tt>A(seqX(i,-1,fix<-1>))</tt> || <tt>A[i::-1]</tt> || <tt>A(i-iota(i+1))</tt><br/><tt>A(-Iota<=i)</tt> <br/> <tt>A(i>=-Iota)</tt>
 +
|-
 +
| last k elements (ascending) || <tt>A(seq(end-k,last))</tt> <br/> <tt>A(seqn(last-k+1,k)) </tt> <br/> <tt>A(seqn(end-k,k)) </tt> || <tt>A(seqX(last-k+1,last+1))</tt><br/><tt>A(seqX(end-k,end))</tt> || <tt>A[-k:]</tt>  || <tt>A(last-k<Iota)</tt> <br/> <tt>A(reverse(last-iota(k)))</tt> ||
 +
|-
 +
|  last k elements <br/> (descending) || <tt>A(seqn(last,k,fix<-1>)) </tt> <br/> <tt>A(aseq(last,end-k,fix<-1>))</tt> || <tt>A(seqX(last,last-k,fix<-1>))</tt> || <tt>A[:-(k+1):-1]</tt>  || <tt>A(last-iota(k))</tt> <br/> <tt>A(last-k<-Iota)</tt> ||
 +
|-
 +
| all but the last k || <tt>A(seq(0,last-k))</tt> <br/> <tt>A(seqn(0,size-k)) </tt> || <tt>A(seqX(0,last-k+1))</tt> || <tt>A[:-k]</tt>  || <tt>A(iota(last+1-k))</tt> <br/> <tt>A(Iota<=last-k)</tt> <br/> <tt>A(Iota<end-k)</tt> || <tt>A(Iota<nowiki>|</nowiki>ubound(last-k))</tt> <br/> <tt>A(Iota<nowiki>|</nowiki>_1<=last-k)</tt>
 +
|-
 +
| all but the last k <br/> (descending) || <tt>A(seq(last-k,0,fix<-1>))</tt> <br/> <tt>A(seqn(last-k,size-k,fix<-1>)) </tt> || <tt>A(seqX(last-k,-1,fix<-1>))</tt> || <tt>A[-(k+1)::-1]</tt>  ||  <tt>A(-Iota<=last-k)</tt> <br/> <tt>A(-Iota<end-k)</tt> <br/> <tt>A(last-k-Iota(last+1-k))</tt> ||
 +
|-
 +
|  subset of the first k <br/> elements with step s || <tt>A(seq(0,k-1,s))</tt> <br/> <tt>A(seqn(0,(k+s-1)/s,s)) </tt> || <tt>A(seqX(0,k,s))</tt> || <tt>A[:k:s]</tt>  || <tt>A(Iota*s<k)</tt> || <tt>A(Iota*s<nowiki>|</nowiki>ubound(k-1))</tt> <br/> <tt>A(Iota*s<nowiki>|</nowiki>_1<k)</tt>
 +
|-
 +
|  subset of the first k <br/> elements with step s <br/> (descending) || <tt>A(seq(k-1,fix<0>,-s))</tt> <br/> <tt>A(seqn(k-1,(k+s-1)/s,-s)) </tt> || <tt>A(seqX(k-1,-1,-s))</tt> || <tt>A(k-1-Iota*s)</tt> <br/> <tt>A[k-1::-s]</tt> || <tt>A(-Iota*s<k)</tt> ||  <tt>A(k-1-Iota*s<nowiki>|</nowiki>lbound<0>)</tt> <br/> <tt>A(k-1-Iota*s<nowiki>|</nowiki>_1>=0)</tt>
 +
|-
 +
|  pick last k elements with <br/> step s (ascending) || <tt>A(seq(last-(k-1)*s,last,s))</tt><br/><tt>A(seqn(last-(k-1)*s,k,s)) </tt> <br/> <tt>A(flip(seqn(last,k,-s))) </tt> || <tt>A(seqX(last-(k-1)*s,end,s))</tt> || <tt>A[-(k-1)*s-1::s]</tt>  || <tt>A(flip(last-iota(k)*s))</tt> <br/> <tt>A(last-(k-1)*s+iota(k)*s)</tt> ||
 +
|-
 +
|  pick last k elements with <br/> step s (descending) || <tt>A(seq(last,end-(k-1)*s,-s))</tt> <br/> <tt>A(seqn(last,k,-s)) </tt> || <tt>A(seqX(last,last-(k-1)*s-1,-s))</tt> || <tt>A[:-(k-1)*s-2:-s]</tt>  || <tt>A(last-iota(k)*s)</tt> ||
 +
|-
 +
|  elements || <tt>A(seq(,))</tt> || <tt>A(seqn(,)) </tt> || <tt>A(seqX(,))</tt> || <tt>A[:]</tt> || <tt>A()</tt>
 +
|-
 +
|  elements || <tt>A(seq(,))</tt> || <tt>A(seqn(,)) </tt> || <tt>A(seqX(,))</tt> || <tt>A[:]</tt> || <tt>A()</tt>
 +
|-
 +
|  elements || <tt>A(seq(,))</tt> || <tt>A(seqn(,)) </tt> || <tt>A(seqX(,))</tt> || <tt>A[:]</tt> || <tt>A()</tt>
 +
|-
 +
|  elements || <tt>A(seq(,))</tt> || <tt>A(seqn(,)) </tt> || <tt>A(seqX(,))</tt> || <tt>A[:]</tt> || <tt>A()</tt>
 +
|-
 +
|}
 +
 +
==Function-based API==
 +
In this section, we will discuss how to design an API of the form:
 +
<pre>
 +
Matrix A;
 +
auto subA = A(foo(...), foo(...));
 +
</pre>
 +
with some proper names for the function <tt>foo</tt> 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:
 +
<pre>
 +
span(start,size[,incr]);
 +
span<SIZE>(start[,incr]);
 +
span<INCR>(start,size);
 +
span<SIZE,INCR>(start);
 +
</pre>
 +
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 <tt>integral_constant<N></tt>, that is currently abbreviated as <tt>ic</tt>.
 +
The previous <tt>span(start,size[,incr]);</tt> example, can then be revisited as follows:
 +
<pre>
 +
span(start, ic<SIZE>)
 +
span(start, ic<SIZE>, incr)
 +
span(start, size, ic<INCR>)
 +
span(start, ic<SIZE>, ic<INCR>)
 +
</pre>
 +
No ambiguity anymore!
 +
Of course, for concision this assumes the user has imported a short <tt>Eigen::ic</tt> in the current scope.
 +
<!-- The long version could be <tt>Eigen::Index_c<size>()</tt>. -->
 +
And for indexing from the end:
 +
<pre>
 +
// 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)
 +
</pre>
 +
 +
====Arithmetic progressions with lower-upper-bounds (range)====
 +
 +
This case is simpler because there is fewer variants. Following the previous ''span'' variants, we coudl propose:
 +
<pre>
 +
range(start, stop)
 +
range(start, stop, incr)
 +
range(start, stop, ic<INCR>)
 +
</pre>
 +
Here, if we agree that only the increment can be specified at compile-time, a template overload like:
 +
<pre>
 +
range<INCR>(start, stop)
 +
</pre>
 +
would also be acceptable.
 +
 +
<div class="mw-collapsible mw-collapsed">
 +
===Variant B: single function name for arithmetic progressions (dismissed)===
 +
<div class="mw-collapsible-content">
 +
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 <tt>aseq</tt> as a contraction for ''arithmetic sequence''.
 +
Then we could have:
 +
<pre>
 +
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>)
 +
</pre>
 +
 +
The main drawbacks is that it implies a weird API for the simple <tt>aseq(start, Size(size))</tt> case, and it is less generalizable as it requires a new identifier in Eigen's namespace for each argument name, whereas <tt>ic<N></tt> can be used everywhere.
 +
</div>
 +
</div>
 +
 +
===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 <tt>seq(a,b,incr)</tt> is '''not''' '<tt>seq(b,a,-incr)</tt> but <tt>seq('a+((b-a)/incr)*incr,a,-incr)</tt>
 +
** As a consequence, we cannot call these arguments ''first'' and ''last''. A more precise wording would be <tt>seq(from,up/down-to,by)</tt>.
 +
 +
===Generalization===
 +
 +
The expressivity power of seq/seqN could be greatly improved by the following additions:
 +
# 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.
 +
# We could add one more overload to seqN with implicit start=0: <tt>seqN(n)=[0,1,2,...,n-1]</tt>. 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 [[#With an explicit size|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:
 +
<pre>
 +
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
 +
</pre>
 +
Actually, we can even go further by recursively using the current API we are designing:
 +
<pre>
 +
A( seq(a,b,d)(seqN(last,n,-1)) ) == A( seq(a,b,d).tail(n).reverse() )
 +
</pre>
 +
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:
 +
<pre>
 +
A( seq(a,b,d).reverse().head(5) )  <-->  A( seq(a,b,d) ).reverse().head(5)
 +
</pre>
 +
 +
This approach also converges towards some aspects of the [[#Variant E: Iota, Pipe, and Sequence-wide function | 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: <tt>anyZ</tt>
 +
* any natural number (0 included) enumerated in increasing order: <tt>anyN</tt>
 +
 +
In some sense <tt>anyZ</tt> and <tt>anyN</tt> represent some unbounded (or half bounded) arithmetic progressions. When passed to a vector or matrix using <tt>operator()</tt>, then the unbounded extremities are implicitly bounded by the size of the underlying vector or matrix. Some first examples with implicit bounds:
 +
<pre>
 +
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
 +
</pre>
 +
 +
Then, one can add explicit inclusive or exclusive bounds using standard comparison operators. The general syntax:
 +
<pre>
 +
a <= phase+stride*anyZ <= b
 +
</pre>
 +
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:
 +
<pre>
 +
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
 +
</pre>
 +
If you want to start at 'a' or finish at 'b' exactly then you can simply pick one of the extremities as the ''phase'':
 +
<pre>
 +
1:  a <= a+stride*anyZ <= b  // [a, a+stride, a+2stride, ...]
 +
2:  a <= b+stride*anyZ <= b  // [..., b-2stride, b-stride, b]
 +
</pre>
 +
The first case corresponds to the more classic <tt>seq(a,b,stride)</tt> that you can more simply mimic with:
 +
<pre>
 +
(a+stride*anyN)<=b
 +
</pre>
 +
It might be surprising that with this syntax the bounds ''a'' and ''b'' does not plays the same role, whereas the syntax <tt>seq(a,b,stride)</tt> seems to imply symmetry.
 +
Actually, this is because, in general, in <tt>seq(a,b,stride)</tt> 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 <tt>seq(b,a,-stride)</tt>, which is wrong in general for the same reason. To be safe, one as to write: <tt>seq('a+((b-a)/stride)*stride,a,-stride)</tt>. With the ''anyN'' syntax, this is more consistant:
 +
<pre>
 +
reverse((a+stride*anyN)<=b)  == seq('a+((b-a)/stride)*stride,a,-stride)
 +
    ==((a-stride*anyN)<=b)
 +
</pre>
 +
Finally, another common use case is taking the last elements from a given index ''a'':
 +
<pre>
 +
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
 +
</pre>
 +
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'':
 +
<pre>
 +
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
 +
</pre>
 +
<!--
 +
<pre>
 +
(0<=anyZ)          == (anyZ>=0)              == anyN
 +
(0<anyZ)            == (anyZ>0)              == (0<anyN)      // aka Z^+
 +
A((a+c*anyN)<=b)    == A(a<=(a+c*anyZ)<=b)    == A(seq(a,b,c))
 +
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
 +
</pre>
 +
-->
 +
A bounded arithmetic sequence can still be manipulated using arithmetic operators, just like the half-bounded <tt>anyN</tt>:
 +
<pre>
 +
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
 +
</pre>
 +
At this stage, one could simply define <tt>anyN</tt> as the following alias: <tt>auto anyN = (anyZ>=0);</tt>.
 +
If so, this means that one could put the bounding comparison operators at any stage of the expression:
 +
<pre>
 +
3:  A( (2*(anyZ>5)+1))<=23 )    // a random, not well thought example!
 +
</pre>
 +
Of course, it would never be allowed to add a bound to an extremity that is already bounded:
 +
<pre>
 +
A(anyN>2)                // compilation error
 +
A(anyN<2)                // OK
 +
A((13-((2*anyZ+1)<13))>4) // compilation error
 +
</pre>
 +
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:
 +
<pre>
 +
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
 +
</pre>
 +
 +
<div class="mw-collapsible mw-collapsed">
 +
==Variant D: first attempt to design an expression-based API (kind of broken)==
 +
 +
<div class="mw-collapsible-content">
 +
 +
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 <tt>iota</tt> function generating an arithmetic progression starting at 0 and with unit increment:
 +
<pre>
 +
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
 +
</pre>
 +
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 <tt>iota<></tt> or <tt>iota()</tt>, 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 <tt>iota(n)/iota<N></tt> versions that generates ''n'' integers from ''0'' to ''n-1''. Then we can apply a scaling factor and an offset to define the <tt>incr</tt> and <tt>start</tt> respectively. The syntax is:
 +
<pre>
 +
[start +] [incr*]iota(size)
 +
</pre>
 +
For compile-time size:
 +
<pre>
 +
[start +] [incr*]iota<SIZE>()    // c++98
 +
[start +] [incr*]iota<SIZE>      // c++14
 +
</pre>
 +
For a compile-time increment, we could use the notion of <tt>integral_constant<N></tt>, that is, as before, currently abbreviated as <tt>ic</tt>:
 +
<pre>
 +
[start +] ic<INCR>()*iota(size)    // c++98
 +
[start +] ic<INCR>*iota(size)      // c++98 with shortcut
 +
</pre>
 +
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:
 +
<pre>
 +
// 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)
 +
</pre>
 +
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 <tt>flip</tt> or <tt>reverse</tt>...
 +
 +
===Using comparison operator to indicate lower-upper-bounds===
 +
Here we start with the default generator <tt>Iota</tt> 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'':
 +
<pre>
 +
// 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)
 +
</pre>
 +
Note that we can define a compile time -1 increment without messing with <tt>ic<-1></tt>.
 +
 +
'''YZ''': Regarding expression <tt>A(2*-Iota)</tt>, This is fishy, Iota means (0, 1, ...) but -Iota means (last, last-1, last-2, ...)? <br/>
 +
'''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) <br/>
 +
'''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 <tt>seq(a, b, incr)</tt> (incr >0), you need to write <tt>a <= a + AnyInt * incr <= b</tt>?
 +
 +
Then explicit bounds (either inclusive or exclusive) can be defined using comparison operators:
 +
<pre>
 +
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))
 +
</pre>
 +
Not too bad, except the 4th line <tt>A(3<=ic<2>*Iota<10)</tt> that is getting messy with all these <tt><</tt> and <tt>></tt> symbols.
 +
 +
The problem with this notation is that if we view the result of <tt>Iota</tt> and <tt>iota(n)</tt> as a pseudo vector, then the meaning of the expression <tt> iota < 10 </tt> should not be different from <tt> v < 10 </tt> where <tt>v</tt> is a regular Eigen vector. If we want to allow boolean mask indexing of the form <tt>v(v >= 3 && v < 5)</tt>, then the type of the expression <tt> v < 10 </tt> is clearly a boolean mask, not another vector formed by elements of <tt>v</tt> that are less than <tt>10</tt>.
 +
 +
=== Using filter and predicate to indicate lower-upper-bounds ===
 +
 +
Here is an alternative syntax for specifying bounds. Here <tt>iota(n)</tt> follows the regular definition of <tt>0, 1, ... n-1</tt>. <tt>iota</tt> without parameters is the infinite sequence <tt>0, 1, 2, ...</tt>
 +
 +
A predicate <tt>P</tt> is a boolean expression that can be evaluated on a number.
 +
 +
An ordered sequence <tt>S</tt>, such as a regular Eigen vector or an <tt>iota</tt> sequence, can be filtered by a predicate <tt>P</tt>, written as <tt>S | P = (s | s in S and P(s))</tt>. <tt>(...)</tt> denotes an ordered sequence here.  <tt>|</tt>, which we call the "such that" operator, is chosen since it is used in set notation <tt>{ x | predicate on x }</tt>.
 +
 +
To form simple predicate, one can use a placeholder like <tt>_1</tt> (alias from std::placeholder, or introduce a new one, or call it <tt>_x</tt> maybe). The expression <tt>_1 < 10</tt> returns a predicate that returns true for anything number less than 10. One can also give names to commonly used predicates such as
 +
<pre>
 +
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>
 +
</pre>
 +
 +
Then the bounded arithmetic progression <tt>seq(a, b, s)</tt> can simply be written as
 +
<pre>
 +
// 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
 +
</pre>
 +
 +
Moreover the such-that operator <tt>| (Sequence v, Predicate p)</tt> is properly defined for any sequence and predicate. In the most general case, <tt>p</tt> is an actual function, so <tt>v | p</tt> needs to evaluate <tt>p</tt> on every element of <tt>v</tt> to create the new sequence. However, with proper overloads, we can recognize expressions such as "10 <= _1 && _1 < n" as <em>bounding predicate</em>. 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
 +
<pre>
 +
  A(2 + iota * s | _1 <= upper | _1 % 3 == 0)
 +
</pre>
 +
Of course, when the predicate is no longer a simple bounding predicate, an actual new sequence needs to be created.
 +
 +
</div>
 +
</div>
 +
 +
== 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 <tt>seq(first, last, incr)</tt>.  However, at the end of the day, neither allows a completely generic <tt>seq</tt> to be written trivially, since the bounding expression needs to take into account the sign of <tt>incr</tt>:
 +
<pre>
 +
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
 +
</pre>
 +
 +
Another possible convention  (e.g., range v3) is that <tt>S | F</tt> means invoking a sequence-wide function <tt>F</tt> on the entire sequence of <tt>S</tt>, relying on the unix pipe connotation.
 +
 +
As before let <tt>iota(n)</tt> be the finite arithmetic progression <tt>0, 1, 2, ..., n-1</tt>. Let <tt>iota</tt> be the infinite arithmetic progression <tt>0, 1, 2, ... </tt>. Then <tt> a + iota * s</tt> defines an infinite AP. To apply a general predicate like the predicate based expression, we write
 +
<pre>
 +
a + iota * s | filter(P)
 +
</pre>
 +
So it's not as simple as <tt>a + iota * s | P</tt>, 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 <tt>ubound</tt> and <tt>lbound</tt> can be changed to work on the whole sequence, so in practice one can still write
 +
<pre>
 +
a + iota * s | ubound(b)
 +
</pre>
 +
 +
Now for any finite or infinite AP, we can define a sequence wide function <tt>until(bound)</tt> that truncate a finite or infinite AP at <tt>bound</tt>, so
 +
<pre>
 +
a + iota * s | until(b) = seq(a, b, s)
 +
</pre>
 +
It automatically takes into consideration the direction of the progression (increasing, or decreasing), so no special handling based on the sign of <tt>incr</tt> 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.
 +
<pre>
 +
(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]
 +
</pre>
 +
 +
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:
 +
 +
<pre>
 +
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
 +
</pre>
 +
 +
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.,
 +
<pre>
 +
finite AP | transform(int -> T) = vec<T>
 +
finite AP | filter(int -> boolean) = vec<int>
 +
</pre>
 +
 +
=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.
  
===Discussions===
+
{| class="wikitable"
 +
|-
 +
!colspan="2"| Integral constant !! For variants A and C
 +
|-
 +
|colspan="2"| <tt>fix<N></tt> || best choice so far
 +
|-
 +
|colspan="2"| <tt>ic<N></tt> || A bit short and unclear the first time you see it, but sounds pretty good.
 +
|-
 +
|colspan="2"| <tt>c<N></tt> || Even shorter! Clearly too polluting.
 +
|-
 +
|colspan="2"| <tt>int_constant<N></tt> || Too verbose!
 +
|-
 +
|colspan="2"| <tt> int_<N></tt> || Used by boost::MPL
 +
|-
 +
|colspan="2"| <tt> integral_c<N></tt> || Used by boost::MPL
 +
|-
 +
|colspan="2"| <tt> int_c<N></tt> || Used by boost::hana
 +
|-
 +
|colspan="2"| <tt>Fixed<N></tt> ||
 +
|-
 +
|colspan="2"| <tt>Cst<N></tt> ||
 +
|-
 +
|colspan="2"| <tt>IndexC<N></tt> || Reference to Eigen::Index
 +
|-
 +
|colspan="2"| <tt>CIndex<N></tt> ||
 +
|-
 +
|colspan="2"| <tt>CstIndex<N></tt> ||
 +
|-
 +
|colspan="2"| <tt>Static<N></tt> || static is a C++ keyword
 +
|-
 +
|colspan="2"| <tt>Stat<N></tt> || Shortcut for Static
 +
|-
 +
|colspan="2"| <tt>(int[N]){}</tt> || A name-free variant, just for fun, does not work with N<1
 +
|-
 +
!colspan="2"| Arithmetic progression/sequence !! For variant A only
 +
|-
 +
! Lower-upper !! Lower-size !!
 +
|-
 +
| <tt>aseq(i,j)</tt> || <tt>aseqn(i,n)</tt> <br/> <tt>aseqN(i,n)</tt> || suffix -n for the size-based version, sounds pretty good.
 +
|-
 +
| <tt>seq(i,j)</tt> || <tt>seqn(i,n)</tt><br/><tt>seqN(i,n)</tt> || without the 'a' prefix, pretty good too, though ''sequence'' alone is a bit too general
 +
|-
 +
| <tt>range(i,j)</tt> || <tt>span(i,n)</tt> ||
 +
|-
 +
| <tt>fromto(i,j)</tt> ||  <tt>fromfor(i,n)</tt> || ''fromto'' looks quite nice, but not ''fromfor''.
 +
|-
 +
| <tt>fromto(i,j)</tt> ||  <tt>fromlen(i,n)</tt> || not really better than ''fromfor''.
 +
|-
 +
| <tt>fap(i,j)</tt> || <tt>fapN(i,n)</tt> ||  ''finite arithmetic progression''
 +
|-
 +
| <tt>range(i,j)</tt> || <tt>rangeN(i,n)</tt> ||
 +
|-
 +
| <tt>span(i,j)</tt> || <tt>spanN(i,n)</tt> ||
 +
|-
 +
!colspan="2"| Full range (aka <tt>:</tt> in Matlab)  !! For variants A and B
 +
|-
 +
|colspan="2"| <tt>All</tt> ||
 +
|-
 +
|colspan="2"| <tt>":"</tt> || catch <tt>char*</tt> at compile-time, with runtime check in debug mode
 +
|-
 +
!colspan="2"| Compile-time 0 !! For variants A, B, C
 +
|-
 +
|colspan="2"| <tt>Zero</tt> ||
 +
|-
 +
|colspan="2"| <tt>"0"</tt> || catch <tt>char*</tt> at compile-time, with runtime check in debug mode
 +
|-
 +
!colspan="2"| Arithmetic progression/sequence !! For variant B only.
 +
|-
 +
|colspan="2"| <tt>aseq(...)</tt> ||
 +
|-
 +
|colspan="2"| <tt>seq(...)</tt> ||
 +
|-
 +
|colspan="2"| <tt>aprog(...)</tt> ||
 +
|-
 +
|colspan="2"| <tt>fap(...)</tt> || finite arithmetic progression
 +
|-
 +
|colspan="2"| <tt>smseq(...)</tt> || strictly monotonic sequence
 +
|-
 +
|colspan="2"| <tt>range(...)</tt> ||
 +
|-
 +
|colspan="2"| <tt>span(...)</tt> ||
 +
|-
 +
!colspan="2"| Size  !! For variant B only
 +
|-
 +
|colspan="2"| <tt>Size<N></tt> || Very clear, but polluting
 +
|-
 +
|colspan="2"| <tt>Length<N></tt> ||
 +
|-
 +
|colspan="2"| <tt>Len<N></tt> ||
 +
|-
 +
|colspan="2"| <tt>Count<N></tt> ||
 +
|-
 +
|colspan="2"| <tt>Card<N></tt> ||
 +
|-
 +
!colspan="2"| Increment  !! For variant B only
 +
|-
 +
|colspan="2"| <tt>Incr<N></tt> || Short and precise
 +
|-
 +
|colspan="2"| <tt>Stride<N></tt> || Already used in Eigen (<tt>Stride<Outer,Inner> <-> Stride<LeadingDim,Incr></tt>)
 +
|-
 +
|colspan="2"| <tt>Step<N></tt> || ''step'' can mean a lot of different things
 +
|-
 +
|}
  
# 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: <pre>mat(span(0,2,last), "", vec>0)</pre> which is why an API like <code>mat.slice(...)</code> is a no go. This is also a crucial limitation of the current <code>Eigen::Block</code>-based API.
+
===Others===
# We use the word 'last' instead of Matlab's 'end' to (1) reduce collision with std::end, and (2) make it more consistent with STL's end as 'last' refer to the last element whereas 'end' refers to the last+1.
+
<!--# We use the word 'last' instead of Matlab's 'end' to (1) reduce collision with std::end, and (2) make it more consistent with STL's end as 'last' refer to the last element whereas 'end' refers to the last+1.-->
 
# Using <code>""</code> for '''all''' might looks odd, but it would avoid introducing a polluting identifier <code>all</code> (nobody wants to write <code>Eigen::all</code>). We have the same issue with <code>last</code> though.
 
# Using <code>""</code> for '''all''' might looks odd, but it would avoid introducing a polluting identifier <code>all</code> (nobody wants to write <code>Eigen::all</code>). We have the same issue with <code>last</code> though.
# c++98 <code>span(first,[step,]last)</code> API:
+
<!--# c++98 <code>range(first,[incr,]last)</code> API:
 
#* Pros:
 
#* Pros:
#** mimics MatLab first:[step:]last API
+
#** mimics MatLab first:[incr:]last API
 
#** enable indexing from the end
 
#** enable indexing from the end
 
#* Cons:
 
#* Cons:
 
#** the argument order is not explicit
 
#** the argument order is not explicit
 
#** we could specify compile time size with <code>span<3>(4,6)</code> or <code>span<3>(4,2,8)</code> but then we have redundant information...
 
#** we could specify compile time size with <code>span<3>(4,6)</code> or <code>span<3>(4,2,8)</code> but then we have redundant information...
 +
-->
 
# [https://en.wikipedia.org/wiki/APL_(programming_language) APL]-like API:
 
# [https://en.wikipedia.org/wiki/APL_(programming_language) APL]-like API:
 
#* Pros:
 
#* Pros:
Line 323: Line 1,239:
 
#* Cons:
 
#* Cons:
 
#** cumbersome?
 
#** cumbersome?
#** needs to specify length instead of first-last
+
#** needs to specify size instead of first-last
# c++11 <code>mat({first,[step,]last})</code> API:
+
<!--# c++11 <code>mat({first,[incr,]last})</code> API:
 
#* Pros:
 
#* Pros:
 
#** very compact
 
#** very compact
#** closely mimics MatLab first:step:last API
+
#** closely mimics MatLab first:incr:last API
 
#** no need to import a possibly polluting name as <code>Eigen::span</code>
 
#** no need to import a possibly polluting name as <code>Eigen::span</code>
 
#** enable indexing from the end
 
#** enable indexing from the end
Line 333: Line 1,249:
 
#** might be confusing with <code>mat({i0,i1,i2})</code> (which is why <code>mat({i0,i1,i2,i3,...})</code> is not allowed)
 
#** might be confusing with <code>mat({i0,i1,i2})</code> (which is why <code>mat({i0,i1,i2,i3,...})</code> is not allowed)
 
#** does not permit to specify compile-time size
 
#** does not permit to specify compile-time size
#** the argument order is not explicit
+
#** the argument order is not explicit-->
# What does not work:
+
<!--# What does not work:
 
<source lang="cpp">
 
<source lang="cpp">
 
mat({5,2,6,7,1});                                // needs to explicitly use a valarray, array, or ...
 
mat({5,2,6,7,1});                                // needs to explicitly use a valarray, array, or ...
 
mat({true,false,false,true,true,false,true});    // needs to explicitly use a valarray, array, or ...
 
mat({true,false,false,true,true,false,true});    // needs to explicitly use a valarray, array, or ...
</source>
+
</source>-->
 
+
==C++03 fallback API==
+
  
==Internal implementation details==
+
=Internal implementation details=
  
===Performance===
+
==Performance==
  
 
One could be tempted to leverage AVX2 gather/scatter instructions, but those are horribly slow. Better emulate them for now.
 
One could be tempted to leverage AVX2 gather/scatter instructions, but those are horribly slow. Better emulate them for now.

Latest revision as of 14:00, 8 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 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.