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))

Inhabiter

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.

Henleigh

Thanks for your answer. I accepted it since I think that, in general, this is a good solution. As I stated in the question, I do not want to change the interface of the function to call. Therefore, my solution is to create the Map I need. This results in waste of memory and cpu time, but compared to the rest of the program it does not matter. – Inhabiter

Thanks! I have thought it some more, and i think i have found another answer that will do exactly what you ask, and not change the function declaration. I will update my answer! – Henleigh

It's been 18 months, but SO has so much useful info. Here new (&a) is the placement version, and thus casting from Map to Matrix object implicitly I believe (please correct me if I am wrong). I guess the resizing you are referring to is that a is not the owner of the data, so cannot assign another matrix to it, using say, assignment operator=. I do this by keeping track of ownership of the data - initially a may not own the data, but before assignment, I resize it, when it is NOT the owner or the current size is not correct. – Georgetown

Answer 2

For anyone still struggling with the problem of passing an Eigen::Map to a function with signature Eigen::Matrix or vice versa, and found the Eigen::Matrix to Eigen::Map implicit casting trick which @Aperture Laboratories suggested didn't work (in my case this gave runtime errors associated with trying to free already released memory, [Mismatched delete / Invalid delete errors when ran with valgrind]),

I suggest using the Eigen::Ref class for function signatures as suggested in the answer given by @ggael here: Passing Eigen::Map<ArrayXd> to a function expecting ArrayXd&

and written in the documentation: http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html#TopicUsingRefClass under the title:

How to write generic, but non-templated function?

For example, for the function specified in the question, changing the signature to

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

means passing either Eigen::Map<Eigen::MatrixXd> orEigen::MatrixXd objects to the function would work seamlessly (see @ggael's answer to Correct usage of the Eigen::Ref<> class for different ways to use Eigen::Ref in function signature).

I appreciate OP said he didn't want to change the function signatures, but in terms of using Eigen::Maps and Eigen::Matrix's interchangeably I found this the easiest and most robust method.

Lumbard