ENH: _lib: add common machinery for low-level callback functions

.gitignore

@@ -113,6 +113,7 @@ benchmarks/scipy
 benchmarks/html
 benchmarks/scipy-benchmarks
 scipy/__config__.py
+scipy/_lib/_ccallback_c.c
 scipy/cluster/_vq.c
 scipy/cluster/_hierarchy.c
 scipy/fftpack/_fftpackmodule.c

benchmarks/benchmarks/integrate.py

@@ -3,6 +3,27 @@
 import numpy as np
 from .common import Benchmark
 
+from scipy.integrate import quad
+
+try:
+    import ctypes
+    import scipy.integrate._test_multivariate as clib_test
+    from scipy._lib import _test_ccallback_cython
+except ImportError:
+    _test_ccallback_cython = None
+
+try:
+    from scipy import LowLevelCallable
+    from_cython = LowLevelCallable.from_cython
+except ImportError:
+    LowLevelCallable = lambda func, data: (func, data)
+    from_cython = lambda *a: a
+
+try:
+    import cffi
+except ImportError:
+    cffi = None
+
 try:
     from scipy.integrate import solve_bvp
 except ImportError:
@@ -61,3 +82,35 @@ def time_gas(self):
         y[0] = 0.5
         y[1] = -0.5
         solve_bvp(self.fun_gas, self.bc_gas, x, y, tol=self.TOL)
+
+
+class Quad(Benchmark):
+    def setup(self):
+        from math import sin
+
+        self.f_python = lambda x: sin(x)
+        self.f_cython = from_cython(_test_ccallback_cython, "sine")
+
+        lib = ctypes.CDLL(clib_test.__file__)
+
+        self.f_ctypes = lib._multivariate_sin
+        self.f_ctypes.restype = ctypes.c_double
+        self.f_ctypes.argtypes = (ctypes.c_int, ctypes.c_double)  # sic -- for backward compat
+
+        if cffi is not None:
+            voidp = ctypes.cast(self.f_ctypes, ctypes.c_void_p)
+            address = voidp.value
+            ffi = cffi.FFI()
+            self.f_cffi = LowLevelCallable(ffi.cast("double (*)(int, double *)", address))
+
+    def time_quad_python(self):
+        quad(self.f_python, 0, np.pi)
+
+    def time_quad_cython(self):
+        quad(self.f_cython, 0, np.pi)
+
+    def time_quad_ctypes(self):
+        quad(self.f_ctypes, 0, np.pi)
+
+    def time_quad_cffi(self):
+        quad(self.f_cffi, 0, np.pi)

doc/source/ccallback.rst

@@ -0,0 +1,30 @@
+============================
+Low-level callback functions
+============================
+
+.. currentmodule:: scipy
+
+Some functions in SciPy take as arguments callback functions, which
+can either be python callables or low-level compiled functions.  Using
+compiled callback functions can improve performance somewhat by
+avoiding wrapping data in Python objects.
+
+Such low-level functions in Scipy are wrapped in `LowLevelCallable`
+objects, which can be constructed from function pointers obtained from
+ctypes, cffi, Cython, or contained in Python `PyCapsule` objects.
+
+.. autosummary::
+   :toctree: generated/
+
+   LowLevelCallable
+
+.. seealso::
+
+   Functions accepting low-level callables:
+
+   `scipy.integrate.quad`, `scipy.ndimage.generic_filter`, `scipy.ndimage.generic_filter1d`,
+   `scipy.ndimage.geometric_transform`
+
+   Usage examples:
+
+   :ref:`ndimage-ccallbacks`, :ref:`quad-callbacks`

doc/source/index.rst

@@ -73,3 +73,4 @@ parameters available for the algorithms.
    special
    stats
    stats.mstats
+   ccallback

doc/source/tutorial/integrate.rst

@@ -312,40 +312,44 @@ does not correspond to
 
 because the order of the polynomial in f2 is larger than two.
 
