The pyquad's discuss from ashkelly

View Code? Open in Web Editor NEW

A lightweight, parallel python wrapper for the GSL integration library

License: GNU General Public License v3.0

Python 5.75% C 93.57% C++ 0.68%

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