Integrating functions that return an array in Python

Question

I have a lot of data to integrate over and would like to find a way of doing it all with just matrices, and would be willing to compromise on accuracy for a performance boost. What I have in mind is something like this:

import numpy
import scipy

a = np.array([1,2,3])

def func(x):
    return x**2 + x

def func2(x):
    global a
    return a*x

def integrand(x):
    return func(x)*func2(x)

integrated = quad(integrand, 0, 1)

So I am trying to integrate each element in the array that comes out of integrand.

I'm aware that there is a possibility of using numpy.vectorize() like this:

integrated = numpy.vectorize(scipy.integrate.quad)(integrand, 0, 1)

but I can't get that working. Is there a way to do this in python?

Solution

Well now that I learnt a bit more python I can answer this question if anyone happens to stable upon it and has the same question. The way to do it is to write the functions as though they are going to take scalar values, and not vectors as inputs. So follow from my code above, what we would have is something like

import numpy as np
import scipy.integrate.quad

a = np.array([1, 2, 3]) # arbitrary array, can be any size

def func(x):
    return x**2 + x

def func2(x, a):
    return a*x

def integrand(x, a):
    return func(x)*func2(x, a)

def integrated(a):
    integrated, tmp = scipy.integrate.quad(integrand, 0, 1, args = (a))
    return integrated

def vectorizeInt():
    global a
    integrateArray = []
    for i in range(len(a)):
        integrate = integrated(a[i])
        integrateArray.append(integrate)
    return integrateArray

Not that the variable which you are integrating over must be the first input to the function. This is required for scipy.integrate.quad. If you are integrating over a method, it is the second argument after the typical self (i.e. x is integrated in def integrand(self, x, a):). Also the args = (a) is necessary to tell quad the value of a in the function integrand. If integrand has many arguments, say def integrand(x, a, b, c, d): you simply put the arguments in order in args. So that would be args = (a, b, c, d).

Answer 1

quadpy (a project of mine) is fully vectorized. Install with

pip install quadpy

and then do

import numpy
import quadpy


def integrand(x):
    return [numpy.sin(x), numpy.exp(x)]  # ,...


res, err = quadpy.quad(integrand, 0, 1)
print(res)
print(err)

[0.45969769 1.71828183]
[1.30995437e-20 1.14828375e-19]

Answer 2

vectorize won't provide help with improving the performance of code that uses quad. To use quad, you'll have to call it separately for each component of the value returned by integrate.

For a vectorized but less accurate approximation, you can use numpy.trapz or scipy.integrate.simps.

Your function definition (at least the one shown in the question) is implemented using numpy functions that all support broadcasting, so given a grid of x values on [0, 1], you can do this:

In [270]: x = np.linspace(0.0, 1.0, 9).reshape(-1,1)

In [271]: x
Out[271]: 
array([[ 0.   ],
       [ 0.125],
       [ 0.25 ],
       [ 0.375],
       [ 0.5  ],
       [ 0.625],
       [ 0.75 ],
       [ 0.875],
       [ 1.   ]])

In [272]: integrand(x)
Out[272]: 
array([[ 0.        ,  0.        ,  0.        ],
       [ 0.01757812,  0.03515625,  0.05273438],
       [ 0.078125  ,  0.15625   ,  0.234375  ],
       [ 0.19335938,  0.38671875,  0.58007812],
       [ 0.375     ,  0.75      ,  1.125     ],
       [ 0.63476562,  1.26953125,  1.90429688],
       [ 0.984375  ,  1.96875   ,  2.953125  ],
       [ 1.43554688,  2.87109375,  4.30664062],
       [ 2.        ,  4.        ,  6.        ]])

That is, by making x an array with shape (n, 1), the value returned by integrand(x) has shape (n, 3). There is one column for each value in a.

You can pass that value to numpy.trapz() or scipy.integrate.simps(), using axis=0, to get the three approximations of the integrals. You'll probably want a finer grid:

In [292]: x = np.linspace(0.0, 1.0, 101).reshape(-1,1)

In [293]: np.trapz(integrand(x), x, axis=0)
Out[293]: array([ 0.583375,  1.16675 ,  1.750125])

In [294]: simps(integrand(x), x, axis=0)
Out[294]: array([ 0.58333333,  1.16666667,  1.75      ])

Compare that to repeated calls to quad:

In [296]: np.array([quad(lambda t: integrand(t)[k], 0, 1)[0] for k in range(len(a))])
Out[296]: array([ 0.58333333,  1.16666667,  1.75      ])

Your function integrate (which I assume is just an example) is a cubic polynomial, for which Simpson's rule gives the exact result. In general, don't expect simps to give such an accurate answer.