# Return a 2D matrix from a C++ function in Python: ctypes versus pybind11

I found two ways in which I can use a custom C++ function (returning a 2D matrix) in Python. One approach that uses the Python module ctypes and another that uses the Python module pybind11.

## Approach 1 using ctypes:

I followed this tutorial (e.g., https://stackoverflow.com/questions/64084033/modern-2020-way-to-call-c-code-from-python) to create a my C++ code, a shared library from the C++ code and then a Python module.

#include <iostream>
#include <cmath>
#include <vector>
#include <cstdio>
using namespace std;

int from_xy(int x, int y, int nside) {
return x + (nside * y);
}

int modNegOperator(int k, int n){
return ((k %= n) < 0) ? k+n : k;
}

std::vector<int> to_xy(int k, int nside) {
int x = modNegOperator(k, nside);
int y = floor(k / nside);
vector<int> res(x, y);
return res;
}

extern "C"
double** make2DMatrix(unsigned nside,  double mx, double my) {

int n_matrix = nside * nside;
double **array2D = 0;
array2D = new double *[n_matrix];

for (int h = 0; h < n_matrix; h++) {
array2D[h] = new double[n_matrix];

for (int w = 0; w < n_matrix; w++) {
// fill in some initial values
// (filling in zeros would be more logic, but this is just for the example)
array2D[h][w] = 0;
}
}

for (int h = 0; h < n_matrix; h++){

std::vector<int> xy_vec = to_xy(h, nside);
int modneg1 = modNegOperator(xy_vec[0] + 1,nside);
int modneg2 = modNegOperator(xy_vec[0] - 1,nside);
int modneg3 = modNegOperator(xy_vec[1] + 1,nside);
int modneg4 = modNegOperator(xy_vec[1] - 1,nside);

int pos1 = from_xy(modneg1, xy_vec[1], nside);
int pos2 = from_xy(modneg2, xy_vec[1], nside);
int pos3 = from_xy(xy_vec[0], modneg3, nside);
int pos4 = from_xy(xy_vec[0], modneg4, nside);

double half_mx = mx / 2;
double half_my = my / 2;

array2D[h][h] = 0;
array2D[h][pos1] = half_mx;
array2D[h][pos2] = half_mx;

array2D[h][pos3] = half_my;
array2D[h][pos4] = half_my;

}

return array2D;
}


I have then created a shared library:

g++ -fPIC --std=c++11 -shared -o libTest.so myfunc.cpp


Created a myLib.py file including:

import ctypes as ct
import sys
import os

dir_path = os.path.dirname(os.path.realpath(__file__))
handle = ct.CDLL(dir_path + "/libTest.so")

handle.make2DMatrix.argtypes = [ctypes.c_uint, ctypes.c_double, ctypes.c_double]
handle.make2DMatrix.restype = ct.POINTER(ct.POINTER(ct.c_double))

def make2DMatrix(nside, mx, my):
return handle.make2DMatrix(nside, mx, my)


which is then loaded as a Python module:

import myLib

mat = myLib.make2DMatrix(4, 0.01, 0.01)


## Approach 2 using pybind11:

This is the edited C++ file, which is adapted for being used in pybind11. This includes the std::vector type, the #include <pybind11/pybind11.h> and #include <pybind11/stl.h>, and

PYBIND11_MODULE(make2DMatrix, m) {
m.doc() = "pybind11 example plugin"; // optional module docstring

m.def("make2DMatrix", &make2DMatrix, "A function that computes a matrix");
}


So, the C++ code to be used with pybind11, looks like:

#include <iostream>
#include <cmath>
#include <vector>
#include <cstdio>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

using namespace std;
namespace py = pybind11;

int from_xy(int x, int y, int nside) {
return x + (nside * y);
}

std::vector<int> to_xy(int k, int nside) {
int x = k%nside;
int y = floor(k / nside);
vector<int> res{x, y};
return res;
}

int modNegOperator(int k, int n){
return ((k %= n) < 0) ? k+n : k;
}

std::vector<std::vector<double>> make2DMatrix(int nside,  double mx, double my) {

int n_matrix = nside * nside;
std::vector<std::vector<double>> array2D(n_matrix, std::vector<double>(n_matrix, 0));

for (int h = 0; h < n_matrix; h++){

std::vector<int> xy_vec = to_xy(h, nside);
int modneg1 = modNegOperator(xy_vec[0] + 1,nside);
int modneg2 = modNegOperator(xy_vec[0] - 1,nside);
int modneg3 = modNegOperator(xy_vec[1] + 1,nside);
int modneg4 = modNegOperator(xy_vec[1] - 1,nside);

int pos1 = from_xy(modneg1, xy_vec[1], nside);
int pos2 = from_xy(modneg2, xy_vec[1], nside);
int pos3 = from_xy(xy_vec[0], modneg3, nside);
int pos4 = from_xy(xy_vec[0], modneg4, nside);

double half_mx = mx / 2;
double half_my = my / 2;

array2D[h][h] = 0;
array2D[h][pos1] = half_mx;
array2D[h][pos2] = half_mx;

array2D[h][pos3] = half_my;
array2D[h][pos4] = half_my;
}

return array2D;
}

PYBIND11_MODULE(make2DMatrix, m) {
m.doc() = "pybind11 example plugin"; // optional module docstring

m.def("make2DMatrix", &make2DMatrix, "A function that computes a matrix");
}


The C++ file is then compiled:

c++ -O3 -Wall -shared -std=c++11 -fPIC $(python3 -m pybind11 --includes) make2DMatrix.cpp -o make2DMatrix$(python3-config --extension-suffix)


And once I open Python, I import the make2DMatrix module:

import make2DMatrix

make2DMatrix.make2DMatrix(4, 0.01, 0.01)


Do you think I could improve both codes for speed? Do you have any other suggestion?