## Cython Journey Part 2: Including Eigen

As I mentioned in my last Cython Journey post, the reason I’m using Cython is to take my speedy C++ code that does a bunch of matrix operations using the Eigen library, and to make them accessible from Python. Now, I have no desire to make an interface for the Eigen library, all I need is to be able to access the `MatrixXf` objects from the Eigen library, and convert them into `numpy.ndarray` objects for manipulation in Python. The logical first step towards this end would then seem to be to successfully compile C++ code that links to the Eigen library.

HOWEVER. After messing around with a bunch of linking problems and the like, I came across code by charanpald that provides an interface for Eigen sparse matrices and sparse matrix operations to Python (https://github.com/charanpald/SpPy). In this code, he links to the Eigen library through a buffer class that extends the `Eigen::SparseMatrix` class and adds some more functionality. There is very possibly a better way to do this, but I also needed to define some extra functions operating on the Eigen matrices, so defining an extension to the `MatrixXf` class worked out well. For this example, though, I’ll just use a barebones class extension to keep the size manageable. Here is my class:
cpp_matrixxfpy.h

```#ifndef MATRIXXFPY_H
#define MATRIXXFPY_H
#include <Eigen/Dense>
using namespace Eigen;

class MatrixXfPy : public MatrixXf {
public:
MatrixXfPy() : MatrixXf() { }
MatrixXfPy(int rows,int cols) : MatrixXf(rows,cols) { }
MatrixXfPy(const MatrixXf other) : MatrixXf(other) { }
};
#endif
```

Like, I said, barebones.
With this now we can go to our .pyx file and create a Cython handle into the `MatrixXf` object.
matrixxfpy.pyx

```cdef extern from "cpp_matrixxfpy.h":
cdef cppclass MatrixXfPy:
MatrixXfPy()
MatrixXfPy(int d1, int d2)
MatrixXfPy(MatrixXfPy other)
int rows()
int cols()
float coeff(int, int)
```

Same as for the `cpp_test` code that we wanted to interface with Python; first thing we do in the .pyx file is redefine the class and variables / functions that we want access to in Python, along with the constructors. In this simple example we’re going to only access three functions from `MatrixXf`, `rows()`, `cols()`, and `coeff(int, int)`. These functions return the number of rows, the number of columns, and the matrix coefficient value at a given index.

So now that we have a Cython handle on the Eigen code, let’s make this functionality available to Python. HOWEVER. Let’s keep in mind the goal, we’re not here to provide a Python wrapper for `MatrixXf` objects. What we really want is just to have Eigen and our C++ code over there, somewhere else, doing all of the calculations etc and then for us to be over in Python and just be able to get a `numpy.ndarray` for Python to play with and plot. Let us set about that now.

Let’s go back to the Test class we defined in Cython Journey Part 1. First thing is first, we need to edit the C++ code to create and return a `MatrixXf` object. So let’s just throw in a simple function:
cpp_test.h

```#ifndef TEST_H
#define TEST_H

#include "cpp_matrixxfpy.h"
using namespace Eigen;

class Test {
public:
int test1;
Test();
Test(int test1);
~Test();
int returnFour();
int returnFive();
Test operator+(const Test& other);
Test operator-(const Test& other);
MatrixXfPy getMatrixXf(int d1, int d2);
};
#endif
```

and define it in cpp_test.cpp as

```#include "cpp_test.h"

Test::Test() {
test1 = 0;
}

Test::Test(int test1) {
this->test1 = test1;
}

Test::~Test() { }

int Test::returnFour() { return 4; }

int Test::returnFive() { return returnFour() + 1; }

Test Test::operator+(const Test& other) {
return Test(test1 += other.test1);
}

Test Test::operator-(const Test& other) {
return Test(test1 -= other.test1);
}

MatrixXfPy Test::getMatrixXf(int d1, int d2) {
MatrixXfPy matrix = (MatrixXf)MatrixXfPy::Ones(d1,d2);
matrix(0,0) = -10.0101003; // some manipulation, to show it carries over
return matrix;
}
```

I put up the whole class again so it’s all in one place.

Alright, that’s done. Now we can look at our .pyx file for `Test`. We’re going to need to make the `MatrixXfPy` class defined just above available to our `Test` class. To do this is easy, all we need to throw in is:
test.pyx

```
include "cpp_matrixxfpy.h"
```

at the top! And while we’re at it, we’re going to be adding some more functionality to test that requires `numpy`, so also throw in:

```
import numpy
```

This additional functionality is going to access the `getMatrix(int,int)`, create an `numpy.ndarray` object, and fill it out with the values from the matrix. Let’s get about doing that. We need to expose the new C++ function to Cython, so add this line to the `cdef cppclass Test` declarations:

```        MatrixXfPy getMatrixXf(int d1, int d2)
```

And now we can add the code that will be accessible from Python. I’ll post the code and then describe it below (although it’s pretty straightforward):

```    def getNDArray(self, int d1, int d2):
cdef MatrixXfPy me = self.thisptr.returnMatrixXf(d1,d2) # get MatrixXfPy object

result = numpy.zeros((me.rows(),me.cols())) # create nd array
# Fill out the nd array with MatrixXf elements
for row in range(me.rows()):
for col in range(me.cols()):
result[row, col] = me.coeff(row, col)
return result
```

I just flat out stripped the form for this from the SpPy project and simplified it, so it’s not going to be optimized for speed in any way, but it gets the job done and it’s very easy to read. First we just set up a `numpy.zeros` matrix to index into, and then we go through each element in the array and fill it in with the appropriate value by calling to the `coeff(int row, int col)` function we defined above. Easy!

Now we can use the same setup file from the last Cython post to build our shared library:

```from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

setup(
name = 'Test',
ext_modules=[
Extension("test",
sources=["test.pyx", "cpp_test.cpp"], # Note, you can link against a c++ library instead of including the source
language="c++"),
],
cmdclass = {'build_ext': build_ext},

)
```

Again reproduced here for hassle-freeness. Now we can compile and test it! At the IPython terminal:

```run setup.py build_ext -i
import test
t = test.pyTest(10)
m = t.getNDArray(4,4)
```

And that’s it! We can see that that last command returns to us a `(4x4)` matrix with the first element set to `-10.0101003`, which is exactly what we did over in our C++ files. Now we’re all set to go ahead and generate matrices in C++, and then ship ’em out over to Python for plotting. Excellent.

Tagged , ,

## N-dimensional Matrices in C++

I’m working on coding up an algorithm in C++, based on a Matlab implementation I wrote, and I’ve been using the Eigen library (http://eigen.tuxfamily.org/) for my linear algebra. It’s been relatively straightforward to use so far, but unfortunately it doesn’t (as far as I can tell) have any built in support for N-dimensional matrices. In fact, it appears that N-dimensional matrices aren’t really addressed by any of the standard C++ linear algebra libraries. So, to get my 4-dimensional matrix, I have to make an array of pointers pointing to an array of pointers pointing to my Eigen MatrixXf objects. I initialize my 4-D matrices using this function:

```MatrixXf*** initialize_4dMatrix(MatrixXf*** matrix, int sizes[], float value=0.0) {
matrix = new MatrixXf** [sizes];
for (int d1 = 0; d1 < sizes; d1++) {
matrix[d1] = new MatrixXf* [sizes];
for (int d2 = 0; d2 < sizes; d2++)
matrix[d1][d2] = new (MatrixXf)(MatrixXf::Ones(sizes,sizes).array() * value);
}
return matrix;
}
```

and call it with:

```int sizes = {1,2,3,4};
MatrixXf ***m = initialize_4dMatrix(m, {1,2,3,4});
```

To then access the matrices, it’s `(*m[d1][d2])(d3,d4);` which isn’t bad for accessing one matrix at a time, but what if I want the jth element of the ith row of every matrix (i.e. element (i,j))? The Python access method `m[:,:,i,:]` (or just `m[:,:,i]`) would be ideal. Unfortunately, I don’t think there’s any way around writing another function to do it:

```    static VectorXf* get_elements_across_matrix_array(MatrixXf** matrix, int size, int element[]) {
VectorXf* vector = new VectorXf(size);
for (int i = 0; i < size; i++)
(*vector)[i] = (*matrix[i])(element,element);
return vector;
}
```

And I would need to write another function that returns a MatrixXf pointer if I wanted to access the ith rows of each matrix.

To delete these matrices I also need to pass into my delete function the sizes of the array. It’s a pain that it can’t dynamically discover the array size, but I don’t think there’s any way around it in C++.

```static void delete_4dMatrix(MatrixXf*** matrix, int sizes[]) {
for (int d1 = 0; d1 < sizes; d1++) {
for (int d2 = 0; d2 < sizes; d2++)
delete matrix[d1][d2];
}
}
```

Additionally, I need to have a separate function for initializing my 3D matrices, and any other dimension that I would need the way I have it set up now. Less than ideal. My guess is that there’s a elegant way to handle this with void pointers, but at the moment I only need 3D and 4D matrices so it’s not worth the effort for me to set it up.

So there are a few obstacles to using N-dimensional matrices with Eigen, but it’s really not a problem once you get a handle on it and how to properly reference different elements and matrices through your (N-2)-dimensional array of pointers.

Tagged , ,