Using pybind11

The package pybind11 is provides an elegant way to wrap C++ code for Python, including automatic conversions for numpy arrays and the C++ Eigen linear algebra library. Used with the cppimport package, this provides a very nice work flow for integrating C++ and Python:

• Edit C++ code

• Run Python code

! pip install pybind11
! pip install cppimport

Clone the Eigen library if necessary - no installation is required as Eigen is a header only library.

! git clone https://github.com/RLovelett/eigen.git

Resources

• `pybind11 <http://pybind11.readthedocs.io/en/latest/>`__

• `cppimport <https://github.com/tbenthompson/cppimport>`__

• `Eigen <http://eigen.tuxfamily.org>`__

Example 1 - Basic usage

%%file ex1.cpp
<%
setup_pybind11(cfg)
%>
#include <pybind11/pybind11.h>
namespace py = pybind11;

PYBIND11_MODULE(ex1, m) {
    m.def("add", [](int a, int b) { return a + b; });
    m.def("mult", [](int a, int b) { return a * b; });
}

Overwriting ex1.cpp

import cppimport
ex1 = cppimport.imp("ex1")

ex1.add(3,4)

7

from ex1 import mult

mult(3,4)

12

ls ex1*so

ex10.cpython-36m-x86_64-linux-gnu.so*  ex1.cpython-36m-x86_64-linux-gnu.so*

Example 2 - Adding doc and named/default arguments

%%file ex2.cpp
<%
setup_pybind11(cfg)
%>
#include <pybind11/pybind11.h>
namespace py = pybind11;

using namespace pybind11::literals;

PYBIND11_MODULE(ex2, m) {
    m.def("add",
          [](int a, int b) { return a + b; },
          "Add two integers.",
          py::arg("a") = 3,
          py::arg("b") = 4);
    m.def("mult",
          [](int a, int b) { return a * b; },
          "Multiply two integers.",
          "a"_a=3,
          "b"_a=4);
}

Overwriting ex2.cpp

import cppimport
ex2 = cppimport.imp("ex2")

help(ex1.add)

Help on built-in function add in module ex1:

add(...) method of builtins.PyCapsule instance
    add(arg0: int, arg1: int) -> int

help(ex2.add)

Help on built-in function add in module ex2:

add(...) method of builtins.PyCapsule instance
    add(a: int = 3, b: int = 4) -> int

    Add two integers.

Example 3 - Split into execution modules for efficient compilation

%%file funcs1.cpp
#include <pybind11/pybind11.h>
namespace py = pybind11;

int add(int a, int b) {
    return a + b;
}

void init_f1(py::module &m) {
    m.def("add", &add);
}

Overwriting funcs1.cpp

%%file funcs2.cpp
#include <pybind11/pybind11.h>
namespace py = pybind11;

int mult(int a, int b) {
    return a + b;
}

void init_f2(py::module &m) {
    m.def("mult", &mult);
}

Overwriting funcs2.cpp

%%file ex3.cpp

<%
setup_pybind11(cfg)
cfg['sources'] = ['funcs1.cpp', 'funcs2.cpp']
%>
#include <pybind11/pybind11.h>
namespace py = pybind11;

void init_f1(py::module &m);
void init_f2(py::module &m);

PYBIND11_MODULE(ex3, m) {
    init_f1(m);
    init_f2(m);
}

Overwriting ex3.cpp

import cppimport
ex3 = cppimport.imp("ex3")

ex3.add(3,4), ex3.mult(3, 4)

(7, 7)

## Example 4 - Using setup.py to create shared libraries

%%file funcs.hpp
#pragma once

int add(int a, int b);
int mult(int a, int b);

Overwriting funcs.hpp

%%file funcs.cpp

#include "funcs.hpp"

int add(int a, int b) {
    return a + b;
}

int mult(int a, int b) {
    return a * b;
}

Overwriting funcs.cpp

%%file ex4.cpp
#include "funcs.hpp"
#include <pybind11/pybind11.h>
namespace py = pybind11;

PYBIND11_MODULE(ex4, m) {
    m.def("add", &add);
    m.def("mult", &mult);
}

