View source for Working notes - Indexing++

The goal of this page is to summarize the different ideas and working plan to (finally!) provide support for flexible row/column indexing in Eigen.  See this [http://eigen.tuxfamily.org/bz/show_bug.cgi?id=329 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 <code>A(:)</code>
* '''singleton''': as in <code>A(3)</code> (already supported!)
* '''index-array''': as in <code>A([3 1 5])</code>
* '''range-based''': as in <code>A(3:9)</code>, but sometimes it is also convenient to specify the ''start+length'' similar to what <code>Eigen::Block</code> offers.
* '''slicing''': as in <code>A(3:2:9)</code>. Same as for range, we sometime also want to specify ''start+length+stride''
* '''negative indices''': as in <code>A(end-3)</code> or <code>A(3:2:end-1)</code>
* '''boolean masking''': as in <code>A(A>0)</code>

All of this with a concise and unambiguous API!
We should also be able to provide compile-time information such as ''lengths'' and maybe ''strides''.

=Some references=

* Some overviews what other programming languages do: [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/#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
* 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?]

=API playground=

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

===Source code===
<div class="toccolours mw-collapsible mw-collapsed">
''Click expand to see the implementation details.''

<div class="mw-collapsible-content">
<source lang="cpp">

#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
}


</source>
</div>
</div>

===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>
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
</pre>

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

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

{| 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
|-
|}

===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.-->
# 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>range(first,[incr,]last)</code> API:
#* Pros:
#** mimics MatLab first:[incr:]last API
#** enable indexing from the end
#* Cons:
#** 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...
-->
# [https://en.wikipedia.org/wiki/APL_(programming_language) 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
<!--# c++11 <code>mat({first,[incr,]last})</code> API:
#* Pros:
#** very compact
#** closely mimics MatLab first:incr:last API
#** no need to import a possibly polluting name as <code>Eigen::span</code>
#** enable indexing from the end
#* Cons:
#** 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
#** the argument order is not explicit-->
<!--# What does not work:
<source lang="cpp">
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 ...
</source>-->

=Internal implementation details=

==Performance==

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