3
$\begingroup$

I need to calculate the integral of this function

def f(z):
    return ((1-2*z)*np.exp(-d/z))/(((1-z)**(2+d))*(z**(2-d)))

Here d is a constant. I am using this function:

 import scipy.integrate as integrate
 return integrate.quad(f,0, 5)

However, I need to calculate this for many different d values. I tried putting all d values in a numpy array but it didn't work. Any idea how to vectorize this?

$\endgroup$
2
$\begingroup$

The function integrate.quad is a python wrapper to the DQAGSE function from QUADPACK. This function uses adaptive quadrature, i.e. it will apply a fixed rule (in this case Gauss-Kronrod) on intervals that it will adaptively refine trying to reach the absolute and/or relative tolerance you requested.

Since the parameter d alters the behaviour of the function, it will also change the iterative interval refinement. Hence, it doesn't make sense to "vectorize".

What you could do is use the Pool class from the multiprocessing module in the standard Python library. Distribute your values of d over multiple cores. Wrap your integral in a function:

import numpy as np
import scipy.integrate as integrate
from multiprocessing import Pool

def f(z,d):
    return ((1-2*z)*np.exp(-d/z))/(((1-z)**(2+d))*(z**(2-d)))

def F_int(d):
    return integrate.quad(f,0.0,0.5,args=(d))


if __name__ == '__main__':
    d_array = np.array([1, 0.1, 2.0])
    pool = Pool(processes=3)
    print(pool.map(F_int, d_array))

The only restriction with this solution is that it doesn't work in interactive modes (so not in Jupyter notebooks nor Qtconsole). See here.

$\endgroup$
  • $\begingroup$ Thanks, I would have to read up on this in detail, but would this be able to handle large sizes of d_array...of the order of a million? $\endgroup$ – dayum Oct 31 '16 at 19:58
  • 1
    $\begingroup$ I just compared a serial version to the parallel version using 4 processes for a vector d with length 1.000.000 (generated randomly but with a fixed seed so the serial and parallel run get the same values): approx 345s vs 101s (on a quad core machine). $\endgroup$ – GertVdE Oct 31 '16 at 20:52
  • $\begingroup$ Got it, thanks. I am currently trying to run it in a module that is not called "main" so I changed if name__=='__main' to =="mymodule" . Also, I removed processes=3 since I dont know how many processes are needed. I'm getting a runtime error saying "Attempt to start a new process before the current process has finished its bootstrapping phase. This probably means that you are on Windows and you have forgotten to use the proper idiom in the main module: if name == 'main': freeze_support() $\endgroup$ – dayum Nov 1 '16 at 0:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.