simple example on how to combine python and numpy with some fast c function using a c++ template function

gistfile1.txt

#include <Python.h>
#include <math.h>
#include <stdexcept>
#include <cstdio>
#include <numpy/arrayobject.h>


template<typename T>
void object_to_numpy1d_nocopy(T* &ptr, PyObject* obj, long long  &count, int& stride=stride_default, int type=NPY_DOUBLE) {
		if(obj == NULL)
			throw std::runtime_error("cannot convert to numpy array");
		if((int)PyArray_NDIM(obj) != 1)
			throw std::runtime_error("array is not 1d");
		long long size = PyArray_DIMS(obj)[0];
		if((count >= 0) && (size != count))
			throw std::runtime_error("arrays not of equal size");
		if(PyArray_TYPE(obj) != type)
			throw std::runtime_error("is not of proper type");
		npy_intp* strides =  PyArray_STRIDES(obj);
		if(stride == -1) {
			stride = strides[0];
		} else {
			if(strides[0] != stride*PyArray_ITEMSIZE(obj)) {
				throw Error("stride is not equal to %d", stride);
			}
		}
		
		ptr = (T*)PyArray_DATA(obj);
		count = size;
}



void histogram1d(const double* const __restrict__ block, const long long block_stride, const double* const weights, const int weights_stride, long long block_length, double* __restrict__ counts, const int counts_length, const double min, const double max) {
	
	const double scale = counts_length / (max-min);;
	for(long long i = 0; i < block_length; i++) {
		const double value = block[i]; //block[i*block_stride];
		if((value > min) & (value < max)) {
			const double scaled = (value - min) * scale;
			const long long index = (long long)(scaled);
			counts[index] += weights == NULL ? 1 : weights[i];
		}
	}
}




PyObject* histogram1d_(PyObject* self, PyObject* args) {
	PyObject* result = NULL;
	PyObject* block, *weights, *counts;
	double min, max;
	if(PyArg_ParseTuple(args, "OOOdd", &block, &weights, &counts, &min, &max)) {
		long long block_length = -1;
		long long counts_length = -1;
		double *block_ptr = NULL;
		int block_stride = -1;
		double *counts_ptr = NULL;
		
		double *weights_ptr = NULL;
		int weights_stride = -1;
		try {
			object_to_numpy1d_nocopy(block_ptr, block, block_length, block_stride);
			object_to_numpy1d_nocopy(counts_ptr, counts, counts_length, weights_stride);
			if(weights != Py_None) {
				object_to_numpy1d_nocopy(weights_ptr, weights, block_length);
			}
			Py_BEGIN_ALLOW_THREADS
			histogram1d(block_ptr, block_stride, weights_ptr, weights_stride, block_length, counts_ptr, counts_length, min, max);
			Py_END_ALLOW_THREADS
			Py_INCREF(Py_None);
			result = Py_None;
		} catch(std::runtime_error e) {
			PyErr_SetString(PyExc_RuntimeError, e.what());
		}
	}
	return result;
}

static PyMethodDef pygavi_functions[] = {
        {"histogram1d", (PyCFunction)histogram1d_, METH_VARARGS, ""},
    { NULL, NULL, 0 }
};

PyMODINIT_FUNC
initgavifast(void)
{
	import_array();
	Py_InitModule("gavifast", pygavi_functions);
}