Sending a C++ array to Python and back (Extending C++ with Numpy)

Question

I am going to send a c++ array to a python function as numpy array and get back another numpy array. After consulting with numpy documentation and some other threads and tweaking the code, finally the code is working but I would like to know if this code is written optimally considering the:

• Unnecessary copying of the array between c++ and numpy (python).
• Correct dereferencing of the variables.
• Easy straight-forward approach.

C++ code:

// python_embed.cpp : Defines the entry point for the console application. //  #include "stdafx.h"  #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include "Python.h" #include "numpy/arrayobject.h" #include<iostream>  using namespace std;  int _tmain(int argc, _TCHAR* argv[]) {     Py_SetProgramName(argv[0]);     Py_Initialize();     import_array()      // Build the 2D array     PyObject *pArgs, *pReturn, *pModule, *pFunc;     PyArrayObject *np_ret, *np_arg;     const int SIZE{ 10 };     npy_intp dims[2]{SIZE, SIZE};     const int ND{ 2 };     long double(*c_arr)[SIZE]{ new long double[SIZE][SIZE] };     long double* c_out;     for (int i{}; i < SIZE; i++)         for (int j{}; j < SIZE; j++)             c_arr[i][j] = i * SIZE + j;      np_arg = reinterpret_cast<PyArrayObject*>(PyArray_SimpleNewFromData(ND, dims, NPY_LONGDOUBLE,          reinterpret_cast<void*>(c_arr)));      // Calling array_tutorial from mymodule     PyObject *pName = PyUnicode_FromString("mymodule");     pModule = PyImport_Import(pName);     Py_DECREF(pName);     if (!pModule){         cout << "mymodule can not be imported" << endl;         Py_DECREF(np_arg);         delete[] c_arr;         return 1;     }     pFunc = PyObject_GetAttrString(pModule, "array_tutorial");     if (!pFunc || !PyCallable_Check(pFunc)){         Py_DECREF(pModule);         Py_XDECREF(pFunc);         Py_DECREF(np_arg);         delete[] c_arr;         cout << "array_tutorial is null or not callable" << endl;         return 1;     }     pArgs = PyTuple_New(1);     PyTuple_SetItem(pArgs, 0, reinterpret_cast<PyObject*>(np_arg));     pReturn = PyObject_CallObject(pFunc, pArgs);     np_ret = reinterpret_cast<PyArrayObject*>(pReturn);     if (PyArray_NDIM(np_ret) != ND - 1){ // row[0] is returned         cout << "Function returned with wrong dimension" << endl;         Py_DECREF(pFunc);         Py_DECREF(pModule);         Py_DECREF(np_arg);         Py_DECREF(np_ret);         delete[] c_arr;         return 1;     }     int len{ PyArray_SHAPE(np_ret)[0] };     c_out = reinterpret_cast<long double*>(PyArray_DATA(np_ret));     cout << "Printing output array" << endl;     for (int i{}; i < len; i++)         cout << c_out[i] << ' ';     cout << endl;      // Finalizing     Py_DECREF(pFunc);     Py_DECREF(pModule);     Py_DECREF(np_arg);     Py_DECREF(np_ret);     delete[] c_arr;     Py_Finalize();     return 0; } 

In CodeReview, there is a fantastic answer: Link...

Answer 1

Try out xtensor and the xtensor-python python bindings.

xtensor is a C++ library meant for numerical analysis with multi-dimensional array expressions.

xtensor provides

• an extensible expression system enabling numpy-style broadcasting (see the numpy to xtensor cheat sheet).
• an API following the idioms of the C++ standard library.
• tools to manipulate array expressions and build upon xtensor.
• bindings for Python, but also R and Julia.

Example of usage

Initialize a 2-D array and compute the sum of one of its rows and a 1-D array.

#include <iostream> #include "xtensor/xarray.hpp" #include "xtensor/xio.hpp"  xt::xarray<double> arr1   {{1.0, 2.0, 3.0},    {2.0, 5.0, 7.0},    {2.0, 5.0, 7.0}};  xt::xarray<double> arr2   {5.0, 6.0, 7.0};  xt::xarray<double> res = xt::view(arr1, 1) + arr2;  std::cout << res; 

Outputs

{7, 11, 14} 

Creating a Numpy-style universal function in C++.

#include "pybind11/pybind11.h" #include "xtensor-python/pyvectorize.hpp" #include <numeric> #include <cmath>  namespace py = pybind11;  double scalar_func(double i, double j) {     return std::sin(i) - std::cos(j); }  PYBIND11_PLUGIN(xtensor_python_test) {     py::module m("xtensor_python_test", "Test module for xtensor python bindings");      m.def("vectorized_func", xt::pyvectorize(scalar_func), "");      return m.ptr(); } 

Python code:

import numpy as np import xtensor_python_test as xt  x = np.arange(15).reshape(3, 5) y = [1, 2, 3, 4, 5] z = xt.vectorized_func(x, y) z 

Outputs

[[-0.540302,  1.257618,  1.89929 ,  0.794764, -1.040465],  [-1.499227,  0.136731,  1.646979,  1.643002,  0.128456],  [-1.084323, -0.583843,  0.45342 ,  1.073811,  0.706945]]