Simple example of PyBind11 landau binding

landau.cpp

// Taken from LCG ROOT MathLib
// License info:
// Authors: Andras Zsenei & Lorenzo Moneta   06/2005

/**********************************************************************
 *                                                                    *
 * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
 *                                                                    *
 *                                                                    *
 **********************************************************************/

// Langaus authors:
//  Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
//  Adapted for C++/ROOT by H.Pernegger (Heinz.Pernegger@cern.ch) and
//  Markus Friedl (Markus.Friedl@cern.ch)
//
//  Adaption for Python by David-Leon Pohl, pohl@physik.uni-bonn.de
//  Adaption for pybind11 by Henry Schreiner, henry.schreiner@cern.ch
//
// Can be manually compiled:
// c++ -O3 -Wall -shared -std=c++11 -fPIC `python -m pybind11 --includes` landau.cpp -o landau`python-config --extension-suffix`

#include <cmath>

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

namespace py = pybind11;
using namespace pybind11::literals;

const double invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)

double landauPDF(const double& x, const double& x0, const double& xi)  // same algorithm is used in GSL
{
	static double p1[5] =
	{ 0.4259894875, -0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253 };
	static double q1[5] =
	{ 1.0, -0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063 };

	static double p2[5] =
	{ 0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211 };
	static double q2[5] =
	{ 1.0, 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714 };

	static double p3[5] =
	{ 0.1788544503, 0.09359161662, 0.006325387654, 0.00006611667319, -0.000002031049101 };
	static double q3[5] =
	{ 1.0, 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675 };

	static double p4[5] =
	{ 0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186 };
	static double q4[5] =
	{ 1.0, 106.8615961, 337.6496214, 2016.712389, 1597.063511 };

	static double p5[5] =
	{ 1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910 };
	static double q5[5] =
	{ 1.0, 156.9424537, 3745.310488, 9834.698876, 66924.28357 };

	static double p6[5] =
	{ 1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109 };
	static double q6[5] =
	{ 1.0, 651.4101098, 56974.73333, 165917.4725, -2815759.939 };

	static double a1[3] =
	{ 0.04166666667, -0.01996527778, 0.02709538966 };

	static double a2[2] =
	{ -1.845568670, -4.284640743 };

	if (xi <= 0)
		return 0;
	double v = (x - x0) / xi;
	double u, ue, us, denlan;
	if (v < -5.5) {
		u = std::exp(v + 1.0);
		if (u < 1e-10)
			return 0.0;
		ue = std::exp(-1 / u);
		us = std::sqrt(u);
		denlan = 0.3989422803 * (ue / us) * (1 + (a1[0] + (a1[1] + a1[2] * u) * u) * u);
	}
	else if (v < -1) {
		u = std::exp(-v - 1);
		denlan = std::exp(-u) * std::sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) / (q1[0] + (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * v);
	}
	else if (v < 1) {
		denlan = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * v) / (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) * v);
	}
	else if (v < 5) {
		denlan = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * v) / (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) * v);
	}
	else if (v < 12) {
		u = 1 / v;
		denlan = u * u * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u) / (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] * u) * u) * u) * u);
	}
	else if (v < 50) {
		u = 1 / v;
		denlan = u * u * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u) / (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] * u) * u) * u) * u);
	}
	else if (v < 300) {
		u = 1 / v;
		denlan = u * u * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u) / (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] * u) * u) * u) * u);
	}
	else {
		u = 1 / (v - v * std::log(v) / (v + 1));
		denlan = u * u * (1 + (a2[0] + a2[1] * u) * u);
	}
	return denlan / xi;
}

double gaussPDF(const double& x, const double& mu, const double& sigma)
{
	return invsq2pi / sigma * exp(- pow((x - mu), 2) / (2 * pow(sigma, 2)));
}

double landauGaussPDF(const double& x, const double& mu, const double& eta, const double& sigma)
{
	// Numeric constants
	double mpshift = 0.; // -0.22278298;     // Landau maximum location shift in original code is wrong, since the shift does not depend on mu only

	// Control constants
	unsigned int np = 100;      // number of convolution steps
	const unsigned int sc = 8;        // convolution extends to +-sc Gaussian sigmas

	// Convolution steps have to be increased if sigma > eta * 5 to get stable solution that does not oscillate, addresses #1
	if (sigma > 3 * eta)
		np *= int(sigma / eta / 3.);
	if (np > 100000)  // Do not use too many convolution steps to save time
		np = 100000;

	// Variables
	double xx;
	double mpc;
	double fland;
	double sum = 0.0;
	double xlow, xupp;
	double step;

	// MP shift correction
	mpc = mu - mpshift;// * eta;

	// Range of convolution integral
	xlow = x - sc * sigma;
	xupp = x + sc * sigma;

	step = (xupp - xlow) / (double) np;

	// Discrete linear convolution of Landau and Gaussian
	for (unsigned int i = 1; i <= np / 2; ++i) {
		xx = xlow + double (i - 0.5) * step;
		fland = landauPDF(xx, mpc, eta) / eta;
		sum += fland * gaussPDF(x, xx, sigma);

		xx = xupp - double (i - 0.5) * step;
		fland = landauPDF(xx, mpc, eta) / eta;
		sum += fland * gaussPDF(x, xx, sigma);
	}

	const double norm = 0.398902310115109;  // normalization of the integral from -inf..intf

	return (step * sum);
}


PYBIND11_MODULE(landau, m) {
    m.def("landauPDF", py::vectorize(landauPDF),
            "x"_a, "x0"_a=0., "xi"_a=1.);

    m.def("gaussPDF", py::vectorize(gaussPDF),
            "x"_a, "mu"_a=0., "sigma"_a=1.);

    m.def("landauGaussPDF", py::vectorize(landauGaussPDF),
            "x"_a, "mu"_a=0., "eta"_a=1., "sigma"_a=1.);
}