I have a Python function of two arguments, f(x,y), which returns a scalar. I also have two np.arrays of possible arguments, xs and ys. I'd like to compute a two-dimensional numpy array of shape (xs.size, ys.size) which holds values of f evaluated at all combinations of arguments drawn from xs and ys. The result should be equal to
np.array([[f(x,y) for y in ys] for x in xs])
I would like to achieve this as efficiently as possible and take advantage from multiple processor cores. f does not have any side-effect. I tried to use numpy.vectorize but did not manage to use it to the desired effect. E.g. numpy.outer does exactly what I want for the special case of f=operator.__mul__.
As an additional caveat, though potentially not immediately relevant, one of the array, say ys, does not contain numbers but objects as returned by scipy.interpolate.interp1d. The function f calculates a definite integral via scipy.integrate.quad of these interpolations. The resulting matrix still only contains numbers.