Foreign Function Interface

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt'ggplot')

In [2]:
import numpy as np

Wrapping functions written in C


  • Write the C header and implementation files
  • Write the Cython .pxd file to declare C function signatures
  • Write the Cython .pyx file to wrap the C functions for Python
  • Write to automate buiding of the Python extension module
  • Run python build_ext --inplace to build the module
  • Import module in Python like any other Python module

C header file

In [3]:
%%file c_math.h

#pragma once
double plus(double a, double b);
double mult(double a, double b);
double square(double a);
double acc(double *xs, int size);

Writing c_math.h

C implementation file

In [4]:
%%file c_math.c
#include <math.h>
#include "c_math.h"

double plus(double a, double b) {
    return a + b;

double mult(double a, double b) {
    return a * b;

double square(double a) {
    return pow(a, 2);

double acc(double *xs, int size) {
    double s = 0;
    for (int i=0; i<size; i++) {
        s += xs[i];
    return s;

Writing c_math.c

Cython "header" file

The .pxd file is similar to a header file for Cython. In other words, we can cimport <filename>.pxd in the regular Cython .pyx files to get access to functions declared in the .pxd files.

In [5]:
%%file cy_math.pxd

cdef extern from "c_math.h":
    double plus(double a, double b)
    double mult(double a, double b)
    double square(double a)
    double acc(double *xs, int size)

Writing cy_math.pxd

Cython "implementation" file

Here is whhere we actually wrap the C code for use in Python. Note especially how we handle passing in of arrays to a C function expecting a pointer to double using typed memoryviews.

In [6]:
%%file cy_math.pyx

cimport cy_math

def py_plus(double a, double b):
    return, b)

def py_mult(double a, double b):
    return cy_math.mult(a, b)

def py_square(double a):
    return cy_math.square(a)

def py_sum(double[::1] xs):
    cdef int size = len(xs)
    return cy_math.acc(&xs[0], size)

Writing cy_math.pyx

Build script

This is a build script for Python, similar to a Makefile

In [7]:

from distutils.core import setup, Extension
from Cython.Build import cythonize
import numpy as np

ext = Extension("cy_math",
                sources=["cy_math.pyx", "c_math.c"],
                extra_compile_args=["-w",  "-std=c99"])

setup(name = "Math Funcs",
      ext_modules = cythonize(ext))


Building an extension module

In [8]:
! python clean
! python -q build_ext --inplace

Compiling cy_math.pyx because it changed.
[1/1] Cythonizing cy_math.pyx
running clean

In [9]:
! ls cy_math*

cy_math.c                     cy_math.pxd cy_math.pyx

Using the extension module in Python

In [10]:
import cy_math
import numpy as np

print(cy_math.py_plus(3, 4))
print(cy_math.py_mult(3, 4))

xs = np.arange(10, dtype='float')


Confirm that we are getting C speedups by comparing with pure Python accumulator

In [11]:
def acc(xs):
    s = 0
    for x in xs:
        s += x
    return s

In [12]:
import cy_math

xs = np.arange(1000000, dtype='float')
%timeit -r3 -n3 acc(xs)
%timeit -r3 -n3 cy_math.py_sum(xs)

3 loops, best of 3: 197 ms per loop
3 loops, best of 3: 1.95 ms per loop


This is similar to C. We will use Cython to wrap a simple funciton.

In [13]:
%%file add.hpp 
#pragma once
int add(int a, int b);

Writing add.hpp

In [14]:
%%file add.cpp
int add(int a, int b) {
    return a+b;

Writing add.cpp

In [15]:
%%file plus.pyx

cdef extern from 'add.cpp':
    int add(int a, int b)
def plus(a, b):
    return add(a, b)

Writing plus.pyx

Note that essentially the only difference from C is language="C++" and the flag -std=c++11

In [16]:
from distutils.core import setup, Extension
from Cython.Build import cythonize

ext = Extension("plus",
                sources=["plus.pyx", "add.cpp"],
                extra_compile_args=["-w", "-std=c++11"])

    ext_modules = cythonize(


In [17]:
python -q build_ext --inplace

Please put "# distutils: language=c++" in your .pyx or .pxd file(s)
Compiling plus.pyx because it changed.
[1/1] Cythonizing plus.pyx
error: invalid argument '-std=c++11' not allowed with 'C/ObjC'
error: command '/usr/bin/clang' failed with exit status 1

In [18]:
import plus, 4)

ImportError                               Traceback (most recent call last)
<ipython-input-18-4d2911b32d13> in <module>()
----> 1 import plus
      3, 4)

ImportError: No module named 'plus'

Wrap an R function from libRmath using ctypes

R comes with a standalone C library of special functions and distributions, as described in the official documentation. These functions can be wrapped for use in Python.

Building the Rmath standalone library

git clone
cd Rmath-julia/src
cd ../..

Functions to wrap

In [ ]:
! grep "\s.norm(" Rmath-julia/include/Rmath.h

In [ ]:
from ctypes import CDLL, c_int, c_double

In [ ]:
ls Rmath-julia/src/*so

In [ ]:
lib = CDLL('Rmath-julia/src/')

def rnorm(mu=0, sigma=1):
    lib.rnorm.argtypes = [c_double, c_double]
    lib.rnorm.restype  = c_double
    return lib.rnorm(mu, sigma)

def dnorm(x, mean=0, sd=1, log=0):
    lib.dnorm4.argtypes = [c_double, c_double, c_double, c_int]
    lib.dnorm4.restype  = c_double
    return lib.dnorm4(x, mean, sd, log)

def pnorm(q, mu=0, sd=1, lower_tail=1, log_p=0):
    lib.pnorm5.argtypes = [c_double, c_double, c_double, c_int, c_int]
    lib.pnorm5.restype  = c_double
    return lib.pnorm5(q, mu, sd, lower_tail, log_p)

def qnorm(p, mu=0, sd=1, lower_tail=1, log_p=0):
    lib.qnorm5.argtypes = [c_double, c_double, c_double, c_int, c_int]
    lib.qnorm5.restype  = c_double
    return lib.qnorm5(p, mu, sd, lower_tail, log_p)

In [ ]:
pnorm(0, mu=2)

In [ ]:
qnorm(0.022750131948179212, mu=2)

In [ ]:
plt.hist([rnorm() for i in range(100)])

In [ ]:
xs = np.linspace(-3,3,100)
plt.plot(xs, list(map(dnorm, xs)))

Using Cython to wrap standalone library

In [ ]:
%%file rmath.pxd

cdef extern from "Rmath-julia/include/Rmath.h":
    double dnorm(double, double, double, int)
    double pnorm(double, double, double, int, int)
    double qnorm(double, double, double, int, int)
    double rnorm(double, double)

In [ ]:
%%file rmath.pyx

cimport rmath

def rnorm_(mu=0, sigma=1):
    return rmath.rnorm(mu, sigma)

def dnorm_(x, mean=0, sd=1, log=0):
    return rmath.dnorm(x, mean, sd, log)

def pnorm_(q, mu=0, sd=1, lower_tail=1, log_p=0):
    return rmath.pnorm(q, mu, sd, lower_tail, log_p)

def qnorm_(p, mu=0, sd=1, lower_tail=1, log_p=0):
    return rmath.qnorm(p, mu, sd, lower_tail, log_p)

In [ ]:
from distutils.core import setup, Extension
from Cython.Build import cythonize

ext = Extension("rmath",
                extra_compile_args=["-w",  "-std=c99", "-DMATHLIB_STANDALONE"],

    ext_modules = cythonize(

In [ ]:
! python build_ext --inplace

In [ ]:
import rmath

plt.hist([rmath.rnorm_() for i in range(100)])

In [ ]:
xs = np.linspace(-3,3,100)
plt.plot(xs, list(map(rmath.dnorm_, xs)))

Cython wrappers are faster than ctypes

In [ ]:
%timeit pnorm(0, mu=2)
%timeit rmath.pnorm_(0, mu=2)


In [ ]:
! pip install fortran-magic

In [ ]:
%load_ext fortranmagic

In [ ]:

subroutine fort_sum(N, s)
    integer*8, intent(in) :: N
    integer*8, intent(out) :: s
    integer*8 i
    s = 0
    do i = 1, N
        s = s + i*i
    end do

In [ ]:

Another example from the documentation

In [ ]:
%%fortran --link lapack

subroutine solve(A, b, x, n)
    ! solve the matrix equation A*x=b using LAPACK
    implicit none

    real*8, dimension(n,n), intent(in) :: A
    real*8, dimension(n), intent(in) :: b
    real*8, dimension(n), intent(out) :: x

    integer :: pivot(n), ok

    integer, intent(in) :: n
    x = b

    ! find the solution using the LAPACK routine SGESV
    call DGESV(n, 1, A, n, pivot, x, n, ok)
end subroutine

In [ ]:
A = np.array([[1, 2.5], [-3, 4]])
b = np.array([1, 2.5])

solve(A, b)

Interfacing with R

In [ ]:
%load_ext rpy2.ipython

In [ ]:
    ggplot(mtcars, aes(x=wt, y=mpg)) + geom_point() + geom_smooth(method=loess)

Converting between Python and R

In [ ]:
%R -o mtcars

mtcars is now a Python dataframe

In [ ]:

We can also pass data from Python to R

In [ ]:
x = np.linspace(0, 2*np.pi, 100)
y = np.sin(x)

In [ ]:
%%R -i x,y
plot(x, y, main="Sine curve in R base graphics")

In [ ]: