External C++ Models

Many libraries in the field of scientific computing are written in C++ to achieve fast code execution and do not have python bindings. The C++ model example shows how you can call C++ code from your python BaseModel class. It is primarily intended for fast implementations of the forward and jacobian method.

This example model uses the following features:

Calling C++ Code: Calls external c++ code using pybind11
Performance Comparison: The file eulerpi/examples/cpp/python_reference_plants.py includes python models implementing the same mapping. You can compare the performance of the different approaches.

Preparation

sudo apt install cmake
sudo apt install pybind11
sudo apt install eigen3-dev
sudo ln -s /usr/include/eigen3/Eigen /usr/include/Eigen

C++ Model Definition

#ifndef MY_CPP_MODEL
#define MY_CPP_MODEL

#include <cmath>

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

// Returned from py::array_t<double>.request()
// struct buffer_info {
//     void *ptr;
//     size_t itemsize;
//     std::string format;
//     int ndim;
//     std::vector<size_t> shape;
//     std::vector<size_t> strides;
// };

namespace py = pybind11;

// Mapping from 
py::array_t<double> forward(py::array_t<double> param) {
    py::buffer_info bufParam = param.request();

    // Storage for the result
    py::array_t<double> 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 */
        1,          /* How many dimensions? */
        { 3 },  /* Number of elements for each dimension */
        { sizeof(double) }  /* Strides for each dimension */
    ));
    py::buffer_info bufRes = result.request();

    double * ptrParam = static_cast<double *>(bufParam.ptr);
    double * ptrRes   = static_cast<double *>(bufRes.ptr);

    double water = ptrParam[0];
    double sun   = ptrParam[1];
    
    double size  = water * sun;
    double green = std::sin(M_PI * water) * std::sin(M_PI * sun);
    double flies = std::exp(water)-0.999;
    
    ptrRes[0] = size;
    ptrRes[1] = green;
    ptrRes[2] = flies;

    return result;
}

// If using Eigen::Ref, prevent copy of data by using RowMajoger storage order which is also used in numpy
// using RowMatrixXd = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
// https://pybind11.readthedocs.io/en/stable/advanced/cast/eigen.html

typedef Eigen::Matrix<double, 3, 2> Matrix32d;
//Jacobian for f:R^n -> R^m is Jf:R^m -> R^n = R^mxn
Matrix32d jacobian(Eigen::Vector2d param) {
    auto res = Matrix32d(3,2);
    res(0,0) = param(1);
    res(0,1) = param(0);
    res(1,0) = M_PI * std::cos(M_PI * param(0)) * std::sin(M_PI * param(1));
    res(1,1) = std::sin(M_PI * param(0)) * M_PI * std::cos(M_PI * param(1));
    res(2,0) = param(0) * std::exp(param(0));
    res(2,1) = 0.0;
    return res;
}

#endif

Wrapping the C++ Code

// pybind11_wrapper.cpp
#include <pybind11/pybind11.h>
#include "cpp_model.hpp"

PYBIND11_MODULE(cpp_model, m) {
    m.doc() = "C++ EPI Model implemented with pybind11"; // Optional module docstring
    m.def("forward", &forward, "A function that calculates the forward pass");
    m.def("jacobian", &jacobian, "A function that calculates the jacobian");
}

The example code is inconsistent in the following way: It uses a normal array for the forward method, but an eigen vector as input and an eigen matrix as output for the jacobian method. This allows to show us how to write wrapper code for normal arrays as well as for eigen objects. On the python side exclusively numpy 1d/2d arrays will be used.

Note

The C++ code is wrapped using pybind11. For more information on how to use pybind11 see: * PyBind11 Documentation: https://pybind11.readthedocs.io/en/stable/ * PyBind11 Eigen Notes: https://pybind11.readthedocs.io/en/stable/advanced/cast/eigen.html

Compilation

cd /eulerpi/examples/cpp/
mkdir build && cd build
cmake ..
make -j

Python Side Model

import numpy as np

from eulerpi.core.models import ArtificialModelInterface, BaseModel
from eulerpi.examples.cpp import cpp_model


class CppPlant(BaseModel, ArtificialModelInterface):
    """A plant model which uses a c++ library with eigen3 to evaluate the forward pass and the gradient
    Param0: Water [0,1]
    Param1: Sun   [0,1]
    Data0: Size [0,2] # the more water and sun the better
    Data1: Health [0,1], to much water is not good, too much sun is not good
    Data2: Sciarid ;)

    """

    param_dim = 2
    data_dim = 3

    PARAM_LIMITS = np.array([[0, 1], [0, 1]])
    CENTRAL_PARAM = np.array([0.5, 0.5])

    def __init__(
        self,
        central_param: np.ndarray = CENTRAL_PARAM,
        param_limits: np.ndarray = PARAM_LIMITS,
        name: str = None,
        **kwargs,
    ) -> None:
        super().__init__(
            central_param,
            param_limits,
            name,
            **kwargs,
        )

    def forward(self, param):
        return cpp_model.forward(param)

    def jacobian(self, param):
        return cpp_model.jacobian(param)

    def generate_artificial_params(self, num_samples: int):
        return np.random.rand(num_samples, 2)

You can use the cpp example model as template for your own C++ Model.