Parallelize calculation of matrix elements using scipy.integrate.quad

Question

I have to evaluate every element of a matrix using a function with a numerical integral (scipy.integrate.quad). The elements of the matrix are pixels of a 5202x3465 gray image.
I have access to a GPU and I would like to evaluate as many elements as possible in parallel, because right now, with linear programming, the entire computation takes more than 24 hours.
Here it's the sample code:

for i in range(0, rows):
    for j in range(0, columns):
        img[i, j] = myFun(constant_args, i, j)

def myFunc(constant_args, i, j):
    new_pixel = quad(integrand, constant_args, i, j)
    ... other calculations ... 
    return new_pixel

I tried to use multiprocessing (as mp) like this:

arows = list(range(0, rows))
acolumns = list(range(0, columns))
with mp.Pool() as pool:
    img = pool.map(myFunc, (constant_args, arows, acolumns))

or with img = pool.map(myFunc(constant_args), (arows, acolumns))
but it gives me:
TypeError: myFunc() missing 2 required positional arguments: 'j' and 'i'

I don't understand how this works from other examples and I don't know the terminology used in the documentation.
I only want to divide that nested loop into subthreads, if someone suggests a different approach I'm all ears.
ps. I tried with numba but it gives errors when interacting with some Scipy libraries

Thank you in advance for your help!

Answer 1

First the error is in your call to the map- operation. It would have to be:

arows = list(range(0, rows))
acolumns = list(range(0, columns))
with mp.Pool() as pool:
    img = pool.map(myFunc, constant_args, arows, acolumns)

However, this might not result in what you are looking for as this just runs through the 3 arguments (which have to be lists). It is not running through combinations of them, especially arows and acolumns. For example, if constant_args has 3 elements, Pool.map would stop after 3 iterations without running through the longer lists arows or acolumns.

First you should make all combinations of rows and column indices

from itertools import product, repeat
comb = list(product(arows, acolumns))

That procudes something like (all possible combinations)

[(1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3), (3, 1), (3, 2), (3, 3)]

Next I would zip these combinations with your constant_args

constant_args = [10, 11]
arguments = list(zip(comb , repeat(constant_args)))

It produces a list of tuples with each tuple containing two elements. The first is your pixel position and the second are your constant_args

[((1, 1), [10, 11]),
 ((1, 2), [10, 11]),
 ((1, 3), [10, 11]),
 ((2, 1), [10, 11]),
 ((2, 2), [10, 11]),
 ((2, 3), [10, 11]),
 ((3, 1), [10, 11]),
 ((3, 2), [10, 11]),
 ((3, 3), [10, 11])]

Now we have to modify your myFuncsomewhat:

def myFunc(pix_id, constant_args):
    new_pixel = quad(integrand, constant_args, pix_id[0], pix_id[1])
    ... other calculations ... 
    return new_pixel

Finally, we use Pool.starmap to work the magic (see here: starmap usage):

with mp.Pool() as pool:
    img = pool.starmap(myFunc, arguments )

What happens is that starmap takes a list of tuples and provides these as inputs for the functions. However, starmap automatically unpacks the list of tuples into single arguments for your funtion. The first argument is pix_id consisting of two elements and the second argument is constant_args.

Answer 2

You can use quadpy (one of my projects). It does vectorized computation, so this will run very fast. Example with an output of shape 2x2:

import quadpy


def f(x):
    return [[x ** 2, x**3], [x**4, x**5]]


val, err = quadpy.quad(f, 0, 1)
print(val)

[[0.33333333 0.25      ]
 [0.2        0.16666667]]

The output of f must be of shape (..., x.shape), and the ... can be any tuple, e.g., (5202, 3465).