6

I wonder how to create a grid (multidimensional array) with numpy mgrid for an unknown number of dimensions (D), each dimension with a lower and upper bound and number of bins:

n_bins =  numpy.array([100 for  d in numpy.arrange(D)])
bounds = numpy.array([(0.,1) for d in numpy.arrange(D)])
grid = numpy.mgrid[numpy.linspace[(numpy.linspace(bounds(d)[0], bounds(d)[1], n_bins[d] for d in numpy.arrange(D)]

I guess above doesn't work, since mgrid creates array of indices not values. But how to use it to create array of values.

Thanks

Aso.agile

Aso Agile
  • 417
  • 1
  • 6
  • 13

1 Answers1

6

You might use

np.mgrid[[slice(row[0], row[1], n*1j) for row, n in zip(bounds, n_bins)]]

import numpy as np
D = 3
n_bins =  100*np.ones(D)
bounds = np.repeat([(0,1)], D, axis = 0)

result = np.mgrid[[slice(row[0], row[1], n*1j) for row, n in zip(bounds, n_bins)]]
ans = np.mgrid[0:1:100j,0:1:100j,0:1:100j]

assert np.allclose(result, ans)

Note that np.ogrid can be used in many places where np.mgrid is used, and it requires less memory because the arrays are smaller.

unutbu
  • 842,883
  • 184
  • 1,785
  • 1,677
  • Thanks @unutbu, it worked fine. Now how to convert this grid to D-dimensional points, i.e. 100**D x D shape? – Aso Agile Nov 24 '12 at 18:23
  • `reshape` indeed gives 100**D D-vectors, but surprisingly they are not unique, which also means not all possible vectors are there. What would be the way to get all D-points in the grid? I actually imagine it as a 100**D table, and have no idea what to do with the additional dimension `mgrid` gives. – Lev Levitsky Nov 25 '12 at 16:59
  • 1
    @LevLevitsky: To get points in the grid, take `result.reshape(D, -1).T`. The `-1` will be replaced by whatever number is necessary to use all the values in the array. Thanks for the question, this is probably more like what AsoAgile was looking for. – unutbu Nov 25 '12 at 19:33
  • @unutbu Thanks! This is exactly what I needed (not sure about the OP; I actually found the question on Google, was quite surprised to see it was so recent). What I'm wondering though is what is the actual meaning of the array returned by `mgrid`? What can it be used for without any transformations? – Lev Levitsky Nov 25 '12 at 19:48
  • 1
    @LevLevitsky: Suppose you have a function like `def f(x): return x[0]**2+x[1]`. Then `grid = np.mgrid[0:5,0:5]` is perfectly set up to evaluate `f(grid)`, which computes `f(x)` for every point on the grid. Note that `np.ogrid` returns smaller arrays which would work here too, by taking advantage of broadcasting. – unutbu Nov 25 '12 at 20:04
  • @unutbu Gosh... This is actually what I need, but I'd have **never** guessed it works this way. Even after reading the docs over for 2 days. `numpy` __really__ takes me some time to get used to. [instant edit] I was going to ask you to clarify this even further, but I think I'm starting to understand how it works now. Many thanks again! P.S. But I can't use it like this for **any** function, right? Only for algebraic stuff that handles arrays, like in your example. – Lev Levitsky Nov 25 '12 at 20:12
  • I guess I need `vectorize`; [this](http://stackoverflow.com/q/5197650/1258041) seems to answer my question. Sorry for spamming you :) Wow, and it's your answer I found. Small world. – Lev Levitsky Nov 25 '12 at 21:14
  • 1
    @LevLevitsky: It is true that you can not use this "trick" with any function. The computation must make sense in terms of numpy arrays. If your computation can not be expressed in terms of arrays, then you are stuck with using the much slower `vectorize` or `frompyfunc` functions. – unutbu Nov 25 '12 at 22:07