Eigen array shift gives strange result-CodePudding

I'm almost sure it's a trivial problem, but I can't understand what wrong is going on. I have a simple array in Eigen that I need to right-shift, so that I can always insert new sample at position 0:

#include <Eigen/Dense>

using namespace Eigen;
typedef VectorXd vec;

// The actual buffer
int N = 10;
vec _x(N);
_x = vec::Constant(0.0);

// This is just for dimonstrative purpose
for (int i = 0; i < N; i  )
    _x(seq(1, N - 1)) = _x(seq(0, N - 2));
    _x(0) = (double) i;
}

What I got is that, at each cycle, it is copying twice each element to be shifted:

--- After cycle 1: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
--- After cycle 2: [1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
--- After cycle 3: [2, 1, 1, 0, 0, 0, 0, 0, 0, 0]
--- After cycle 4: [3, 2, 2, 1, 1, 0, 0, 0, 0, 0]
--- ...
--- After cycle 10: [9, 8, 8, 7, 7, 6, 6, 5, 5, 4]

I'm pretty sure I'm messing something around with Eigen. What's wrong with this approach?

CodePudding user response：

It was tricky and trivial at the same time. The problem is into the aliasing that Eigen vectors undergo when performing block assignment.

The solution is to compute .eval() at the right-hand operator in order to let Eigen evaluate the new array first, and only then assign it to the left-hand block:

_x(seq(1, N - 1)) = _x(seq(0, N - 2)).eval();

Hope this can be useful for someone

CodePudding user response：

As you found out yourself, you have an aliasing problem -- which can be solved by adding an eval() at the end of the expression. However, this (usually) requires an extra memory allocation. If you want to move the elements in-place, you can use std::copy_backward:

std::copy_backward(_x.begin(), _x.end()-1, _x.end());

But this is also not optimal, since it likely does not vectorize well.

The actual solution to your problem should be to re-think how you can collect your data without always having to copy everything.

For example, have a bigger array which you start filling backwards from position end()-10, and always just use a .segment() of that array.

If you only rarely need to actually use the array, use some kind of ring-buffer which you only copy to an Eigen structure when needed.

CodePudding user response：

First, as other say, aliasing is the problem, don't get caught in the "easy syntax" of the library. Expression templates can lead to bad code. Also mutation and aliasing are incompatible.

My solution would be to forget about this feature and go back to the basics of algorithms. If you want to shift elements, well std::shift_ the elements then:

#include <Eigen/Dense>
#include <iostream>

int main() {
    using namespace Eigen;
    typedef VectorXd vec;

    // The actual buffer
    int N = 10;
    vec _x(N);

    // This is just for dimonstrative purpose
    for (int i = 0; i < N; i  ) {
    //    _x(seq(1, N - 1)) = _x(seq(0, N - 2));
        std::shift_right(_x.begin(), _x.end(), 1);
        _x(0) = static_cast<double>(i);

        std::cout << '[' << _x << ']' << std::endl;
    }

}

https://godbolt.org/z/rnhsPxdsE

more idiomatic:

    for(int i = 0; i != N; i  ) {
        *(std::shift_right(_x.begin(), _x.end(), 1) - 1) = static_cast<double>(i);
    }