Create a n-D grid; meshgrid alternative

Hey,
quite often I want to evaluate a function on a grid in a n-D space. What I end up doing (and what I really dislike) looks something like this:
x = np.linspace(0, 5, 20) M1, M2 = np.meshgrid(x, x) X = np.column_stack([M1.flatten(), M2.flatten()]) X.shape # (400, 2)
fancy_function(X)
I don't think I ever used `meshgrid` in any other way. Is there a better way to create such a grid space?
I wrote myself a little helper function:
def gridspace(linspaces): return np.column_stack([space.flatten() for space in np.meshgrid(*linspaces)])
But maybe something like this should be part of numpy?
Best, Stefan

I just drafted different versions of the `gridspace` function: https://tmp23.tmpnb.org/user/1waoqQ8PJBJ7/notebooks/2015-05%20gridspace.ipyn...
Beste Grüße, Stefan
On Sun, May 10, 2015 at 1:40 PM, Stefan Otte stefan.otte@gmail.com wrote:
Hey,
quite often I want to evaluate a function on a grid in a n-D space. What I end up doing (and what I really dislike) looks something like this:
x = np.linspace(0, 5, 20) M1, M2 = np.meshgrid(x, x) X = np.column_stack([M1.flatten(), M2.flatten()]) X.shape # (400, 2)
fancy_function(X)
I don't think I ever used `meshgrid` in any other way. Is there a better way to create such a grid space?
I wrote myself a little helper function:
def gridspace(linspaces): return np.column_stack([space.flatten() for space in np.meshgrid(*linspaces)])
But maybe something like this should be part of numpy?
Best, Stefan

On Sun, May 10, 2015 at 7:05 AM, Stefan Otte stefan.otte@gmail.com wrote:
I just drafted different versions of the `gridspace` function:
https://tmp23.tmpnb.org/user/1waoqQ8PJBJ7/notebooks/2015-05%20gridspace.ipyn...
The link seems to be broken...
Jaime

On Sun, May 10, 2015 at 4:40 AM, Stefan Otte stefan.otte@gmail.com wrote:
Hey,
quite often I want to evaluate a function on a grid in a n-D space. What I end up doing (and what I really dislike) looks something like this:
x = np.linspace(0, 5, 20) M1, M2 = np.meshgrid(x, x) X = np.column_stack([M1.flatten(), M2.flatten()]) X.shape # (400, 2)
fancy_function(X)
I don't think I ever used `meshgrid` in any other way. Is there a better way to create such a grid space?
I wrote myself a little helper function:
def gridspace(linspaces): return np.column_stack([space.flatten() for space in np.meshgrid(*linspaces)])
But maybe something like this should be part of numpy?
Isn't what you are trying to build a cartesian product function? There is a neat, efficient implementation of such a function in StackOverflow, by our own pv.:
http://stackoverflow.com/questions/1208118/using-numpy-to-build-an-array-of-...
Perhaps we could make this part of numpy.lib.arraysetops? Isthere room for other combinatoric generators, i.e. combinations, permutations... as in itertools?
Jaime

On Sun, May 10, 2015 at 4:40 AM, Stefan Otte stefan.otte@gmail.com wrote:
Hey,
quite often I want to evaluate a function on a grid in a n-D space. What I end up doing (and what I really dislike) looks something like this:
x = np.linspace(0, 5, 20) M1, M2 = np.meshgrid(x, x) X = np.column_stack([M1.flatten(), M2.flatten()]) X.shape # (400, 2)
fancy_function(X)
I don't think I ever used `meshgrid` in any other way. Is there a better way to create such a grid space?
I feel like our "house style" has moved away from automatic flattening, and would maybe we should be nudging people towards something more like
# using proposed np.stack from pull request #5605 X = np.stack(np.meshgrid(x, x), axis=-1) assert X.shape == (20, 20, 2) fancy_function(X) # vectorized to accept any array with shape (..., 2)
-n

