0

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?

Mike
  • 21
  • 8
  • 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 Nov 05 '19 at 00: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 Nov 05 '19 at 00:48
  • That is, a different integral for the values of `r` and `c` (broadcasted or not)? – hpaulj Nov 05 '19 at 01: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 Nov 05 '19 at 01:08
  • I think you found the best link on the topic. I too found it while answering a similar question about `solve_ivp`, https://stackoverflow.com/questions/54991056/scipy-integrate-solve-ivp-vectorized. `scipy` does not have all the wiz-bang features that MATLAB provides. – hpaulj Nov 05 '19 at 01:14
  • Ok, thank you. Is there anything outside scipy that could help? I’d be gobsmacked if the python community hasn’t come up with a solution for this issue, as the difference in speed is phenomenal. – Mike Nov 05 '19 at 07:32

2 Answers2

1

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

ev-br
  • 24,968
  • 9
  • 65
  • 78
1

quadpy (a project of mine) has vectorized computation.

Nico Schlömer
  • 53,797
  • 27
  • 201
  • 249
  • 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? – Mike Nov 22 '19 at 11:08
  • Instead of `scipy.integrate.quad`, just use `quapy.quad`. – Nico Schlömer Nov 22 '19 at 14:26