It seems to me that it would make more sense to have the indices implemented as Arrays instead of Matrices, and the corresponding indexing for Array types would also be needed.

I was going to file a RFE for this today, but it's nice to see that you've already started on a patch!

Hi Michael,

this is a very good starting point! Of course I have many remarks though.

The first is that there is of course a large overlap between the current "Block" API (block, head, tail, segment, corner*, middle*, left*, right*) and what can be done with ranges. Personally I consider this as a plus because, depending on the context, one is more convenient than the other. Ranges are also more flexible because they include a step to extract what people sometimes call "slices", and they can be directly combined with col/row extraction and/or indices based extraction. On the other hand, "ranges" are currently slower than "blocks" because they lost direct access, vectorization and unrolling abilities.

Unrolling requires to know the size at compile-time. I guess it is fine to give up on this because this would lead to a cumbersome API:

Range(start,end)
Range<Size>(start)

without clear gain compared to what offers the current block based API.


DirectAccess is however important to preserve. This is possible when the nested expression has DirectAccess and when we only have single indices and/or ranges. In this case, MatrixViewSubscript should be based on MapBase<> with proper strides computed from the steps.

Knowing the "steps" at compile time is also very important, especially when step==1 for vectorization. I'd suggest to have the step as a template parameter of the class IndexRange defaulting to 1, and to add convenient functions to create ranges:

template<int Step> IndexRange<Step> Range(int s, int e)
{ return IndexRange<Step>(s,e); }

IndexRange<> Range(int s, int e) { return IndexRange<>(s,e); }

IndexRange<Dynamic> Range(int s, int e, int i) { return IndexRange<Dynamic>(s,e,i); }


In most cases, users would only have to deal with these Range() functions:

Range(1,5)
Range(1,5,2)
Range<2>(1,5)


The advantage of inheriting MapBase<> is that all the vectorization logic is already implemented.


Regarding the implementation of MatrixViewSubscript, I'm pretty sure that all the numerous variants could be factorized using some simple meta programming at a lower level.


MatrixViewIndex sounds quite dangerous since its result depend on the storage order.

I second this!

I could use it for following case:
I store a mesh of n Vector3f in a VectorXf of size 3*n. Parts of my code need to select a subset of these vertices (i.e. those in a certain bounding box), and then multiply a VectorXf constructed from them ({v1.x,v1.y,v1.z, v2.x,v2.y,v2.z, ...}) with a matrix.

Having the ability to use an index array to form a view on the original VectorXf would simplify my task quite a bit.

I would suggest to also support using a std::vector<int> in addition to a VectorXi to store the indices. This seems natural to me and would allow using iterators to construct it (e.g. I use a boost::filter_iterator<... bounding_box_filter ...> to create the vector).

What's the status on this feature? This is a central feature in MATLAB and Armadillo:

non-contiguous views (added in version 3.0):

X.elem( vector_of_indices )
X( vector_of_indices )           (added in version 3.810)

X.cols( vector_of_column_indices )
X.rows( vector_of_row_indices )

X( vector_of_row_indices, vector_of_column_indices )
X.submat( vector_of_row_indices, vector_of_column_indices )

I would appreciate if this went into Eigen as well. We use this kind of indexed views a lot.