On 2015-05-10 14:46:12, Jaime Fernández del Río jaime.frio@gmail.com wrote:
Isn't what you are trying to build a cartesian product function? There is a neat, efficient implementation of such a function in StackOverflow, by our own pv.:
http://stackoverflow.com/questions/1208118/using-numpy-to-build-an-array-of-...
And a slightly faster version just down that page ;)
Stéfan

Hello,
indeed I was looking for the cartesian product.
I timed the two stackoverflow answers and the winner is not quite as clear:
n_elements: 10 cartesian 0.00427 cartesian2 0.00172 n_elements: 100 cartesian 0.02758 cartesian2 0.01044 n_elements: 1000 cartesian 0.97628 cartesian2 1.12145 n_elements: 5000 cartesian 17.14133 cartesian2 31.12241
(This is for two arrays as parameters: np.linspace(0, 1, n_elements)) cartesian2 seems to be slower for bigger.
I'd really appreciate if this was be part of numpy. Should I create a pull request?
Regarding combinations and permutations: I could be convenient to have as well.
Cheers, Stefan

I'm totally in favor of the 'gridspace(linspaces)' version, as you probably end up wanting to create grids of other things than linspaces (e.g. a logspace grid, or a grid of random points etc.).
It should be called somewhat different though. Maybe 'cartesian(arrays)'?
Best, Johannes
Quoting Stefan Otte (2015-05-10 16:05:02)
I just drafted different versions of the `gridspace` function: https://tmp23.tmpnb.org/user/1waoqQ8PJBJ7/notebooks/2015-05%20gridspace.ipyn...
Beste Grüße, Stefan
On Sun, May 10, 2015 at 1:40 PM, Stefan Otte stefan.otte@gmail.com wrote:
Hey,
quite often I want to evaluate a function on a grid in a n-D space. What I end up doing (and what I really dislike) looks something like this:
x = np.linspace(0, 5, 20) M1, M2 = np.meshgrid(x, x) X = np.column_stack([M1.flatten(), M2.flatten()]) X.shape # (400, 2)
fancy_function(X)
I don't think I ever used `meshgrid` in any other way. Is there a better way to create such a grid space?
I wrote myself a little helper function:
def gridspace(linspaces): return np.column_stack([space.flatten() for space in np.meshgrid(*linspaces)])
But maybe something like this should be part of numpy?
Best, Stefan
NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

On Tue, May 12, 2015 at 1:17 AM, Stefan Otte stefan.otte@gmail.com wrote:
Hello,
indeed I was looking for the cartesian product.
I timed the two stackoverflow answers and the winner is not quite as clear:
n_elements: 10 cartesian 0.00427 cartesian2 0.00172 n_elements: 100 cartesian 0.02758 cartesian2 0.01044 n_elements: 1000 cartesian 0.97628 cartesian2 1.12145 n_elements: 5000 cartesian 17.14133 cartesian2 31.12241
(This is for two arrays as parameters: np.linspace(0, 1, n_elements)) cartesian2 seems to be slower for bigger.
On my system, the following variation on Pauli's answer is 2-4x faster than his for your test cases:
def cartesian4(arrays, out=None): arrays = [np.asarray(x).ravel() for x in arrays] dtype = np.result_type(*arrays)
n = np.prod([arr.size for arr in arrays]) if out is None: out = np.empty((len(arrays), n), dtype=dtype) else: out = out.T
for j, arr in enumerate(arrays): n /= arr.size out.shape = (len(arrays), -1, arr.size, n) out[j] = arr[np.newaxis, :, np.newaxis] out.shape = (len(arrays), -1)
return out.T
I'd really appreciate if this was be part of numpy. Should I create a pull request?
There hasn't been any opposition, quite the contrary, so yes, I would go ahead an create that PR. I somehow feel this belongs with the set operations, rather than with the indexing ones. Other thoughts?
Also for consideration: should it work on flattened arrays? or should we give it an axis argument, and then "broadcast on the rest", a la generalized ufunc?
Jaime

