pass_pycallback

.gitignore

_myclib.cpp
_myclib.so
build/
cython_debug/

00README.rst

pass_pycallback

This example shows how to pass python callbacks to an external C-library using Cython. It is assumed that the C-library allows the user to pass a void* to hold user state.

For demonstration purposes there is no C library, only the header file myclib.h containing a single function grad() (but you may imagine that you would link against an external library)

Status

Currently works.

Installation

To run the example you need Python, Cython and a C compiler.

$ python setup.py build_ext -i
$ python myclib_main.py

License

The code is in the public domain.

Author

Björn Dahlgren, January 2015.

README.rst

00README.rst

_myclib.pyx

# distutils: language = c++
# distutils: extra_compile_args = ['-O0', '-g']

import numpy as np
cimport numpy as cnp

from cpython.object cimport PyObject

from myclib_wrapper cimport py_grad

def grad(f, cnp.ndarray[cnp.float64_t, ndim=1] x, double h):
    cdef cnp.ndarray[cnp.float64_t, ndim=1] dfdx = np.zeros(x.size)
    py_grad(x.size, &x[0], h, &dfdx[0], <PyObject *>f)
    return dfdx

myclib.h

#ifdef __cplusplus
extern "C" {
#endif

typedef double (*RhsFn_Rn_R1)(int, const double * const, void *);

void grad(RhsFn_Rn_R1 f, int n, const double * const x, double h, 
          double * const out, void *user_data){
    // compute the gradient
    int i;
    double * xtmp = (double*)malloc(sizeof(double)*n);
    memcpy(xtmp, x, sizeof(double)*n);
    for (i=0; i<n; ++i){
        xtmp[i] += h;
        out[i] = (f(n, xtmp, user_data) - f(n, x, user_data))/h;
        xtmp[i] -= h;
    }
    free(xtmp);
}

#ifdef __cplusplus
}
#endif

myclib_main.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np

from _myclib import grad

def f_vector(x):
    """ f(x, y) = x**2/y + y**2
           dfdx = 2*x/y
           dfdy = -x**2/y**2 + 2*y
    """
    assert len(x) == 2
    return x[0]**2/x[1] + x[1]**2

def test_grad():
    x = np.array([3.0, 2.0])
    h = 1e-6
    result = grad(f_vector, x, h)
    assert x[0] == 3.0
    assert x[1] == 2.0
    ref = [
        (f_vector([x[0]+h, x[1]]) - f_vector(x)) / h,
        (f_vector([x[0], x[1]+h]) - f_vector(x)) / h
    ]

    dfdx = lambda x, y: 2*x/y
    dfdy = lambda x, y: -x**2/y**2 + 2*y

    assert abs(ref[0] - dfdx(*x)) < 1e-4
    assert abs(ref[1] - dfdy(*x)) < 1e-4
    assert abs(result[0] - ref[0]) < 1e-14
    assert abs(result[1] - ref[1]) < 1e-14


if __name__ == '__main__':
    test_grad()

myclib_wrapper.h

#include <Python.h>
#include <numpy/arrayobject.h>
#include "myclib.h"

double py_f(int n, const double * const x, void * user_data) {
    PyObject * f_cb = (PyObject *) user_data;
    npy_intp dims[1];
    dims[0] = n;
    PyObject * x_ = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, (void *)x);
    PyObject * args_ = Py_BuildValue("(O)", x_);
    PyObject * result_;
    if (PyCallable_Check(f_cb))
        result_ = PyEval_CallObject(f_cb, args_);
    else
        PyErr_SetString(PyExc_RuntimeError, "f() not callable");

    Py_DECREF(args_);
    Py_DECREF(x_);

    if (result_ == NULL)
        return 0.0/0.0; // Return NaN

    double result = PyFloat_AsDouble(result_);
    Py_DECREF(result_);
    return result;
}

void py_grad(int n, const double * const x, double h, 
             double * const out, PyObject * user_data)
{
    import_array(); // Enable PyArray_* functions
    grad(&py_f, n, x, h, out, (void *)user_data);
}

myclib_wrapper.pxd

from cpython.object cimport PyObject

cdef extern from "myclib_wrapper.h":
    void py_grad(int, const double * const, double, double * const, PyObject *)

setup.py

from distutils.core import setup
from Cython.Build import cythonize

if __name__ == '__main__':
    setup(
        name = 'myclib',
        ext_modules = cythonize("_myclib.pyx", gdb_debug=True)
    )