Overwriting ex4.cpp

import os

if not os.path.exists('./pybind11'):
    ! git clone https://github.com/pybind/pybind11.git

%%file setup.py
import os, sys

from distutils.core import setup, Extension
from distutils import sysconfig

cpp_args = ['-std=c++14']

ext_modules = [
    Extension(
    'ex4',
        ['funcs.cpp', 'ex4.cpp'],
        include_dirs=['pybind11/include'],
    language='c++',
    extra_compile_args = cpp_args,
    ),
]

setup(
    name='ex4',
    version='0.0.1',
    author='Cliburn Chan',
    author_email='cliburn.chan@duke.edu',
    description='Example',
    ext_modules=ext_modules,
)

Overwriting setup.py

%%bash

python3 setup.py build_ext -i

running build_ext
building 'ex4' extension
gcc -pthread -B /opt/conda/compiler_compat -Wl,--sysroot=/ -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -Ipybind11/include -I/opt/conda/include/python3.6m -c funcs.cpp -o build/temp.linux-x86_64-3.6/funcs.o -std=c++14
gcc -pthread -B /opt/conda/compiler_compat -Wl,--sysroot=/ -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -Ipybind11/include -I/opt/conda/include/python3.6m -c ex4.cpp -o build/temp.linux-x86_64-3.6/ex4.o -std=c++14
g++ -pthread -shared -B /opt/conda/compiler_compat -L/opt/conda/lib -Wl,-rpath=/opt/conda/lib -Wl,--no-as-needed -Wl,--sysroot=/ build/temp.linux-x86_64-3.6/funcs.o build/temp.linux-x86_64-3.6/ex4.o -o /home/jovyan/work/sta-663-2020/notebooks/ex4.cpython-36m-x86_64-linux-gnu.so

cc1plus: warning: command line option ‘-Wstrict-prototypes’ is valid for C/ObjC but not for C++
cc1plus: warning: command line option ‘-Wstrict-prototypes’ is valid for C/ObjC but not for C++

import ex4

ex4.add(3,4), ex4.mult(3,4)

(7, 12)

Example 5 - Using STL containers

%%file ex5.cpp
<%
setup_pybind11(cfg)
%>
#include "funcs.hpp"
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

#include <vector>
namespace py = pybind11;

double vsum(const std::vector<double>& vs) {
    double res = 0;
    for (const auto& i: vs) {
        res += i;
    }
    return res;
}

std::vector<int> range(int start, int stop, int step) {
    std::vector<int> res;
    for (int i=start; i<stop; i+=step) {
        res.push_back(i);
    }
    return res;
}


PYBIND11_MODULE(ex5, m) {
    m.def("vsum", &vsum);
    m.def("range", &range);
}

Overwriting ex5.cpp

import cppimport
ex5 = cppimport.imp("ex5")

ex5.vsum(range(10))

45.0

ex5.range(1, 10, 2)

[1, 3, 5, 7, 9]

Using cppimport

The cppimport package allows you to specify several options. See Github page

Use of cppimport.imp

Note that cppimport.imp only needs to be called to build the shared library. Once it is called, the shared library is created and can be sued. Any updates to the C++ files will be detected by cppimport and it will automatically trigger a re-build.

Example 6: Vectorizing functions for use with numpy arrays

Example showing how to vectorize a square function. Note that from here on, we don’t bother to use separate header and implementation files for these code snippets, and just write them together with the wrapping code in a code.cpp file. This means that with cppimport, there are only two files that we actually code for, a C++ code.cpp file and a python test file.

%%file ex6.cpp
<%
cfg['compiler_args'] = ['-std=c++14']
setup_pybind11(cfg)
%>

#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

namespace py = pybind11;
double square(double x) {
    return x * x;
}

PYBIND11_MODULE(ex6, m) {
    m.doc() = "pybind11 example plugin";
    m.def("square", py::vectorize(square), "A vectroized square function.");
}

Overwriting ex6.cpp

import cppimport

ex6 = cppimport.imp("ex6")
ex6.square([1,2,3])

array([1., 4., 9.])

import ex6

