In Python everything is an object and so too in boost python. It is amazingly easy to pass a function defined in Python to the C++ layer.

In the following I set a file that can handle also numpy arrays. Let's first focus on the signature of the main function:

bp::list doSomething(np::ndarray input, PyObject *pyobj , PyObject *pyobj2)

This function wants 3 inputs: a numpy array, and 2 Python objects (functions). It returns a Python list. Looking deeply into the implementation of doSomething we can see what it actually does:

1. The code creates 2 empty arrays (zz and fzz) of the same shape as the input, one of the same type as the input one and a second with float members
2. The array is simply copied into zz the first array by using the pointer to the data of the zz boost numpy array
3. The second array contains sqrt values of the element of the input. Again the input is set by memcpy-ing the value in the data array of fzz

The 2 Python functions are used differently: we discard the return value of the first function, if there. The return of the second function is actually stored in one of the outputs of doSomething. Clearly the second function must return a number (a float) otherwise the whole thing won't work.

The interesting thing is now that you can set element-wise calculation at the C++ level with a function defined in Python!

import prova
import numpy as np
 
a = np.asarray([i for i in range(3*4*2)])
a = a.reshape([3,4,2])
print (a)
 
def print_element(input):
	print ("f: {0}".format(input))
 
prova.doSomething(a, print_element, None)
 
c = []
def append_to_list(input, shouldPrint=False):
	c.append(input)
	if shouldPrint:
		print ("{0} appended to list {1}".format(input, c))
 
def element_wise_algebra(input, shouldPrint=True):
	ret = input - 7
	if shouldPrint:
		print ("element_wise {0}".format(ret))
	return ret
 
prova.doSomething(a, append_to_list, None)
#print ("this is c: {0}".format(c))
 
b = prova.doSomething(a, None, element_wise_algebra)
#print (a)
print (b[5])

  #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
 
  #include <iostream>
  #include <cmath>
 
  #include <boost/python.hpp>
  #include <boost/python/numpy.hpp>
  #include "boost/tuple/tuple.hpp"
 
  #if defined(_WIN32) || defined(_WIN32) || defined(__WIN32__) || defined(_WIN64)
  #include <windows.h>
  // this trick only if compiler is MSVC
  __if_not_exists(uint8_t) { typedef __int8 uint8_t; }
  __if_not_exists(uint16_t) { typedef __int8 uint16_t; }
  #endif
 
  namespace bp = boost::python;
  namespace np = boost::python::numpy;
 
bp::list doSomething(np::ndarray input, PyObject *pyobj , PyObject *pyobj2) {
 
	boost::python::object output(boost::python::handle<>(boost::python::borrowed(pyobj)));
	int isOutput = !(output == boost::python::api::object());
 
	boost::python::object calculate(boost::python::handle<>(boost::python::borrowed(pyobj2)));
	int isCalculate = !(calculate == boost::python::api::object());
 
	int number_of_dims = input.get_nd();
	int dim_array[3];
 
	dim_array[0] = input.shape(0);
	dim_array[1] = input.shape(1);
	if (number_of_dims == 2) {
		dim_array[2] = -1;
	}
	else {
		dim_array[2] = input.shape(2);
	}
 
	/**************************************************************************/
	np::ndarray zz = zeros(3, dim_array, (int)0);
	np::ndarray fzz = zeros(3, dim_array, (float)0);
	/**************************************************************************/
 
	int * A = reinterpret_cast<int *>(input.get_data());
	int * B = reinterpret_cast<int *>(zz.get_data());
	float * C = reinterpret_cast<float *>(fzz.get_data());
 
	//Copy data and cast
	for (int i = 0; i < dim_array[0]; i++) {
		for (int j = 0; j < dim_array[1]; j++) {
			for (int k = 0; k < dim_array[2]; k++) {
				int index = k + dim_array[2] * j + dim_array[2] * dim_array[1] * i;
				int val = (*(A + index));
				float fval = sqrt((float)val);
				std::memcpy(B + index, &val, sizeof(int));
				std::memcpy(C + index, &fval, sizeof(float));
				// if the PyObj is not None evaluate the function 
				if (isOutput)	
					output(fval);
				if (isCalculate) {
					float nfval = (float)bp::extract<float>(calculate(val));
					if (isOutput)
						output(nfval);
					std::memcpy(C + index, &nfval, sizeof(float));
				}
			}
		}
	}
 
 
	bp::list result;
 
	result.append<int>(number_of_dims);
	result.append<int>(dim_array[0]);
	result.append<int>(dim_array[1]);
	result.append<int>(dim_array[2]);
	result.append<np::ndarray>(zz);
	result.append<np::ndarray>(fzz);
 
	return result;
 
}
 
 
BOOST_PYTHON_MODULE(prova)
{
	np::initialize();
 
	//To specify that this module is a package
	bp::object package = bp::scope();
	package.attr("__path__") = "prova";
 
	np::dtype dt1 = np::dtype::get_builtin<uint8_t>();
	np::dtype dt2 = np::dtype::get_builtin<uint16_t>();
 
	def("doSomething", doSomething);
}