0

I'm converting some code from MATLAB to python, and am struggling to get a function that takes array arguments (as parameters) to integrate using Scipy.

I've reduced the code to a basic example that produces the same error in Scipy, while the equivalent MATLAB code functions as expected.

I'm trying to pass a row vector parametric argument of length M and a column vector parametric argument of length N to the function being integrated over a further, separate, integration argument, with the expectation that my integrated output will have the shape MxN.

The following python code produces this error:


File "C:\Anaconda3\lib\site-packages\scipy\integrate\quadrature.py", line 196, in quadrature
    if err < tol or err < rtol*abs(val):

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

from

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.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), vec_func=True)

The equivalent MATLAB code (which outputs an MxN shaped result) is

c = [1, 2];
r = transpose([2, 1]);

out = integral(@(p)intgd(p, r, c), -pi/2, pi/2, 'ArrayValued', true);

function intgd = intgd(p, r, c)

    c_bcr = repmat(c, length(r), 1);
    r_bcc = repmat(r, 1, length(c));

    A = ones(length(r), length(c));

    s = sin(r_bcc) - sin(p);
    A(s ~= 0) = sin(c_bcr(s ~= 0).*s(s ~= 0))./(c_bcr(s ~= 0).*s(s ~= 0));
    intgd = 1./fun(p, c).^2.*(A.^2);
end

function d = fun(p, c)

    p1 = zeros(1, length(c));
    mask = sqrt(pi./(2*c)) < 1;
    p1(mask) = acos(sqrt(pi./(2*c(mask))));

    d = zeros(1, length(c));
    mask = abs(p) <= p1; 
    d(mask) = 1./(pi./(2*c(mask).^2) + cos(p));
    mask = and(abs(p) > p1, abs(p) <= pi/2);
    d(mask) = 1./(pi./(2*c(mask).^2) + ((cos(p1(mask)) - cos(p))./2));

end

The MATLAB output out from the above is

[6.58727018139280,  0.963083280848789;
6.78600314283299,   1.05994693990888]

I'm unsure how scipy.integrate.quadrature is dealing with the dimensions of objects being passed through it, but the idea is that it should be producing the same MxN output.

I'm aware that numpy has its own built in broadcasting that generally avoids the need for explicit broadcasts as shown here, but I'm uncertain about how it deals with arrays where M=N, as in this case, so I've kept it explicit. Any pointers on this secondary issue would also be welcome.

Mike
  • 21
  • 8
  • `integrate.quadrature` expects a scalar function, not an array one. The `vec_func` just tells it that the function will accept multiple `p` values (and return a like sized array). – hpaulj Aug 05 '19 at 00:32
  • https://stackoverflow.com/questions/54991056/scipy-integrate-solve-ivp-vectorized, https://stackoverflow.com/questions/54893336/numerical-quadrature-of-scalar-valued-function-with-vector-input-using-scipy – hpaulj Aug 05 '19 at 01:30
  • Ok, thanks for the suggestions. I'll have a look and report back. Do you think ```nquad``` is what is needed, based on the description and code? – Mike Aug 05 '19 at 15:41
  • I've also reviewed the documentation for ```nquad``` based on the comment response to https://stackoverflow.com/questions/54893336/numerical-quadrature-of-scalar-valued-function-with-vector-input-using-scipy , but it seems from https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.nquad.html that this is not really what I'm looking for; ```nquad``` is a ```quad``` wrapper that allows integration over multiple integration arguments - it returns a ```float``` so that cannot produce the desired output. – Mike Aug 19 '19 at 20:47
  • I may not have outlined the question precisely. In my example, I want to integrate over the integration argument ```p``` while passing ```r``` and ```c``` as parametric (vector-valued) arguments that are not integrated over, but are simply passed into the function, with each value passed concurrently for each integration over ```p```. The output is then a vector with a size related to ```r``` and ```c```. This is what the MATLAB code achieves. – Mike Aug 19 '19 at 20:59
  • I've been looking into ```quadpy``` (https://github.com/nschloe/quadpy) as well but so far unable to understand how this may relate to my problem. – Mike Aug 20 '19 at 16:55

0 Answers0