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))
It can be used to let Eigen interface without any overhead with non-Eigen data structures, such as plain C arrays or structures from other libraries. By default, it assumes that the data is laid out contiguously in memory. You can however override this by explicitly specifying inner and outer strides.
Resizing. The current size of a matrix can be retrieved by rows(), cols() and size(). These methods return the number of rows, the number of columns and the number of coefficients, respectively. Resizing a dynamic-size matrix is done by the resize() method.
The default in Eigen is column-major.
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::Map
s and Eigen::Matrix
's interchangeably I found this the easiest and most robust method.
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With