Ah, fantastic, thanks Per! I'd still be interested to hear from the core devs as to why this isn't the default, both with meshgrid and mgrid...
Juan. On 9 Mar 2017, 6:29 PM +1100, per.brodtk...@ffi.no, wrote: > Hi, Juan. > > Meshgrid can actually give what you want, but you must use the options: > copy=False and indexing=’ij’. > > In [7]: %timeit np.meshgrid(np.arange(512), np.arange(512)) > 1000 loops, best of 3: 1.24 ms per loop > > In [8]: %timeit np.meshgrid(np.arange(512), np.arange(512), copy=False) > 10000 loops, best of 3: 27 µs per loop > > In [9]: %timeit np.meshgrid(np.arange(512), np.arange(512), copy=False, > indexing='ij') > 10000 loops, best of 3: 23 µs per loop > > Best regards > Per A. Brodtkorb > > From: NumPy-Discussion [mailto:numpy-discussion-boun...@scipy.org] On Behalf > Of Juan Nunez-Iglesias > Sent: 9. mars 2017 04:20 > To: Discussion of Numerical Python > Subject: Re: [Numpy-discussion] Why do mgrid and meshgrid not return > broadcast arrays? > > Hi Warren, > > ogrid doesn’t solve my problem. Note that my code returns arrays that would > evaluate as equal to the mgrid output. It’s just that they are copied in > mgrid into a giant array, instead of broadcast: > > > In [176]: a0, b0 = np.mgrid[:5, :5] > > In [177]: a1, b1 = th.broadcast_mgrid((np.arange(5), np.arange(5))) > > In [178]: a0 > Out[178]: > array([[0, 0, 0, 0, 0], > [1, 1, 1, 1, 1], > [2, 2, 2, 2, 2], > [3, 3, 3, 3, 3], > [4, 4, 4, 4, 4]]) > > In [179]: a1 > Out[179]: > array([[0, 0, 0, 0, 0], > [1, 1, 1, 1, 1], > [2, 2, 2, 2, 2], > [3, 3, 3, 3, 3], > [4, 4, 4, 4, 4]]) > > In [180]: a0.strides > Out[180]: (40, 8) > > In [181]: a1.strides > Out[181]: (8, 0) > > > > On 9 Mar 2017, 2:05 PM +1100, Warren Weckesser <warren.weckes...@gmail.com>, > wrote: > > > > On Wed, Mar 8, 2017 at 9:48 PM, Juan Nunez-Iglesias <jni.s...@gmail.com> > wrote: > I was a bit surprised to discover that both meshgrid nor mgrid return fully > instantiated arrays, when simple broadcasting (ie with stride=0 for other > axes) is functionally identical and happens much, much faster. > > > Take a look at ogrid: > https://docs.scipy.org/doc/numpy/reference/generated/numpy.ogrid.html > Warren > > > I wrote my own function to do this: > > > > > > def broadcast_mgrid(arrays): > > shape = tuple(map(len, arrays)) > > ndim = len(shape) > > result = [] > > for i, arr in enumerate(arrays, start=1): > > reshaped = np.broadcast_to(arr[[...] + [np.newaxis] * (ndim - i)], > > shape) > > result.append(reshaped) > > return result > > > > > > For even a modest-sized 512 x 512 grid, this version is close to 100x > > faster: > > > > > > In [154]: %timeit th.broadcast_mgrid((np.arange(512), np.arange(512))) > > 10000 loops, best of 3: 25.9 µs per loop > > > > In [156]: %timeit np.meshgrid(np.arange(512), np.arange(512)) > > 100 loops, best of 3: 2.02 ms per loop > > > > In [157]: %timeit np.mgrid[:512, :512] > > 100 loops, best of 3: 4.84 ms per loop > > > > > > Is there a conscious design decision as to why this isn’t what > > meshgrid/mgrid do already? Or would a PR be welcome to do this? > > > > Thanks, > > > > Juan. > > > > _______________________________________________ > > NumPy-Discussion mailing list > > NumPy-Discussion@scipy.org > > https://mail.scipy.org/mailman/listinfo/numpy-discussion > > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > https://mail.scipy.org/mailman/listinfo/numpy-discussion > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > https://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion
