Re: [Numpy-discussion] Why do mgrid and meshgrid not return broadcast arrays?

2017-03-09 Thread 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:numpy-discussion-boun...@scipy.org] 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, 
per.brodtk...@ffi.no<mailto: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)
1 loops, best of 3: 27 µs per loop

In [9]: %timeit np.meshgrid(np.arange(512), np.arange(512), copy=False, 
indexing='ij')
1 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<mailto:warren.weckes...@gmail.com>>, wrote:


On Wed, Mar 8, 2017 at 9:48 PM, Juan Nunez-Iglesias 
<jni.s...@gmail.com<mailto: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)))
1 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<mailto:NumPy-Discussion@scipy.org>
https://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org<mailto:NumPy-Discussion@scipy.org>
https://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org<mailto: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

Re: [Numpy-discussion] Why do mgrid and meshgrid not return broadcast arrays?

2017-03-08 Thread 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, 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)
> 1 loops, best of 3: 27 µs per loop
>
> In [9]: %timeit np.meshgrid(np.arange(512), np.arange(512), copy=False, 
> indexing='ij')
> 1 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)))
> > 1 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

Re: [Numpy-discussion] Why do mgrid and meshgrid not return broadcast arrays?

2017-03-08 Thread 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)
1 loops, best of 3: 27 µs per loop

In [9]: %timeit np.meshgrid(np.arange(512), np.arange(512), copy=False, 
indexing='ij')
1 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<mailto:warren.weckes...@gmail.com>>, wrote:



On Wed, Mar 8, 2017 at 9:48 PM, Juan Nunez-Iglesias 
<jni.s...@gmail.com<mailto: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)))
1 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<mailto:NumPy-Discussion@scipy.org>
https://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org<mailto: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

Re: [Numpy-discussion] Why do mgrid and meshgrid not return broadcast arrays?

2017-03-08 Thread 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 , 
wrote:
>
>
> > On Wed, Mar 8, 2017 at 9:48 PM, Juan Nunez-Iglesias  
> > 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)))
> > > 1 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

Re: [Numpy-discussion] Why do mgrid and meshgrid not return broadcast arrays?

2017-03-08 Thread Warren Weckesser 
On Wed, Mar 8, 2017 at 9:48 PM, Juan Nunez-Iglesias 
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)))*
> *1 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