(How) Can Scipy integrate functions with array valued arguments efficiently (without loops)?

Question

I want to integrate a function that takes array arguments using an efficient (vectorised/parallelised) method.

I can get the desired output using unwanted loops in the code as in the following (reduced complexity) example:

import numpy as np
from scipy import integrate

c = np.array([1, 2])
r = np.array([2, 1])

def fun(p, c):
    p1 = np.array(np.zeros(c.shape))
    mask = np.array(np.sqrt(np.pi/(2*c)) < 1)
    p1[mask] = np.arccos(np.sqrt(np.pi/(2*c[mask])))

    d = np.zeros(c.shape)
    mask = np.abs(p) <= p1
    d[mask] = 1/(np.pi/(2*c[mask]**2) + np.cos(p))
    mask = np.logical_and(np.abs(p) > p1, np.abs(p) <= np.pi/2)
    d[mask] = 1/(np.pi/(2*c[mask]**2) +
                     ((np.cos(p1[mask]) - np.cos(p))/2))
    return(d)

def intgd(p, r, c):
    A = np.ones((np.size(r), np.size(c)))

    s = np.sin(r) - np.sin(p)
    A[s != 0] = np.sin(c[s != 0]*s[s != 0])/(c[s != 0]*s[s != 0])

    return 1/fun(p, c)**2*(A**2)

res = np.zeros((np.size(r), np.size(c)))
for ii in range(0, np.size(r)):
    for jj in range(0, np.size(c)):
        res[ii, jj], err = integrate.quad(intgd, -np.pi/2, np.pi/2,
                                           epsabs=1e-10, limit=100,
                                           args=(r[ii], c[jj]))

However, my real function needs to handle much larger array inputs, which leads to unacceptably lengthy calculation duration.

I have tried variations along the following, and gained the knowledge that (as noted in comments on this question), the vec_func=True option for scipy.integrate.quadrature does not in fact enable vector-valued arguments to be passed as parameters into the function being integrated. [Aside: this makes it quite different to the MATLAB integral function, for which the option ArrayValued, true does enable that functionality, which results in much faster, apparently parallelised, evaluation of the integral.]

import numpy as np
from scipy import integrate

c = np.array([1, 2], ndmin=2)
r = np.array([2, 1])
r = r[:, np.newaxis]

def fun(p, c):
    p1 = np.zeros(c.shape)
    mask = np.array(np.sqrt(np.pi/(2*c)) < 1, ndmin=2)
    p1[mask] = np.arccos(np.sqrt(np.pi/(2*c[mask])))

    d = np.zeros(c.shape)
    mask = np.abs(p) <= p1
    d[mask] = 1/(np.pi/(2*c[mask]**2) + np.cos(p))
    mask = np.logical_and(np.abs(p) > p1, np.abs(p) <= np.pi/2)
    d[mask] = 1/(np.pi/(2*c[mask]**2) +
                     ((np.cos(p1[mask]) - np.cos(p))/2))
    return(d)

def intgd(p, r, c):
    A = np.ones((np.size(r), np.size(c)))

    c_bcr = np.broadcast_to(c, (np.size(r), np.size(c)))
    r_bcc = np.broadcast_to(r, (np.size(r), np.size(c)))

    s = np.sin(r_bcc) - np.sin(p)
    A[s != 0] = np.sin(c_bcr[s != 0]*s[s != 0])/(c_bcr[s != 0]*s[s != 0])

    return 1/fun(p, c)**2*(A**2)

res, err = integrate.quadrature(intgd, -np.pi/2, np.pi/2,
                                         args=(r, c), tol=1e-10, vec_func=True)

How can I use Scipy to integrate array-argument functions without resorting to loops?

Can you give a simple description of what the the 'loops' you are trying to avoid? What variable(s)? If I read quadrature right, it integrates a scalar valued function. vec_func=True means it can handle multiple values of the integration variable at once, returning a value for each. But it is still scalar integration, not multidimensional. nquad is for multidimensional integration. — hpaulj
– hpaulj, Commented Nov 5, 2019 at 0:43
Thanks. It’s looping over the values of the input array arguments I want to avoid. I think nquad is for multiple integrals where you want to evaluate more than one integration variable; that’s not the case here - there is only one integration variable: p, in the example. — Mike
– Mike, Commented Nov 5, 2019 at 0:48
That is, a different integral for the values of r and c (broadcasted or not)? — hpaulj
– hpaulj, Commented Nov 5, 2019 at 1:07
In the example, p is the variable integrated over while r and c are passed in as additional arguments. As r and c are arrays of values, they are looped through in a nested for loop to produce an output that has a shape of (np.size(r), np.size(c)) containing the result of the integration over p for each combination of r and c. I don’t want to have to loop through r and c, I want an integration function or method that recognises these variables as array-valued arguments and evaluates the integration over p for all combinations of r and c in an efficient, parallelised manner. — Mike
– Mike, Commented Nov 5, 2019 at 1:08
I think you found the best link on the topic. I too found it while answering a similar question about solve_ivp, stackoverflow.com/questions/54991056/…. scipy does not have all the wiz-bang features that MATLAB provides. — hpaulj
– hpaulj, Commented Nov 5, 2019 at 1:14

ev-br · Accepted Answer · 2019-11-05 13:33:22Z

1

Vectorized quad_vec will be available in scipy 1.4, when it's released.

answered Nov 5, 2019 at 13:33

ev-br

26.3k9 gold badges68 silver badges84 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

Mike Over a year ago

Thank you for the heads up, that will be very welcome!

Nico Schlömer · Accepted Answer · 2019-11-05 21:59:11Z

1

quadpy (a project of mine) has vectorized computation.

answered Nov 5, 2019 at 21:59

Nico Schlömer

59.7k35 gold badges216 silver badges291 bronze badges

2 Comments

Mike Over a year ago

Thank you for the suggestion @Nico , I have looked at quadpy and it looks like a really good library. For me, though, I think its learning curve may be a bit too steep, perhaps because I'm not really a mathematician. I found it difficult to see how to apply it - I can see one needs to have an understanding of the geometry of the solution domain but really this is beyond my capability. Are you able to suggest how it could be applied to the problem set out above?

Nico Schlömer Over a year ago

Instead of scipy.integrate.quad, just use quapy.quad.

Collectives™ on Stack Overflow

(How) Can Scipy integrate functions with array valued arguments efficiently (without loops)?

2 Answers 2

1 Comment

2 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

2 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related