MATLAB find() / Numpy nonzero idioms for Eigen
Asked Answered
B

3

5

Chances are this is a very stupid question but I spent a pretty absurd amount of time looking for it on the documentation, to no avail.

in MATLAB, the find() function gives me an array with the indices of nonzero elements. Numpy's np.nonzero function does something similar.

How do I do this in the C++ Eigen library? I have a Boolean array of

typedef <bool, 10, 1> foobar = MatrixA < MatrixB;

so far. Thanks!

Bedford answered 27/4, 2013 at 7:1 Comment(1)
IMHO, Unlike MATLAB or Numpy, Eigen doesn't really need a find() function' since you can easily write it yourself in C++ using loops or whatnot. If you did that in MATLAB or Numpy, there'd be a big performance penalty.Stymie
S
2

It is reasonable to expect Eigen to have a find() function. Unfortunately, Eigen doesn't have one, or even a less than operator for matrices. Fortunately, the problem isn't too difficult. Here is one solution to the problem. I am using vector to store the Column Major indices of elements > 0. You could use VectorXf if you prefer that. Use this on B - A (B-A > 0 is the same as evaluating B>A). I'm using the stl for_each() function.

#include<algorithm>
#include<vector>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

class isGreater{
    public:
    vector<int>* GT;
    isGreater(vector<int> *g){GT = g;}
    void operator()(float i){static int it = 0; if(i>0)GT->push_back(it); it++;}
};
int main(int argc,char **argv){
    MatrixXf P = MatrixXf::Random(4,5);
    vector<int> GT;
    for_each(P.data(),P.data()+P.rows()*P.cols(),isGreater(&GT));
    cout<<P<<endl;
    for(int i=0;i<GT.size();++i)cout<<GT[i]<<" ";
    cout<<GT.size()<<endl;
    return 0;
}
Sapindaceous answered 1/7, 2013 at 14:54 Comment(0)
S
8

Not sure if this is part of your question, but to construct the appropriate element-wise inequality result you must first cast your matrices to arrays:

MatrixXd A,B;
...
Matrix<bool,Dynamic,Dynamic> C = A.array()<B.array();

Now C is the same size as A and B and C(i,j) = A(i,j) < B(i,j).

To find all of the indices (assuming column-major order) of the true entries, you can use this compact c++11 routine---as described by libigl's conversion table:

VectorXi I = VectorXi::LinSpaced(C.size(),0,C.size()-1);
I.conservativeResize(std::stable_partition(
  I.data(), I.data()+I.size(), [&C](int i){return C(i);})-I.data());

Now I is C.nonZeros() long and contains indices of the true entries in C. These two lines essentially implement find.

Sherard answered 2/2, 2015 at 15:33 Comment(1)
this implementation is very smart but if one wants to emulate [i j]=find(A) where i,j are index of row and column where a matrix A is true, how can this be done?Neolith
S
2

It is reasonable to expect Eigen to have a find() function. Unfortunately, Eigen doesn't have one, or even a less than operator for matrices. Fortunately, the problem isn't too difficult. Here is one solution to the problem. I am using vector to store the Column Major indices of elements > 0. You could use VectorXf if you prefer that. Use this on B - A (B-A > 0 is the same as evaluating B>A). I'm using the stl for_each() function.

#include<algorithm>
#include<vector>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

class isGreater{
    public:
    vector<int>* GT;
    isGreater(vector<int> *g){GT = g;}
    void operator()(float i){static int it = 0; if(i>0)GT->push_back(it); it++;}
};
int main(int argc,char **argv){
    MatrixXf P = MatrixXf::Random(4,5);
    vector<int> GT;
    for_each(P.data(),P.data()+P.rows()*P.cols(),isGreater(&GT));
    cout<<P<<endl;
    for(int i=0;i<GT.size();++i)cout<<GT[i]<<" ";
    cout<<GT.size()<<endl;
    return 0;
}
Sapindaceous answered 1/7, 2013 at 14:54 Comment(0)
P
1

This might work for you and others who check this out. In order to set elements of a matrix m based on the condition on another matrix A, you can use this notation:

m = (A.array() != 0).select(1, m);

This command replaces those elements in matrix m that have non-zero corresponding elements in A, with one.

Playa answered 29/1, 2016 at 18:3 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.