I am using mpmaths incomplete gamma function with negative z, and complex integration bounds a, b. See documentation: docs
This can't be done with scipy's incomplete gamma function. I benched marked mpmath against Mathematica and got the same results. Unfortunately the mpmaths function can only evaluate a single value at a time, which means I have to loop over the array which contains the values of z,a,b.
I believe vectorization is hard because for negative z one has to evaluate the integral recursively in parts (see this post and the respective wiki: incomplete gamma function in python?)
Does anyone know a module that can handle this? Or if there is no such thing if it is even possible?
So in essence:
from mpmath import gammainc
values = np.random.normal(0, 100, 1000) + 1j*np.random.normal(0, 100, 1000)
res_mp = np.array([gammainc(cc, 4, 10) for cc in values])
I've actually implemented something that I believe works to some extend:
import numpy as np
def new_quad_routine(func, c, a, b,x_list, w_list):
c_1 = (b-a)/2.0
c_2 = (b+a)/2.0
eval_points = c_1*x_list+c_2
func_evals = func(eval_points, c)
return c_1 * np.sum(func_evals * w_list[:,np.newaxis], axis=0)
def new_quad_gauss_7(func, c, a, b):
"""
integrates a complex function with bounds, a, b with an array of arguments c
call for instance with the gamma function:
new_quad_gauss_7(gamma_integrator, 4, 10)
"""
x_gauss = np.array([-0.949107912342759, -0.741531185599394, -0.405845151377397, 0, 0.405845151377397, 0.741531185599394, 0.949107912342759])
w_gauss = np.array([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830050505119, 0.279705391489277,0.129484966168870])
return new_quad_routine(func, c, a, b, x_gauss, w_gauss)
def gamma_integrator(t, c):
return t[:, np.newaxis]**c*np.exp(-t[:,np.newaxis])
def gammainc_function(c, a, b):
if type(c) is not np.ndarray:
raise ValueError("Need a numpy array for function argument!")
return new_quad_gauss_7(gamma_integrator, c-1, a, b)
This is based on the numerical integrator in this post.
It is pretty fast:
values=np.repeat(-2+3j-1, 100)
In [143]: %timeit new_quad_gauss_7(gamma_integrator, 4, 10, values)
107 µs ± 574 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [144]: %timeit [gammainc(cc, 4, 10) for cc in values]
224 ms ± 2.86 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
I don't know how good it is though. So the question slightly changes to how can I benchmark and or improve this? Furthermore, this seems to work for fixed bounds. However, I do not know how one would implement an integral with infinite bounds on either end.
