ashkelly / pyquad

A lightweight, parallel python wrapper for the GSL integration library

pyquad

A drop-in replacment for scipy.integrate.quad which is much faster for repeat integrals over a parameter space.

The library provides a thin, parallel wrapper for the GNU Scientific Library (GSL) integration routines.

Examples

The code below is a python example of integrating over a grid of parameters using scipy.integrate.quad,

import numpy as np
import scipy.integrate


def test_integrand_func(x, alpha, beta, i, j, k, l):
    return x * alpha * beta + i * j * k


grid = np.random.random((10000000, 2))

res = np.zeros(grid.shape[0])
for i in range(res.shape[0]):
    res[i] = scipy.integrate.quad(
        test_integrand_func, 0, 1, (grid[i, 0], grid[i, 1], 1.0, 1.0, 1.0, 1.0)
    )[0]

this can be replaced with,

import numpy as np
import pyquad


def test_integrand_func(x, alpha, beta, i, j, k, l):
    return x * alpha * beta + i * j * k


grid = np.random.random((10000000, 2))

res, err = pyquad.quad_grid(test_integrand_func, 0, 1, grid, (1.0, 1.0, 1.0, 1.0))

which reduces the runtime significantly. For an example of the performance see the benchmarks below.

Benchmarks

We first compare a test integral in both pyquad and scipy,

and then we look in more detail at the scaling of pyquad with an increased thread count,

These benchmarks were carried out on cosma7 which has 28 cores (2x Intel Xeon Gold 5120 CPU @ 2.20GHz). Perfect scaling was never to be expected, given the problem becomes completely memorybound with a high core count.

Installing

To get started using the package you can use pip to download wheels for linux or osx,

pip install pyquad --user

or you can clone the repository,

git clone https://github.com/AshKelly/pyquad.git

and then go into the repository and run the setup file,

cd pyquad
python setup.py install

Requirements

The package requires that numpy is already installed and we require a C compiler to build from source.

Running the tests

The tests are currently incredibly primitive and just do a variety of answer testing by comparing to scipy.integrate.quad. These can be run with,

pytest tests

inside the pyquad folder. You will need to install pytest and scipy for this (pip install pytest scipy --user) i

History

I started to write components of this to help speed up some integrals for PyAutoLens. As it happened, a few other people in my department were also interested. I attempted to neaten up the API and roll it out for easy use.

If you want to contribute or want any extra feature implementing please just get in touch via email or a pull request.

Authors

• Ashley Kelly ([email protected])
• Arnau Quera ([email protected])

Citations

Please use the following citation:

@software{kelly:2020,
  author       = {Ashley J. Kelly},
  title        = {pyquad},
  month        = jul,
  year         = 2020,
  publisher    = {Zenodo},
  version      = {0.6.4},
  doi          = {10.5281/zenodo.3936959},
  url          = {https://doi.org/10.5281/zenodo.3936959}
}

License

This project is licensed under the GPL v3.0 License - see the LICENSE.md file for details

pyquad's Issues

Exceptionss that can be raised in Python

When I integrate a very elliptical mass profile, I get the following GSL error:

/home/jammy/PycharmProjects/VirtualEnvs/PyAuto/bin/python3 /home/jammy/PycharmProjects/PyAuto/PyAutoLens/test_autolens/numerics/deflections/sersic.py
gsl: pyquad/integration/qags.c:562: ERROR: integral is divergent, or slowly convergent
Default GSL error handler invoked.

Process finished with exit code 134 (interrupted by signal 6: SIGABRT)

Unfortunately, I can't bypass this exception using a try / except in Python, as the error is within pyquad and crashes the code.

Is there a neat way to stop this from crashing the code or returning a Python exception?

Sum inside an integral.

Hello,

It is possible passing and array as argument to a function? Something like that:

import numpy as np
import pyquad


def test_integrand_func(x, alpha, beta, i, j, k, l):
    a = np.sum(i)
    return x * alpha * beta + a * j * k



grid = np.random.random((10, 2))

i = np.array([1,2,3])

res, err = pyquad.quad_grid(test_integrand_func, 0, 1, grid, (i,1,1,1))

This is just a toy example. My problem actually is that i have a sum inside an integral, and i have to evaluate this integral over a grid of values (x,y). For now I'm using scipy.quad, but for a large grid the code is going very slow. Bellow I have pasted a more realistic example of I want to do.

import numpy as np
from scipy.integrate import quad


def fun(tau, x, y, alpha, beta):
    """
    x: int
    y: int
    alpha: array
    beta: array
    """
    x_til = x/alpha
    y_til = y/beta
    aux = (x_til**2 + y_til**2)
    exp = np.exp(- tau**2*aux)
    
    return (tau**2)*exp.sum()



grid = np.random.random((10, 2))

alpha = np.array([1,2,3])
beta = np.array([4,5,6])

res = np.empty_like(grid)

for i in range(len(grid)):
    res[i] = quad(fun, 0, 1, args=(grid[i][0],grid[i][1], alpha, beta))

Any help is welcome, and thanks in advance.

Cheers.