0

The following code concatenates a vector of ones to a matrix:

using Eigen::VectorXd;
using Eigen::MatrixXd;

MatrixXd cbind1(const Eigen::Ref<const MatrixXd> X) {
  const unsigned int n = X.rows();
  MatrixXd Y(n, 1 + X.cols());
  Y << VectorXd::Ones(n), X;
  return Y;
}

The function copies the contents of X to Y. How do I define the function so that it avoids doing the copy and returns a reference to a matrix containing VectorXd::Ones(n), X?

Thanks.

user3294195
  • 1,748
  • 1
  • 19
  • 36
  • That seems to me impossible. Since you do not want to concatenate a vector to **Original Matrix X**, the only difference to make is a deep copy or shallow copy , or something like that. – xtluo Nov 19 '16 at 02:30
  • I've edited my answer to your [previous question](http://stackoverflow.com/questions/40647219) on that matter showing that they are actually duplicates. – ggael Nov 20 '16 at 20:14

1 Answers1

3

If you had followed and read ggael's answer to your previous question (and it appears you did, as you accepted his answer), you would have read this page of the docs. By modifying the example slightly, you could have written as part of an MCVE:

#include <iostream>
#include <Eigen/Core>

using namespace Eigen;

template<class ArgType>
struct ones_col_helper {
    typedef Matrix<typename ArgType::Scalar,
        ArgType::SizeAtCompileTime,
        ArgType::SizeAtCompileTime,
        ColMajor,
        ArgType::MaxSizeAtCompileTime,
        ArgType::MaxSizeAtCompileTime> MatrixType;
};

template<class ArgType>
class ones_col_functor
{
    const typename ArgType::Nested m_mat;

public:
    ones_col_functor(const ArgType& arg) : m_mat(arg) {};

    const typename ArgType::Scalar operator() (Index row, Index col) const {
        if (col == 0) return typename ArgType::Scalar(1);
        return m_mat(row, col - 1);
    }
};

template <class ArgType>
CwiseNullaryOp<ones_col_functor<ArgType>, typename ones_col_helper<ArgType>::MatrixType>
cbind1(const Eigen::MatrixBase<ArgType>& arg)
{
    typedef typename ones_col_helper<ArgType>::MatrixType MatrixType;
    return MatrixType::NullaryExpr(arg.rows(), arg.cols()+1, ones_col_functor<ArgType>(arg.derived()));
}

int main()
{
    MatrixXd mat(4, 4);
    mat << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16;

    auto example = cbind1(mat);

    std::cout << example << std::endl;
    return 0;
}
Community
  • 1
  • 1
Avi Ginsburg
  • 10,323
  • 3
  • 29
  • 56
  • My compiler is unable to identify the type `Index` in `one_col_functor`. Is it part of the `Eigen` namespace? – user3294195 Nov 20 '16 at 18:22
  • Yes, Eigen::Index – Avi Ginsburg Nov 20 '16 at 18:25
  • It says no type named `Index` in namespace `Eigen`. Isn't it usually something like `Eigen::VectorXd::Index`? – user3294195 Nov 20 '16 at 18:26
  • Are you using Eigen 3.3? – Avi Ginsburg Nov 20 '16 at 18:27
  • I'm using `RcppEigen`, which uses Eigen 3.2.9. – user3294195 Nov 20 '16 at 18:28
  • I just added `namespace Eigen { typedef std::ptrdiff_t Index; }`. That seems to do the trick. – user3294195 Nov 20 '16 at 18:38
  • Actually the answer to the previous question also solve that one, as seen in my [edited answer](http://stackoverflow.com/a/40664512/1641621). – ggael Nov 20 '16 at 20:16
  • BTW, the line `const ArgType &m_mat` is not correct in case you nest a temporary expression. It must be `const typename ArgType::Nested m_mat`. – ggael Nov 21 '16 at 14:12
  • @ggael So the `NestByRefBit` determines if a copy should be made, which in the case of an expression (temporary or not) resolves to false so a copy is made, while in the case of a matrix is true so it's "stored" as a reference? – Avi Ginsburg Nov 21 '16 at 15:07
  • Yes exactly. This is handled internally by `internal:ref_selector`. Of course, if you nest a temporary matrix, then you are screwed. For instance `auto X = cbind1((A*B).eval());` is wrong. – ggael Nov 23 '16 at 08:49