Hey,
here is an ipython notebook with benchmarks of all implementations (scroll to the bottom for plots):
https://github.com/sotte/ipynb_snippets/blob/master/2015-05%20gridspace%20-%...
Overall, Jaime's version is the fastest.
On Tue, May 12, 2015 at 2:01 PM Jaime Fernández del Río < jaime.frio@gmail.com> wrote:
On Tue, May 12, 2015 at 1:17 AM, Stefan Otte stefan.otte@gmail.com wrote:
Hello,
indeed I was looking for the cartesian product.
I timed the two stackoverflow answers and the winner is not quite as clear:
n_elements: 10 cartesian 0.00427 cartesian2 0.00172 n_elements: 100 cartesian 0.02758 cartesian2 0.01044 n_elements: 1000 cartesian 0.97628 cartesian2 1.12145 n_elements: 5000 cartesian 17.14133 cartesian2 31.12241
(This is for two arrays as parameters: np.linspace(0, 1, n_elements)) cartesian2 seems to be slower for bigger.
On my system, the following variation on Pauli's answer is 2-4x faster than his for your test cases:
def cartesian4(arrays, out=None): arrays = [np.asarray(x).ravel() for x in arrays] dtype = np.result_type(*arrays)
n = np.prod([arr.size for arr in arrays]) if out is None: out = np.empty((len(arrays), n), dtype=dtype) else: out = out.T for j, arr in enumerate(arrays): n /= arr.size out.shape = (len(arrays), -1, arr.size, n) out[j] = arr[np.newaxis, :, np.newaxis] out.shape = (len(arrays), -1) return out.T
I'd really appreciate if this was be part of numpy. Should I create a pull request?
There hasn't been any opposition, quite the contrary, so yes, I would go ahead an create that PR. I somehow feel this belongs with the set operations, rather than with the indexing ones. Other thoughts?
Also for consideration: should it work on flattened arrays? or should we give it an axis argument, and then "broadcast on the rest", a la generalized ufunc?
Jaime
-- (__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

Hey,
I just created a pull request: https://github.com/numpy/numpy/pull/5874
Best, Stefan
On Tue, May 12, 2015 at 3:29 PM Stefan Otte stefan.otte@gmail.com wrote:
Hey,
here is an ipython notebook with benchmarks of all implementations (scroll to the bottom for plots):
https://github.com/sotte/ipynb_snippets/blob/master/2015-05%20gridspace%20-%...
Overall, Jaime's version is the fastest.
On Tue, May 12, 2015 at 2:01 PM Jaime Fernández del Río < jaime.frio@gmail.com> wrote:
On Tue, May 12, 2015 at 1:17 AM, Stefan Otte stefan.otte@gmail.com wrote:
Hello,
indeed I was looking for the cartesian product.
I timed the two stackoverflow answers and the winner is not quite as clear:
n_elements: 10 cartesian 0.00427 cartesian2 0.00172 n_elements: 100 cartesian 0.02758 cartesian2 0.01044 n_elements: 1000 cartesian 0.97628 cartesian2 1.12145 n_elements: 5000 cartesian 17.14133 cartesian2 31.12241
(This is for two arrays as parameters: np.linspace(0, 1, n_elements)) cartesian2 seems to be slower for bigger.
On my system, the following variation on Pauli's answer is 2-4x faster than his for your test cases:
def cartesian4(arrays, out=None): arrays = [np.asarray(x).ravel() for x in arrays] dtype = np.result_type(*arrays)
n = np.prod([arr.size for arr in arrays]) if out is None: out = np.empty((len(arrays), n), dtype=dtype) else: out = out.T for j, arr in enumerate(arrays): n /= arr.size out.shape = (len(arrays), -1, arr.size, n) out[j] = arr[np.newaxis, :, np.newaxis] out.shape = (len(arrays), -1) return out.T
I'd really appreciate if this was be part of numpy. Should I create a pull request?
There hasn't been any opposition, quite the contrary, so yes, I would go ahead an create that PR. I somehow feel this belongs with the set operations, rather than with the indexing ones. Other thoughts?
Also for consideration: should it work on flattened arrays? or should we give it an axis argument, and then "broadcast on the rest", a la generalized ufunc?
Jaime
-- (__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
participants (5)
-
Jaime Fernández del Río
-
Johannes Kulick
-
Nathaniel Smith
-
Stefan Otte
-
Stefan van der Walt