(In reply to comment #4)
> What's the status on this feature?

It appears that nobody is actively working on this at the moment, maybe mostly because a clean solution is blocked by bug 99. We could implement an ad-hoc solution with no guarantee that the API may change later.

I increased the priority of this bug.

The biggest difference from Armadillo and MATLAB is in my experience the lack of indexed views in Eigen. That's the only missing killer-feature. I vote for an ad-hoc solution that might change in the future. Software is iterative, so I don't care if the API changes later on.

In our use we would also appreciate this feature. I think an ad-hoc implementation would be fine, one could use a new method name to avoid stepping on other parts of Eigen if it sounds too dangerous to overload the () operator for this.

As long as it is clearly stated to be experimental, we are OK with it changing later.

How is this bug going?

I'm working with FEM, so I solve part of a matrix and then add the results to a vector. In Armadillo I could do this:

vector(vector_of_indices) += matrix(vector_of_indices, vector_of_indices).solve(vector(vector_of_indices));

In Eigen I currently have to split this into like 20 rows of code where I actually have to create a new (dynamic) matrix and copy the entries according to vector_of_indices in it. The same for the right-hand-side vector. Then I have to loop over the results and add it to the resulting vector. This is not as performant as it could be with this feature added.

I know it's no use to nag but I just wanted to give you a real world example of how this could improve performance and simplicity of some projects.

(In reply to comment #8)
> How is this bug going?

Sorry to disappoint, but I'm afraid nobody is working on it yet.

> vector(vector_of_indices) += matrix(vector_of_indices,
> vector_of_indices).solve(vector(vector_of_indices));

In your case, you will merely safe one matrix-copy, because the decomposition must be on a copy anyways (you will save the copy from your temporary). 
As copying is O(N^2) but the decomposition is O(N^3), you will hardly notice a difference.
Also, solving the system will be much more efficient on an (aligned) temporary vector -- you might again save a copy operation here (but again copying is O(N), solving is O(N^2)).

The most important improvement is that you save your 20 rows of code (and I do think that is important).

> I know it's no use to nag [...]

It's actually like the second best way, to force a feature to be implemented.
The best way would be to implement it and provide a patch :)

We need to decide in what form vector_of_indices can be provided. The API-wise simplest solution would be to accept anything that has a size() and operator[] implemented. But we should ensure that operator[] returns any kind of integer types.

I agree that we should accept anything with a size() and operator[] returning integers. We should probably also specialize for operator[] returning bools to select the respective columns/rows.

To this end, without C++11 I guess that the only robust solution would be to check for either value_type or Scalar. Maybe we could only check for value_type and make value_type an alias of Scalar in Eigen (this was suggested to leverage more STL compatibility).

A simple test without requiring value_type would be:
template<typename T>
void assert_integer_type(const T&) { /* static assertion here */};

Usage:
template<class Indexes>
... foo(const Indexes& ind) {
  if(ind.size()>0) assert_integer(ind[0]);
  ...

I'm making this depend on Bug 360, anyways.

Another thing: An optimization we'd want to have is to check if Indexes itself has a fixed size, e.g., if we pass an std::array or an Array<int,Size,1> as index vector we should return an expression which itself has a fixed size. (This is definitely not required for the first versions ...)

The problem is that we probably want to specialize wrt integers or bools at compile time.

Yes, you are right. Checking for Scalar and value_type seems a reasonable solution for that. If we implement a traits class for that, this could also provide a SizeAtCompileTime value.

Any progress on this feature somewhere in the repository?

Now that Eigen 3.3 is almost there, this feature will be my top 2 priority for Eigen 3.4 (the top 1 being bug fixes of course). I've recently made a 20 lines example implementing such a feature:

http://eigen.tuxfamily.org/dox-devel/TopicCustomizing_NullaryExpr.html#title1

Thanks for this. My experience with CwiseNullaryOp for this sub-view capability was sub-optimal because it was not possible to change the entries of sub-matrix. In the end, I implemented both non-const and const views of matrices/vectors via derivation from MatrixBase class. I had to implement `operator=(const MatrixBase<Derived>& other)`. Then I created functions `rows(mat, irows)`, `cols(mat, icols)`, and `submatrix(mat, irows, icols)` that returned the row-, col-, or submatrix-view of a matrix `mat`. In this way, it was possible to do, for example:

~~~
rows(mat, irows).fill(0.0); // or
submatrix(mat, irows, icols) = ma1 + scalar*mat2; // assuming consistent dimensions
~~~

*** Bug 1334 has been marked as a duplicate of this bug. ***

I've started a proof of concept for the API, there: http://eigen.tuxfamily.org/index.php?title=Working_notes_-_Indexing%2B%2B

As briefly noted yesterday in IRC, I'd prefer {a,b,c} to be equivalent to Matlab [a,b,c]. Mapping {a,b,c} to anything like a:b:c could lead to confusion, especially if (implicit) conversions happens (or does not happen) somewhere. E.g., passing std::array<int>{a,b,c} or std::vector<int>{a,b,c} should be equivalent to passing {a,b,c}, IMO.

If we need something like Eigen::Range(a, b) for Matlab's a:b, we would automatically keep C++03 compatibility.

Regarding API, especially for the 3-parameter Range function/object, it will be difficult to give an intuitive parameter meaning/order 
  Range(start, step, end)
  Range(size, start, end)
  Range(size, start, step)
  Range(start, end, step)
(and many more, e.g. step and size could be template parameters instead), could make sense and since all parameters are integers, we can't really do compile time checks to avoid ambiguity (without getting very verbose).

Regarding end (or last): Matlab/Octave supports almost arbitrary expressions involving end, e.g.:
>> A=1:10; A(1:20/end:ceil(sqrt(end)))
ans = 1  3

I.e., essentially end is directly substituted by the corresponding dimension of A, which I guess we won't be able to fully port to C++

Are there plans to support these views as lvalues? It would be great to have the following behavior:

~~~
VectorXd v = {1, 2, 3}; // support to initializer_list would be great
v.rows({0, 2}) = {100, 300};     // v = [100, 2, 300]
VectorXd w = v.rows({2, 1, 0});  // w = [300, 2, 100]
~~~

yes, that's the plan, but without the ".rows":

~~~
VectorXd v = {1, 2, 3};     // support to initializer_list would be great
v({0, 2}) = {100, 300};     // v = [100, 2, 300]
VectorXd w = v({2, 1, 0});  // w = [300, 2, 100]
~~~

Great! What if one wants to get a view of some rows or columns of a matrix. How would the syntax look like?

~~~
MatrixXd A = MatrixXd::Random(10, 10);
MatrixXd B = A({0, 1});      // not very clear/intuitive what this would mean
MatrixXd C = A.rows({0, 1}); // C is 2x10
MatrixXd D = A.cols({0, 1}); // D is 10x2
~~~

(In reply to Allan Leal from comment #23)
> Great! What if one wants to get a view of some rows or columns of a matrix.
> How would the syntax look like?

Something like this:
  MatrixXd C = A(Eigen::all, {0,1});  // Matlab: C = A(:, [1 2]);
  MatrixXd D = A({0,1}, Eigen::all);  // Matlab: D = A([1 2], :);


The final syntax and keyword names are still debatable, of course.

This is all very great guys. If possible, it would be nice to have, in addition to your proposal: 

~~~
  MatrixXd C = A(Eigen::all, {0,1});  // Matlab: C = A(:, [1 2]);
  MatrixXd D = A({0,1}, Eigen::all);  // Matlab: D = A([1 2], :);
~~~

a syntax that makes a parallel to other similar methods:

~~~
  MatrixXd C = A.indexedRows({0, 1}); // in parallel with topRows, bottomRows.
  MatrixXd D = A.indexedCols({0, 1}); // in parallel with leftCols, rightCols.
~~~

But I understand if the common understanding is to not extend to much the interface of MatrixBase class. After all, one can always write a plugin to such class and extend it with other methods.

Because this is a central feature for Eigen's user, most of the discussions occurred on the public ML. For the record, the implementation is very well advanced and fully usable there:

https://bitbucket.org/ggael/eigen-flexidexing

You can see the usages in the respective unit test:

https://bitbucket.org/ggael/eigen-flexidexing/src/4e045c1fd2920470817aeb3de50a2d1e7b101a84/test/indexed_view.cpp?at=default&fileviewer=file-view-default#indexed_view.cpp-79

API and code reviews are very welcomed. You can see the diff here:

https://bitbucket.org/ggael/eigen-flexidexing/branches/compare/ggael/eigen-flexidexing:default%0Deigen/eigen:default#diff

I forgot to update this entry, but the fork mentioned in the previous comment has been merged in de default branch a long time ago. It only remains to update the manual pages. Currently,  you can look at operator()() doc:

http://eigen.tuxfamily.org/dox-devel/classEigen_1_1DenseBase.html#a0b44220621cd59a75cd0f48cc499518f

and unit tests for examples:

https://bitbucket.org/eigen/eigen/src/8d1ccfd9c5a0bf1a662bd0fe35858f4e527d2bac/test/indexed_view.cpp?at=default&fileviewer=file-view-default

This feature seems mature enough, but it is not yet available in 3.3.5. Are there any plans in releasing a new tag soon with indexed views supported? This is needed to ensure proper version of Eigen is used and also when using some package manager to install eigen (e.g., conda, conan, etc.).

This likely won't get backported to the 3.3.x branch. You need to use the devel-branch or wait for 3.4 for this.

-- GitLab Migration Automatic Message --

This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.

You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/329.