Why do mgrid and meshgrid not return broadcast arrays?

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Why do mgrid and meshgrid not return broadcast arrays?

Juan Nunez-Iglesias
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.

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
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Why do mgrid and meshgrid not return broadcast arrays?

Warren Weckesser-2


On Wed, Mar 8, 2017 at 9:48 PM, Juan Nunez-Iglesias <[hidden email]> 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.



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
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion



_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Why do mgrid and meshgrid not return broadcast arrays?

Juan Nunez-Iglesias
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 <[hidden email]>, wrote:


On Wed, Mar 8, 2017 at 9:48 PM, Juan Nunez-Iglesias <[hidden email]> 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.



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
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion


_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Why do mgrid and meshgrid not return broadcast arrays?

Per.Brodtkorb

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:[hidden email]] 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 <[hidden email]>, wrote:

 

 

On Wed, Mar 8, 2017 at 9:48 PM, Juan Nunez-Iglesias <[hidden email]> 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.

 

 

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
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion

 

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion


_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Why do mgrid and meshgrid not return broadcast arrays?

Juan Nunez-Iglesias
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, [hidden email], 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:[hidden email]] 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 <[hidden email]>, wrote:

 

 

On Wed, Mar 8, 2017 at 9:48 PM, Juan Nunez-Iglesias <[hidden email]> 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.

 

 

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
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion

 

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Why do mgrid and meshgrid not return broadcast arrays?

Per.Brodtkorb

The reason for returning copies from meshgrid as default instead of views into to input arrays, was to not break backwards compatibility.

The old meshgrid returned copied arrays, which is safe if you need to write to those arrays.

If you use copy=False, a view into the original arrays are returned in order to

conserve memory, but will likely return non-contiguous  arrays.  Furthermore, more than one element of a broadcast array

may refer to a single memory location.

 

Per A

 

From: NumPy-Discussion [mailto:[hidden email]] On Behalf Of Juan Nunez-Iglesias
Sent: 9. mars 2017 08:34
To: Discussion of Numerical Python
Subject: Re: [Numpy-discussion] Why do mgrid and meshgrid not return broadcast arrays?

 

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, [hidden email], 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 [[hidden email]] 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 <[hidden email]>, wrote:

 

 

On Wed, Mar 8, 2017 at 9:48 PM, Juan Nunez-Iglesias <[hidden email]> 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.

 

 

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
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion

 

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion


_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/numpy-discussion
Loading...