Foreign Function Interface



In [1]:

    
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')



In [2]:

    
import numpy as np

Wrapping functions written in C

Steps

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 setup.py to automate buiding of the Python extension module
Run python setup.py 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

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 cy_math.plus(a, 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 `setup.py`

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



In [7]:

    
%%file setup.py

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"],
                libraries=["m"],
                extra_compile_args=["-w",  "-std=c99"])

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









    



Writing setup.py

Building an extension module



In [8]:

    
! python setup.py clean
! python setup.py -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.cpython-35m-darwin.so 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))
print(cy_math.py_square(3))

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

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

C++

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]:

    
%%file setup.py 
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"])

setup(
    ext_modules = cythonize(
            ext,
            language="c++",        
      ))









    



Overwriting setup.py



In [17]:

    
%%bash
python setup.py -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

plus.plus(3, 4)









    



---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-18-4d2911b32d13> in <module>()
----> 1 import plus
      2 
      3 plus.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 https://github.com/JuliaLang/Rmath-julia.git
cd Rmath-julia/src
make
cd ../..

Functions to wrap



In [ ]:

    
! grep "\s.norm(" Rmath-julia/include/Rmath.h



In [ ]:

    
from ctypes import CDLL, c_int, c_double



In [ ]:

    
%%bash
ls Rmath-julia/src/*so



In [ ]:

    
lib = CDLL('Rmath-julia/src/libRmath-julia.so')

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)])
pass



In [ ]:

    
xs = np.linspace(-3,3,100)
plt.plot(xs, list(map(dnorm, xs)))
pass

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 [ ]:

    
%%file setup.py 
from distutils.core import setup, Extension
from Cython.Build import cythonize

ext = Extension("rmath",
                sources=["rmath.pyx"],
                include_dirs=["Rmath-julia/include"],
                library_dirs=["Rmath-julia/src"],
                libraries=["Rmath-julia"],
                runtime_library_dirs=["Rmath-julia/src"],
                extra_compile_args=["-w",  "-std=c99", "-DMATHLIB_STANDALONE"],
                extra_link_args=[],
               )

setup(
    ext_modules = cythonize(
            ext
    ))



In [ ]:

    
! python setup.py build_ext --inplace



In [ ]:

    
import rmath

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



In [ ]:

    
xs = np.linspace(-3,3,100)
plt.plot(xs, list(map(rmath.dnorm_, xs)))
pass

`Cython` wrappers are faster than `ctypes`



In [ ]:

    
%timeit pnorm(0, mu=2)
%timeit rmath.pnorm_(0, mu=2)

Fortran



In [ ]:

    
! pip install fortran-magic



In [ ]:

    
%load_ext fortranmagic



In [ ]:

    
%%fortran

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
end



In [ ]:

    
fort_sum(10)

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 [ ]:

    
%%R
library(ggplot2)
suppressPackageStartupMessages(
    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 [ ]:

    
mtcars.head(n=3)

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 [ ]: