Hello,

I'm working on converting some Matlab code to Eigen. In several places in the code I have a vector of indices, and I need to get the corresponding elements from another vector. For example, if v = [10; 20; 30; 40; 50] and i = [0; 1; 4] then I need to be able to get the value of v(i) (in this case [10; 20; 50]). There doesn't seem to be any built-in functions in Eigen for doing this, so I wrote the following function:

Code: Select all

// v_idx(res, vec, idx):  res = vec(idx)
template <typename D1, typename D2, typename D3>
void v_idx(MatrixBase<D3> EIGEN_REF_TO_TEMPORARY _res,
           MatrixBase<D1> const& vec, MatrixBase<D2> const& idx) {
  const ArrayWrapper<D1> a = vec.array();
  const ArrayWrapper<D2> x = idx.array();
  MatrixBase<D3> &res = const_cast< MatrixBase<D3>& >(_res);
  res.derived().resize(x.size());
  ArrayWrapper<D3> r = res.array();
  for (size_t i=0; i < x.size(); i++) {
    r(i) = a(x(i));
  }
}

So now I can call v_idx(x, v, i) and x will be set to [10; 20; 50]. This works fine, but it's a little cumbersome to use, because I have to pass in a result vector as one of the arguments every time. When I convert more complicated expressions from the Matlab code, I end up having to add several temporary variables, which I don't want to do if I don't have to.

Is there any way to rewrite this function so that it returns a vector without introducing unnecessary temporaries? I guess it would involve using expression templates somehow, but I'm having trouble understanding that part of the documentation...

Also, how would I write a function to replace the elements of v(i) with another vector? For example, if I write something like v.idx(i) = w (where w = [100; 200; 300]) then that should make v = [100; 200; 30; 40; 300].

Thanks!

ccrawford wrote:Is there any way to rewrite this function so that it returns a vector without introducing unnecessary temporaries?

I can point you to the ReturnByValue class in Eigen 3. There are quite a few examples in the code. Here is a rough sketch of what you have to do...

The idea is, that you write a struct inheriting ReturnByValue. That struct implements a function evalTo(..) which is used during the assignment of that struct to another Matrix or Array. The actual implementation of your function will reside within evalTo(..). The other important thing is that you need to provide a typedef that tells the compiler in which type the object should be evaluated when it is used as a temporary - e.g. when passing it to something like std::cout. So where does the struct come from? You initial function will create and return it as a proxy object.

In the end it will look similar to this

Code: Select all

// this guy can be assigned to matrices!!
template <typename D1, typename D2>
struct vector_idx_returnbyvalue_type : public ReturnByValue< ... >
{
 // ...
};

template <typename D1, typename D2>
vector_idx_returnbyvalue_type
v_idx(MatrixBase<D1> const& vec, MatrixBase<D2> const& idx) {
  return vector_idx_returnbyvalue_type(vec,idx);
}


As said before, for the rest I am kindly asking you to take a look at the source code. Just come back if you have problems.

Regards,
Hauke

When I asked for an explanation of ReturnByValue, Benoit kindly wrote me the following example which I found very helpful: http://listengine.tuxfamily.org/lists.t ... 00216.html .

Incidentally, I remember that we had a couple of people already asking about the functionality that ccrawford mentions, so it would be very useful to have an implementation.

At the same time, in this particular example, I think if this is best implemented as a new expression class, enabling it to be completely lazy. (The downside of that, though, would be the indirect memory access and potentially poor memory locality.)

There already have been some discussions on the list about adding such an expression, a few months ago.

Thanks Hauke and jitseniesen for your help. After a couple of days wrangling with compiler errors and segfaults, I came up with the following implementation, which seems to be working:

Code: Select all

template <typename D1, typename D2> struct v_idx_retval;

namespace Eigen {
  template <typename D1, typename D2> struct ei_traits<v_idx_retval<D1, D2> > {
    typedef Matrix<typename D1::Scalar, D2::RowsAtCompileTime, 1, D1::Options, D2::MaxRowsAtCompileTime, 1> ReturnType;
  };
}

template <typename D1, typename D2>
  struct v_idx_retval : public ReturnByValue<v_idx_retval<D1, D2> > {
 v_idx_retval(D1 const& vec, D2 const& idx) : m_vec(vec), m_idx(idx) {
  }
  inline size_t rows() const { return m_idx.size(); }
  inline size_t cols() const { return 1; }
  template<typename Dest> inline void evalTo(Dest& dst) const {
    dst.resize(m_idx.size());
    for (size_t i=0; i<m_idx.size(); i++) {
      dst.coeffRef(i) = m_vec(m_idx.coeff(i));
    }
  }
 protected:
  const D1& m_vec;
  const D2& m_idx;
};

template <typename D1, typename D2>
  v_idx_retval<D1, D2> v_idx2(MatrixBase<D1> const& vec,
               MatrixBase<D2> const& idx) {
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(D1)
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(D2)
    return v_idx_retval<D1, D2>(vec.derived(), idx.derived());
}