How can I create an n-dimensional grid in numpy to evaluate a function for arbitrary n?

Question

I'm trying to create a naive numerical integration function to illustrate the benefits of Monte Carlo integration in high dimensions. I want something like this:

def quad_int(f, mins, maxs, numPoints=100):
    '''
    Use the naive (Riemann sum) method to numerically integrate f on a box 
    defined by the mins and maxs.

    INPUTS:
    f         - A function handle. Should accept a 1-D NumPy array 
        as input.
    mins      - A 1-D NumPy array of the minimum bounds on integration.
    maxs      - A 1-D NumPy array of the maximum bounds on integration.
    numPoints - An integer specifying the number of points to sample in 
        the Riemann-sum method. Defaults to 100.
    '''
    n = len(mins)

    # Create a grid of evenly spaced points to evaluate f on
    # Evaluate f at each point in the grid; sum all these values up
    dV = np.prod((maxs-mins/numPoints))
    # Multiply the sum by dV to get the approximate integral

I know my dV is going to cause problems with numerical stability, but right now what I'm having trouble with is creating the domain. If the number of dimensions was fixed, it would be easy enough to just use np.meshgrid like this:

# We don't want the last value since we are using left-hand Riemann sums
x = np.linspace(mins[0],maxs[0],numPoints)[:-1]
y = np.linspace(mins[1],maxs[1],numPoints)[:-1]
z = np.linspace(mins[2],maxs[2],numPoints)[:-1]

X, Y, Z = np.meshgrid(x,y,z)
tot = 0
for index, x in np.ndenumerate(X):
    tot += f(x, Y[index], Z[index])

Is there an analogue to np.meshgrid that can do this for arbitrary dimensions, maybe accept a tuple of arrays? Or is there some other way to do Riemann sums in higher dimensions? I've thought about doing it recursively but can't figure out how that would work.

Answer 1

You could use a list comprehension to generate all of the linspaces, and then pass that list to meshgrid with a * (to convert the list to a tuple of arguments).

XXX = np.meshgrid(*[np.linspace(i,j,numPoints)[:-1] for i,j in zip(mins,maxs)])

XXX is now a list of n arrays (each n dimensional).

I'm using straight forward Python list and argument operations.

np.lib.index_tricks has other index and grid generation functions and classes that might be of use. It's worth reading just to see how things can be done.

---

A neat trick used in various numpy functions when indexing arrays of unknown dimension, is to construct as list of the desired indices. It can include slice(None) where you'd normally see :. Then convert it to a tuple and use it.

In [606]: index=[2,3]
In [607]: [slice(None)]+index
Out[607]: [slice(None, None, None), 2, 3]
In [609]: Y[tuple([slice(None)]+index)]
Out[609]: array([ 0. ,  0.5,  1. ,  1.5])
In [611]: Y[:,2,3]
Out[611]: array([ 0. ,  0.5,  1. ,  1.5])

They use a list where they need to change elements. Converting to a tuple isn't always needed

index=[slice(None)]*3
index[1]=0
Y[index] # same as Y[:,0,:]