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

Question

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?