I am a bit worried that someone might use broadcast arrays with more then 32 arguments, then it would be a regression, I am not really worried for the new functions.

@seberg In fact, the test failure I'm getting here is for exactly that reason. So I guess I'll switch broadcast_arrays back to the more verbose implementation.

Since it is very simple and if the speed difference is big, we could think of making a fast pat. Though maybe it is a bad idea considering all the funny little things with subclasses (at least if subok=True). and the general maintanence overhead of it (and if just that we need to duplicate all tests)...

Thanks for keeping the subclasses! (after all the time @seberg and I spent on it...). One question (which I had earlier myself but didn't address): is the calculation of the shape much better/faster than using np.broadcast? If so, should that be adapted to use the routine here? If not, could np.broadcast be used here?

I wouldn't bother trying to mess around with fast paths until someone finds
a situation where broadcast_arrays is actually a bottleneck (which seems
kinda unlikely to me honestly). Let's save fastpath complexity for the
cases where it matters.

I think the idea is that we want to steer users away from np.broadcast b/c
the way it returns an iterator is confusing and unnecessary as compared to
broadcast_arrays. So better if we avoid depending on it in new code.
(Should we formally deprecate np.broadcast?)

I don't think there's much point in adapting np.broadcast to use the new
code either -- it'd be pretty complicated.

On Tue, Dec 16, 2014 at 9:08 PM, Marten van Kerkwijk <
notifications@github.com> wrote:

Thanks for keeping the subclasses! (after all the time @seberg
https://github.com/seberg and I spent on it...). One question (which I
had earlier myself but didn't address): is the calculation of the shape
much better/faster than using np.broadcast? If so, should that be adapted
to use the routine here? If not, could np.broadcast be used here?

—
Reply to this email directly or view it on GitHub
#5371 (comment).

Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh
http://vorpus.org

Hmm, just been recommending np.broadcast twice in the context of astropy (e.g., liberfa/erfa-wrapping-experiments#9 (comment)) -- to replace something very similar to broadcast_shape here, as it was quite significantly faster (and it somewhat mattered, as we were calling C code and quickly dominated by setup overhead). Indeed,

a = np.ones((3,1,3,4))
b = np.ones((3,1,4))
c = np.ones((3,1,1,4))
%timeit broadcast_shape(a,b,c)
# 100000 loops, best of 3: 9.58 µs per loop
%timeit np.broadcast(a,b,c).shape
# 1000000 loops, best of 3: 764 ns per loop

Hooray for data! If it does matter then it matters.

...looking at that thread, it sounds like the overall operation you're
talking about is 4.something ms? Does ~9 us matter?

On Tue, Dec 16, 2014 at 9:31 PM, Marten van Kerkwijk <
notifications@github.com> wrote:

Hmm, just been recommending np.broadcast twice in the context of astropy
(e.g., liberfa/erfa-wrapping-experiments#9 (comment)
liberfa/erfa-wrapping-experiments#9 (comment))
-- to replace something very similar to broadcast_shape here, as it was
quite significantly faster (and it somewhat mattered, as we were calling C
code and quickly dominated by setup overhead). Indeed,

a = np.ones((3,1,3,4))
b = np.ones((3,1,4))
c = np.ones((3,1,1,4))
%timeit broadcast_shape(a,b,c)

100000 loops, best of 3: 9.58 µs per loop

%timeit np.broadcast(a,b,c).shape

1000000 loops, best of 3: 764 ns per loop

—
Reply to this email directly or view it on GitHub
#5371 (comment).

Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh
http://vorpus.org

The above said, in most cases I think people would indeed just want the final shape, not the broadcast object returned by np.broadcast. This might suggest exposing broadcast_shape.

Also a small note: again in the context of astropy, yet another implementation of the same thing could not be replaced by np.broadcast because it needed to pass in a known shape which did not have an array yet. It might make broadcast_shape more general if tuples were interpreted as shape rather than as 1-d arrays.

@njsmith - the 4 ms was for a 10000-element array. But for the simpler functions and fewer elements it does matter. We did quite elaborate timing tests; eg., astropy/astropy#3141 (comment). For small arrays, total time is 20 µs, so, yes, the additional 10 µs does matter a little.

I think that can be interpreted as a request for my suggestion (1) in this
email?
http://thread.gmane.org/gmane.comp.python.numeric.general/59443/focus=59494

(Well, that talks about broadcast_arrays, but the generalization to
broadcast_shape is obvious.)

On Tue, Dec 16, 2014 at 9:37 PM, Marten van Kerkwijk <
notifications@github.com> wrote:

The above said, in most cases I think people would indeed just want the
final shape, not the broadcast object returned by np.broadcast. This
might suggest exposing broadcast_shape.

Also a small note: again in the context of astropy, yet another
implementation of the same thing could not be replaced by np.broadcast
because it needed to pass in a known shape which did not have an array yet.
It might make broadcast_shape more general if tuples were interpreted as
shape rather than as 1-d arrays.

—
Reply to this email directly or view it on GitHub
#5371 (comment).

Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh
http://vorpus.org

yes, looks (very) simiar...

Does that 20 µs matter though? (Genuine question, not trying to denigrate
your use case. Is this something where you perform the 20 µs operation a
million times in a loop so it adds up?)

On Tue, Dec 16, 2014 at 9:40 PM, Marten van Kerkwijk <
notifications@github.com> wrote:

@njsmith https://github.com/njsmith - the 4 ms was for a 10000-element
array. But for the simpler functions and fewer elements it does matter. We
did quite elaborate timing tests; eg., astropy/astropy#3141 (comment)
astropy/astropy#3141 (comment). For
small arrays, total time is 20 µs, so, yes, the additional 10 µs does
matter a little.

—
Reply to this email directly or view it on GitHub
#5371 (comment).

Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh
http://vorpus.org

For the purposes of discussion, I just pushed a fastpath for broadcast_shape that works for <32 dimensions (without new tests yet).

Interestingly, I needed to use the old iterator, though. This implementation with the new iterator almost works:

def broadcast_shape(*args):
    return np.nditer(args, flags=['multi_index', 'zerosize_ok']).shape

However, it fails with ValueError: Iterator is past the end if any of the arrays have size zero. That isn't a problem for np.broadcast(...).shape.

We can also add a fastpath for broadcast_arrays, too, if desired (I have original implementation in my first commit).

The other option is to keep a single path for broadcast_arrays by porting it to Cython. That would probably make up for the 10x speed factor difference. Also, it would let us treat tuples as shapes in broadcast_shape without an appreciable slow down.

@shoyer, how about just using np.broadcast(...).shape? Since it is already written in C, it should be hard to beat and avoids duplication of code.

@njsmith - the erfa routines underlie some of the coordinate and time transformations we do in astronomy, so are quite basic, and are used a lot. That said, if one really were calling them lots of times, it would probably be possible to optimize them by using larger arrays instead. Anyway, for that particular PR I mostly felt it was good to reduce many lines to one -- we try hard not to duplicate numpy functionality. But within numpy of course, it is good to think what is the best approach.

@mhvk I would use np.broadcast(...).shape only, but that imposes the 32 argument limit -- which apparently would break the existing API for broadcast_arrays.

or one could fix the 32 argument limit of the iterator, someone want to finish this:
https://github.com/juliantaylor/numpy/tree/no-max-args (turns out fixing einsum is tricky)

Well, the current code is slow for small arrays, whether this matters in any real world applications, I don't know. I have used np.broadcast(...).shape, too. Possibly even within numpy source, but I agree using the broadcasted iterator from python is probably almost never of any use at all. Or you are such an advanced user that you could even mess with nditer. On the other hand, I would like to remind everyone that broadcast_arrays has its own problems in the sense that 0-stride arrays can easily confuse new users if they are not aware of their behaviour.

OK, fixing @juliantaylor's tree would obviously be best, but in the meantime how about the following to ensure one is both fast and does not break the API?:

def broadcast_shape2(*args):
    if len(args) <= 32:
        return np.broadcast(*args).shape

    b = np.broadcast(*args[:32])
    for i in range(32, len(args), 31):
        b = np.broadcast(b, *args[i:min(i+31, len(args))])
    return b.shape

a = [np.ones((4,1)) for i in range(50)]
%timeit broadcast_shape(*a)
# 10000 loops, best of 3: 80 µs per loop
%timeit broadcast_shape2(*a)
# 100000 loops, best of 3: 11.8 µs per loop

(here, broadcast_shape is the implementation of the current PR, as in my earlier trials).

On Wed, Dec 17, 2014 at 12:39 AM, Julian Taylor notifications@github.com
wrote:

or one could fix the 32 argument limit of the iterator, someone want to
finish this:
https://github.com/juliantaylor/numpy/tree/no-max-args (turns out fixing
einsum is tricky)

I glanced over your code, and you have a few variable sized arrays here and
there, e.g in npyiter_get_common_dtype in nditer_constr.c that I am afraid
MSVC would frown upon.

In any case, if the trouble is with einsum, can this not be partially
fixed, starting with the nditer?

I still have trouble imagining a real world application where the 32
argument limit is a real limitation...

Jaime

(__/)
( O.o)
( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes
de dominación mundial.

@mhvk Good idea for simplifying broadcast_shape -- I just pushed it in the latest commit.

I like the idea of making broadcast_shape handle tuples as shapes, but even doing [isinstance(arg, tuple) for arg in args] slows down broadcast_shape by a factor of 70% or so.

@shoyer - had forgotten how forgiving python is with indices outside of a list's range... nice. Seems OK to leave the tuple is shape thing for another PR (would also need to be careful that it wouldn't break the API).

@mhvk In that case, we should not expose broadcast_shape in the public API, because it currently casts the input to ndarrays with np.asarray, so (2, 3) gets converted into an ndarray with shape (2,).

@shoyer - agreed that for now it is best not to expose it, as it currently is mostly a work-around the 32-argument limit of np.broadcast. But very useful here!

Instead of handling tuples as shapes, let's put them in a separate kwarg.
isinstance is almost always a code smell.
On 17 Dec 2014 15:45, "Marten van Kerkwijk" notifications@github.com
wrote:

@shoyer https://github.com/shoyer - agreed that for now it is best not
to expose it, as it currently is mostly a work-around the 32-argument limit
of np.broadcast. But very useful here!

—
Reply to this email directly or view it on GitHub
#5371 (comment).

Anyone have opinions on my question from the first post: should we allow broadcasting with negative numbers in the shape?

Given that nditer supports it and that it could be useful, my inclination is to document it and add tests to verify that it works as expected.

The fact that after 30 seconds thought I can't figure out what negative
numbers would mean in this situation makes me inclined not to support it.

For reshape you know what the total array size is, so an ndim array has
ndim-1 degrees of freedom in its shape. Here the final array size could be
anything, so... I don't get it.

I'm generally against supporting things just because we can. That way lies
confusion and difficult maintenance and annoying deprecation cycles.
On 17 Dec 2014 21:36, "Stephan Hoyer" notifications@github.com wrote:

Anyone have opinions on my question from the first post: should we allow
broadcasting with negative numbers in the shape?

Given that nditer supports it and that it could be useful, my inclination
is to document it and add tests to verify that it works as expected.

—
Reply to this email directly or view it on GitHub
#5371 (comment).

@njsmith In the context of nditer, apparently negative numbers in the shape mean: size of the original dimension (if it's present) or 1 otherwise. So, for example, np.broadcast_to(x, (-1, -1, -1)) ends up being roughly equivalent to np.atleast_3d, except it follows broadcasting rules, so the extra dimensions always end up inserted at the front.

We are >< that close to finishing this...

The readonly flag can be set separately if needed?

Yeah, sorry for delaying things in the first place. It can be added later, but if we say it makes sense, I would think that it probably makes sense to make it default to True (for the new function only unfortunatly)? And in that case we need to decide before the release, but we can just make a blocker issue to not forget it.

How about we just default it to readonly now (which does seem safer) and
then if we decide to change our minds later, then oh well it's easy because
going readonly -> readwrite is backwards compatible.
On 14 Feb 2015 22:05, "seberg" notifications@github.com wrote:

Yeah, sorry for delaying things in the first place. It can be added later,
but if we say it makes sense, I would think that it probably makes sense to
make it default to True (for the new function only unfortunatly)? And in
that case we need to decide before the release, but we can just make a
blocker issue to not forget it.

—
Reply to this email directly or view it on GitHub
#5371 (comment).

Readonly default sounds reasonable to me.

Should we even include a readonly keyword argument? readonly=False seems so niche that we might as well insist users set the writeable flag directly.

The only justification I can see for readonly is as a temporary argument to ease the transition to readonly-only at some point in the future. For example, we could default to readonly=False in broadcast_arrays, and then issue a warning if it hasn't been set to True.

Maybe that makes sense, I am still wondering if there is some nice way to deprecate broadcast_arrays being writeable, but I can't think of one easily, or can we have a "writable but warn" flag without much hassle?

We already implemented writeable-but-warn for diagonal and for some record
array thing. So it's a trivial matter of just using it, or maybe
resurrecting it from recent git history.

And leaving out the readonly flag altogether for the new api sounds fine to
me.
On Feb 20, 2015 5:40 AM, "seberg" notifications@github.com wrote:

Maybe that makes sense, I am still wondering if there is some nice way to
deprecate broadcast_arrays being writeable, but I can't think of one
easily, or can we have a "writable but warn" flag without much hassle?

—
Reply to this email directly or view it on GitHub
#5371 (comment).

Ok, lets move this finally. Since I guess it is a bit beyond this (though if you are happy to help with it, please do!) to do a deprecation for broadcast_arrays, we don't need to do it in this PR. But there seems to be an agreement that there is no reason to have a writable output for now.

Could you make the output readonly (and hack broadcast_arrays to not change for now)? That would be great! I think it might be easiest to make a private function that has a readonly keyword arg for this purpose, so that there is an argument, but maybe the user doesn't need to know about it.

Then we can finally put this in. Again, I am sorry for delaying it all the time....

@seberg Sounds good, I'll get on that.

OK, I have pushed another commit that makes the result of broadcast_to readonly, while keeping the results of broadcast_arrays writeable.

I have not added the warn-on-write flag for broadcast_arrays for two reasons:

1. The recent debate about breaking backwards compatibility for np.diagonal is still unresolved. This issue feels somewhat similar.
2. As far as I can tell from the original PR (Clean up diagonal #280), the warn-on-write flag can currently only be set from the C-side. I don't really know C, so I'm not in a good position to add the appropriate hack to let it be set from the Python side.

OK, great. I will look at it later (hopefully), if I don't, wake me up some time ;). I think we are much safer here then with diagonal to be honest, though it might be a misconception because I simply can't see anyone writing to a non-contiguous multi-dimensional array without risking insane bugs ;).

Sounds like this is ready. I'll merge this evening if no one complains or beats me to it.

@shoyer Could you squash the commits into one?

Thanks @shoyer.

Thanks a lot!

Leaving this here in case someone needs to support older versions of numpy. I timed it against the version of broadcast_to @shoyer included in xray. It's much less efficient, of course, and involves a memory allocation, but for many applications the difference is negligible.

def naive_broadcast_to(a, shape):
    "Broadcast an array to a specified shape. Returns a view."
    return np.broadcast_arrays(a, np.empty(shape, dtype=np.bool))[0]

a = np.random.rand(100, 1, 100)
np.testing.assert_array_equal(naive_broadcast_to(a, (100, 1000, 100)), broadcast_to(a, (100, 1000, 100)))
%timeit naive_broadcast_to(a, (100, 10000, 100))
%timeit broadcast_to(a, (100, 10000, 100))

10000 loops, best of 3: 56.7 µs per loop
The slowest run took 4.38 times longer than the fastest. This could mean that an intermediate result is being cached 
100000 loops, best of 3: 7.69 µs per loop