ex6.square([2,4,6])

array([ 4., 16., 36.])

Example 7: Using numpy arrays as function arguments and return values

Example showing how to pass numpy arrays in and out of functions. These numpy array arguments can either be generic py:array or typed py:array_t<double>. The properties of the numpy array can be obtained by calling its request method. This returns a struct of the following form:

struct buffer_info {
    void *ptr;
    size_t itemsize;
    std::string format;
    int ndim;
    std::vector<size_t> shape;
    std::vector<size_t> strides;
};

Here is C++ code for two functions - the function twice shows how to change a passed in numpy array in-place using pointers; the function sum shows how to sum the elements of a numpy array. By taking advantage of the information in buffer_info, the code will work for arbitrary n-d arrays.

%%file ex7.cpp
<%
cfg['compiler_args'] = ['-std=c++11']
setup_pybind11(cfg)
%>

#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

namespace py = pybind11;

// Passing in an array of doubles
void twice(py::array_t<double> xs) {
    py::buffer_info info = xs.request();
    auto ptr = static_cast<double *>(info.ptr);

    int n = 1;
    for (auto r: info.shape) {
      n *= r;
    }

    for (int i = 0; i <n; i++) {
        *ptr++ *= 2;
    }
}

// Passing in a generic array
double sum(py::array xs) {
    py::buffer_info info = xs.request();
    auto ptr = static_cast<double *>(info.ptr);

    int n = 1;
    for (auto r: info.shape) {
      n *= r;
    }

    double s = 0.0;
    for (int i = 0; i <n; i++) {
        s += *ptr++;
    }

    return s;
}

PYBIND11_MODULE(ex7, m) {
    m.doc() = "auto-compiled c++ extension";
    m.def("sum", &sum);
    m.def("twice", &twice);
}

Overwriting ex7.cpp

%%file test_code.py
import cppimport
import numpy as np

ex7 = cppimport.imp("ex7")

if __name__ == '__main__':
    xs = np.arange(12).reshape(3,4).astype('float')
    print(xs)
    print("np :", xs.sum())
    print("cpp:", ex7.sum(xs))

    print()
    ex7.twice(xs)
    print(xs)

Overwriting test_code.py

%%bash

python test_code.py

[[ 0.  1.  2.  3.]
 [ 4.  5.  6.  7.]
 [ 8.  9. 10. 11.]]
np : 66.0
cpp: 66.0

[[ 0.  2.  4.  6.]
 [ 8. 10. 12. 14.]
 [16. 18. 20. 22.]]

Example 8: More on working with numpy arrays (optional)

This example shows how to use array access for numpy arrays within the C++ function. It is taken from the pybind11 documentation, but fixes a small bug in the official version. As noted in the documentation, the function would be more easily coded using py::vectorize.

%%file ex8.cpp
<%
cfg['compiler_args'] = ['-std=c++11']
setup_pybind11(cfg)
%>

#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

namespace py = pybind11;

py::array_t<double> add_arrays(py::array_t<double> input1, py::array_t<double> input2) {
    auto buf1 = input1.request(), buf2 = input2.request();

    if (buf1.ndim != 1 || buf2.ndim != 1)
        throw std::runtime_error("Number of dimensions must be one");

    if (buf1.shape[0] != buf2.shape[0])
        throw std::runtime_error("Input shapes must match");

    auto result = py::array(py::buffer_info(
        nullptr,            /* Pointer to data (nullptr -> ask NumPy to allocate!) */
        sizeof(double),     /* Size of one item */
        py::format_descriptor<double>::value, /* Buffer format */
        buf1.ndim,          /* How many dimensions? */
        { buf1.shape[0] },  /* Number of elements for each dimension */
        { sizeof(double) }  /* Strides for each dimension */
    ));

    auto buf3 = result.request();

    double *ptr1 = (double *) buf1.ptr,
           *ptr2 = (double *) buf2.ptr,
           *ptr3 = (double *) buf3.ptr;

    for (size_t idx = 0; idx < buf1.shape[0]; idx++)
        ptr3[idx] = ptr1[idx] + ptr2[idx];

    return result;
}

