Guvectorize does not understand how to broadcast

Numba Support: What is this error message?
1

Running this code

import numpy as np
from numba import guvectorize, int32, float32

@guvectorize([(int32, int32, int32[:], int32[:], int32, float32[:],
               float32[:, :], float32[:, :])],
             '(), (), (n), (n), (), (k), (l, p) -> (l, p)', nopython=True)
def calculate_and_insert_residuals(chunkposx, chunkposy, posx, posy, no_pixels,
                                   gaussian_parms, data, residual_data):
    peak = gaussian_parms[0]
    xbar = gaussian_parms[1] 
    ybar = gaussian_parms[2] 
    smaj = gaussian_parms[3]
    smin = gaussian_parms[4]
    theta = gaussian_parms[5]

    cos_theta = np.cos(theta)
    sin_theta = np.sin(theta)
    log2 = np.log(2.0)

    for i in range(posx.shape[0]):
        dx = posx[i] - xbar
        dy = posy[i] - ybar

        gauss_val = peak * np.exp(-log2 * (
            ((cos_theta * dx + sin_theta * dy) / smin) ** 2 +
            ((cos_theta * dy - sin_theta * dx) / smaj) ** 2))

        residual_data[posx[i] + chunkposx, posy[i] + chunkposy] += \
            data[posx[i] + chunkposx, posy[i] + chunkposy] - gauss_val

n_islands = 5 
n_positions = 10 

chunkposx = np.arange(0, 100, step=20, dtype=np.int32)
chunkposy = np.arange(0, 100, step=20, dtype=np.int32)
posx = np.random.randint(0, 20, size=(n_islands, n_positions), dtype=np.int32)
posy = np.random.randint(0, 20, size=(n_islands, n_positions), dtype=np.int32)
no_pixels = np.array([n_positions] * n_islands, dtype=np.int32)
gaussian_parms = 10 * np.random.rand(n_islands, 6).astype(np.float32)

data = np.zeros((100, 100), dtype=np.float32)
residual_data = np.zeros_like(data)

calculate_and_insert_residuals(chunkposx, chunkposy, posx, posy, no_pixels,
                               gaussian_parms, data, residual_data)

gives

ValueError: output operand requires a reduction along dimension -1, 
but the reduction is not enabled. 
The dimension size of 1 does not match the expected output shape.

while the equivalent Numpy code does work:

import numpy as np

def calculate_and_insert_residuals(chunkposx, chunkposy, posx, posy, no_pixels,
                                   gaussian_parms, data, residual_data):
    peak = gaussian_parms[0]
    xbar = gaussian_parms[1]
    ybar = gaussian_parms[2]
    smaj = gaussian_parms[3]
    smin = gaussian_parms[4]
    theta = gaussian_parms[5]

    cos_theta = np.cos(theta)
    sin_theta = np.sin(theta)
    log2 = np.log(2.0)

    for i in range(posx.shape[0]):
        dx = posx[i] - xbar
        dy = posy[i] - ybar

        gauss_val = peak * np.exp(-log2 * (
            ((cos_theta * dx + sin_theta * dy) / smin) ** 2 +
            ((cos_theta * dy - sin_theta * dx) / smaj) ** 2))

        residual_data[posx[i] + chunkposx, posy[i] + chunkposy] += \
            data[posx[i] + chunkposx, posy[i] + chunkposy] - gauss_val

for i in range(n_islands):
    calculate_and_insert_residuals(chunkposx[i], chunkposy[i], posx[i], posy[i], no_pixels[i],
                                   gaussian_parms[i], data, residual_data)

Why is that?

2

Could it be related to this?

“Broadcasting is not supported yet. Calling a guvectorize function in a scenario where broadcasting is needed may result in incorrect behavior. Numba will attempt to detect those cases and raise an exception.”
Creating NumPy universal functions — Numba 0+untagged.871.g53e976f.dirty documentation

3

Not sure…this is mentioned in the paragraph on Dynamic generalized universal functions which is not applicable here. It seems to refer to calling a function with a guvectorize decorator from a jitted function.

Besides that, the unique selling point of the guvectorize decorator is that it provides for broadcasting, so that should work.

4

Changing the signature from

to

'(), (), (n), (n), (), (k), (l, p), (l, p)'

i.e. no output, fixes this problem.