-Faster integration using Ctypes
--------------------------------
-
-A user desiring reduced integration times may pass a C function pointer through
-`ctypes` to `quad`, `dblquad`, `tplquad` or `nquad` and it will be integrated
-and return a result in Python.  The performance increase here arises from two
-factors.  The primary improvement is faster function evaluation, which is
-provided by compilation.  This can also be achieved using a library like Cython
-or F2Py that compiles Python.  Additionally we have a speedup provided by the
-removal of function calls between C and Python in :obj:`quad` - this cannot be
-achieved through Cython or F2Py.  This method will provide a speed increase of
-~2x for trivial functions such as sine but can produce a much more noticeable
-increase (10x+) for more complex functions.  This feature then, is geared
-towards a user with numerically intensive integrations willing to write a
-little C to reduce computation time significantly.
-
-`ctypes` integration can be done in a few simple steps:
-
-1.) Write an integrand function in C with the function signature 
-``double f(int n, double args[n])``, where ``args`` is an array containing the
-arguments of the function f.  
+.. _quad-callbacks:
+
+Faster integration using low-level callback functions
+-----------------------------------------------------
+
+A user desiring reduced integration times may pass a C function
+pointer through `scipy.LowLevelCallable` to `quad`, `dblquad`,
+`tplquad` or `nquad` and it will be integrated and return a result in
+Python.  The performance increase here arises from two factors.  The
+primary improvement is faster function evaluation, which is provided
+by compilation of the function itself.  Additionally we have a speedup
+provided by the removal of function calls between C and Python in
+:obj:`quad`.  This method may provide a speed improvements of ~2x for
+trivial functions such as sine but can produce a much more noticeable
+improvements (10x+) for more complex functions.  This feature then, is
+geared towards a user with numerically intensive integrations willing
+to write a little C to reduce computation time significantly.
+
+The approach can be used, for example, via `ctypes` in a few simple steps:
+
+1.) Write an integrand function in C with the function signature
+``double f(int n, double *x, void *user_data)``, where ``x`` is an
+array containing the point the function f is evaluated at, and ``user_data``
+to arbitrary additional data you want to provide.
 
 .. code-block:: c
 
-   //testlib.c
-   double f(int n, double args[n]){
-       return args[0] - args[1] * args[2]; //corresponds to x0 - x1 * x2
+   /* testlib.c */
+   double f(int n, double *x, void *user_data) {
+       double c = *(double *)user_data;
+       return c + x[0] - x[1] * x[2]; /* corresponds to c + x - y * z */
    }
 
 2.) Now compile this file to a shared/dynamic library (a quick search will help
 with this as it is OS-dependent). The user must link any math libraries,
 etc. used.  On linux this looks like::
 
-    $ gcc -shared -o testlib.so -fPIC testlib.c
+    $ gcc -shared -fPIC -o testlib.so testlib.c
 
 The output library will be referred to as ``testlib.so``, but it may have a 
 different file extension. A library has now been created that can be loaded
@@ -354,23 +358,30 @@ into Python with `ctypes`.
 3.) Load shared library into Python using `ctypes` and set ``restypes`` and
 ``argtypes`` - this allows Scipy to interpret the function correctly:
 
->>> import ctypes
->>> from scipy import integrate
->>> lib = ctypes.CDLL('/**/testlib.so')   # Use absolute path to testlib
->>> func = lib.f   # Assign specific function to name func (for simplicity)
->>> func.restype = ctypes.c_double
->>> func.argtypes = (ctypes.c_int, ctypes.c_double)
+.. code:: python
+
+   import os, ctypes
+   from scipy import integrate, LowLevelCallable
+
+   lib = ctypes.CDLL(os.path.abspath('testlib.so'))
+   lib.f.restype = ctypes.c_double
+   lib.f.argtypes = (ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.c_void_p)
+
+   c = ctypes.c_double(1.0)
+   user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)
+
+   func = LowLevelCallable(lib.f, user_data)
 
-Note that the ``argtypes`` will always be ``(ctypes.c_int, ctypes.c_double)``
-regardless of the number of parameters, and ``restype`` will always be
-``ctypes.c_double``.
+The last ``void *user_data`` in the function is optional and can be omitted
+(both in the C function and ctypes argtypes) if not needed. Note that the
+coordinates are passed in as an array of doubles rather than a separate argument.
 
 4.) Now integrate the library function as normally, here using `nquad`:
 
 >>> integrate.nquad(func, [[0, 10], [-10, 0], [-1, 1]])
-(1000.0, 1.1102230246251565e-11)
+(1200.0, 1.1102230246251565e-11)
 
-And the Python tuple is returned as expected in a reduced amount of time.  All 
+The Python tuple is returned as expected in a reduced amount of time.  All 
 optional parameters can be used with this method including specifying
 singularities, infinite bounds, etc.