Search code examples
c++arraysmatlabsliceeigen3

Is there a better way to implement matlab's Logical Indexing using Eigen/C++?


I've been looking all over for a more "eigen" way of implementing the functionality of Matlab's Logical indexing. Here's the best I could come up with. (Focusing on an int array here for simplicity)

//an attempt at matlab-style Logical Indexing
//equivalent to the matlab:
//  original = [1,2,3,4]
//  subset = original(original < 3)

using namespace Eigen;
using std::cout;
using std::endl;

IOFormat OctaveFmt(StreamPrecision, 0, ", ", " ", "", "", "[", "]");

ArrayXi original(4);
original << 1,2,3,4;
cout<<"Original with bad values:"<<endl
    <<original.format(OctaveFmt)<<endl;

Array<bool, Dynamic,1> selections = original < 3;
cout<<"One if it's a good value:"<<endl
    <<selections.format(OctaveFmt)<<endl;

std::vector<int> picked;
for(int i = 0; i < selections.size(); i++ )
{
    if(selections(i))
    {
        picked.push_back(original(i));
    }
}

//put the vector values back into an eigen array
ArrayXi theGoodStuff = Map<ArrayXi, Unaligned>
                       (picked.data(), picked.size());

cout<<"Just the good stuff:"<<endl
    <<theGoodStuff.format(OctaveFmt)<<endl;

Here's the output I get:

Original with bad values:
[1  2  3  4]
One if it's a good value:
[1  1  0  0]
Just the good stuff:
[1  2]

Does anyone know how to do this in a more 'eigen' way, or just a faster way than looping through the arrays?


Solution

  • This is an old one but....Eigen 3.4 makes it way cleaner to do this. Sadly I can't speak to the performance. I'm mainly interested readability. I made this:

    class logical
    {
    private:
        const Index new_size;
        Array<Index, Dynamic, 1> old_inds;
    
    public:
        logical(const Array<bool, Dynamic, 1> &keep) : new_size(keep.count()), old_inds(new_size)
        {
            for (Index i = 0, j = 0; i < keep.size(); i++)
                if (keep(i))
                    old_inds(j++) = i;
        }
        Index size() const { return new_size; }
        Index operator[](Index new_ind) const { return old_inds(new_ind); }
    };
    

    Which is used like:

    logical inds(ind_to_keep);
    Y_slice = Y(inds); // vector
    H_slice = H(inds, all); // remove some rows
    R_slice = R(inds, inds); // remove the same rows as columns
    

    Eigen help is at the bottom of this page: https://eigen.tuxfamily.org/dox-devel/group__TutorialSlicingIndexing.html