PYBIND11_MODULE(ex8, m) {
    m.def("add_arrays", &add_arrays, "Add two NumPy arrays");
}

Overwriting ex8.cpp

import cppimport
import numpy as np

code = cppimport.imp("ex8")

xs = np.arange(12)
print(xs)

print(code.add_arrays(xs, xs))

[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0.  2.  4.  6.  8. 10. 12. 14. 16. 18. 20. 22.]

Example 9: Using the C++ eigen library to calculate matrix inverse and determinant

Example showing how Eigen vectors and matrices can be passed in and out of C++ functions. Note that Eigen arrays are automatically converted to/from numpy arrays simply by including the pybind/eigen.h header. Because of this, it is probably simplest in most cases to work with Eigen vectors and matrices rather than py::buffer or py::array where py::vectorize is insufficient.

Note: When working with matrices, you can make code using eigen more efficient by ensuring that the eigen Matrix and numpy array have the same data types and storage layout, and using the Eigen::Ref class to pass in arguments. By default, numpy stores data in row major format while Eigen stores data in column major format, and this incompatibility triggers a copy which can be expensive for large matrices. There are basically 3 ways to make pass by reference work:

1. Use Eigen reference with arbitrary storage order

Eigen::Ref<MatrixType, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>> 2. Use Eigen row order matrices

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> 3. Create numpy arrays with column order

np.array(data, order='F')

This is an advanced topic that you can explore in the docs.

%%file ex9.cpp
<%
cfg['compiler_args'] = ['-std=c++11']
cfg['include_dirs'] = ['./eigen']
setup_pybind11(cfg)
%>

#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>

#include <Eigen/LU>

namespace py = pybind11;

// convenient matrix indexing comes for free
double get(Eigen::MatrixXd xs, int i, int j) {
    return xs(i, j);
}

// takes numpy array as input and returns double
double det(Eigen::MatrixXd xs) {
    return xs.determinant();
}

// takes numpy array as input and returns another numpy array
Eigen::MatrixXd inv(Eigen::MatrixXd xs) {
    return xs.inverse();
}

PYBIND11_MODULE(ex9, m) {
    m.doc() = "auto-compiled c++ extension";
    m.def("inv", &inv);
    m.def("det", &det);
}

Overwriting ex9.cpp

import cppimport
import numpy as np

code = cppimport.imp("ex9")

A = np.array([[1,2,1],
              [2,1,0],
              [-1,1,2]])

print(A)
print(code.det(A))
print(code.inv(A))

[[ 1  2  1]
 [ 2  1  0]
 [-1  1  2]]
-3.0
[[-0.66666667  1.          0.33333333]
 [ 1.33333333 -1.         -0.66666667]
 [-1.          1.          1.        ]]

Example 10: Using pybind11 with openmp (optional)

Here is an example of using OpenMP to integrate the value of \(\pi\) written using pybind11.

%%file ex10.cpp
/*
<%
cfg['compiler_args'] = ['-std=c++11', '-fopenmp']
cfg['linker_args'] = ['-lgomp']
setup_pybind11(cfg)
%>
*/
#include <cmath>
#include <omp.h>
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

namespace py = pybind11;

// Passing in an array of doubles
void twice(py::array_t<double> xs) {
    py::gil_scoped_acquire acquire;

    py::buffer_info info = xs.request();
    auto ptr = static_cast<double *>(info.ptr);

    int n = 1;
    for (auto r: info.shape) {
      n *= r;
    }

    #pragma omp parallel for
    for (int i = 0; i <n; i++) {
        *ptr++ *= 2;
    }
}

PYBIND11_MODULE(ex10, m) {
  m.doc() = "auto-compiled c++ extension";
  m.def("twice", [](py::array_t<double> xs) {
      /* Release GIL before calling into C++ code */
      py::gil_scoped_release release;
      return twice(xs);
    });
}

Overwriting ex10.cpp

import cppimport
import numpy as np

code = cppimport.imp("ex10")
xs = np.arange(10).astype('double')
code.twice(xs)
xs

array([0., 2., 2., 3., 4., 5., 6., 7., 8., 9.])