Using Eigen::Map<Eigen::MatrixXd> as function argument of type Eigen::MatrixXd

Question

In short, the question is how to pass a

Eigen::Map<Eigen::MatrixXd>

object to a function which expects a

Eigen::MatrixXd

object.

---

Longer story:

I have this C++ function declaration

void npMatrix(const Eigen::MatrixXd &data, Eigen::MatrixXd &result);

together with this implementation

void npMatrix(const Eigen::MatrixXd &data, Eigen::MatrixXd &result)
{
//Just do s.th. with arguments
std::cout << data << std::endl;

result(1,1) = -5;
std::cout << result << std::endl;
}

I want to call this function from python using numpy.array as arguments. To this end, I use a wrapper function written in c++

void pyMatrix(const double* p_data, const int dimData[],
                              double* p_result, const int dimResult[]);

which takes a pointer to data, the size of the data array, a pointer to result, and the size of the result array. The data pointer points to a const patch of memory, since data is not to be altered while the patch of memory reserved for result is writeable. The implementation of the function

void pyMatrix(const double *p_data, const int dimData[], double *p_result, const int dimResult[])
{
Eigen::Map<const Eigen::MatrixXd> dataMap(p_data, dimData[0], dimData[1]);
Eigen::Map<Eigen::MatrixXd> resultMap(p_result, dimResult[0], dimResult[1]);

resultMap(0,0) = 100;

npMatrix(dataMap, resultMap);
}

defines a Eigen::Map for data and result, respectively. A Eigen::Map allows to access raw memory as a kind of Eigen::Matrix. The dataMap is of type

<const Eigen::MatrixXd>

since the associated memory is read only; resultMap in contrast is of type

<Eigen::MatrixXd>

since it must we writeable. The line

resultMap(0,0) = 100;

shows, that resultMap is in deed writeable. While passing dataMap to the npMatrix() where a const Eigen::MatrixXd is expected works, I could not find a way to pass resultMap in the same way. I am sure, the trouble comes from the fact, that the first argument of npMatrix is const, and the second is not. A possible solution I found is to define

Eigen::MatrixXd resultMatrix = resultMap;

and pass this resutlMatrix to npMatrix(). However, I guess, this creates a copy and hence kills the nice memory mapping of Eigen::Map. So my question is.

Is there a way to pass a Eigen:Map to a function which expects a non-const Eigen::MatrixXd instead?

As a side note: I could change npMatrix to expect a Eigen::Map, but since in the real project, functions are already there and tested, I would rather not temper with them.

To complete the question, here is the python file to call pyMatrix()

import ctypes as ct
import numpy as np
import matplotlib.pyplot as plt

# Load libfit and define input types
ct.cdll.LoadLibrary("/home/wmader/Methods/fdmb-refactor/build/pyinterface/libpyfit.so")
libfit = ct.CDLL("libpyfit.so")

libfit.pyMatrix.argtypes = [np.ctypeslib.ndpointer(dtype=np.float64, ndim=2),
                                                     np.ctypeslib.ndpointer(dtype=np.int32, ndim=1),
                                                     np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, flags='WRITEABLE'),
                                                     np.ctypeslib.ndpointer(dtype=np.int32, ndim=1)
                                                     ]

data = np.array(np.random.randn(10, 2), dtype=np.float64, order='F')
result = np.zeros_like(data, dtype=np.float64, order='F')

libfit.pyMatrix(data, np.array(data.shape, dtype=np.int32),
                              result, np.array(result.shape, dtype=np.int32))

Answer 1

Pass it as plain pointer to data, and Eigen::Map it there. Alternatively, use template <typename Derived> and the like, found in http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html My personal Choice is the first though, as it is better to have code that doesn't expose all the stubbornness of every API you have used. Also, you won't lose compatibility neither with eigen ,nor with any other kind of library that you (or anyone else) may use later.

There is also another trick i found out, which can be used in numerous occasions:

Eigen::MatrixXd a; //lets assume a data pointer like double* DATA that we want to map //Now we do new (&a) Eigen::Map<Eigen::Matrix<Double,Eigen::Dynamic,Eigen::Dynamic>> (DATA,DATA rows,DATA cols);

This will do what you ask, without wasting memory. I think it is a cool trick, and a will behave as a matrixXd, but i haven't tested every occasion. It has no memory copy. However, you might need to resize a to the right size before assigning. Even so, the compiler will not immediately allocate all memory at the time you request the resize operation, so there won't be big useless memory allocations either!

Be careful! Resizing operations might reallocate the memory used by an eigen matrix! So, if you ::Map a memory but then you perform an action that resizes the matrix, it might be mapped to a different place in memory.