WIP: gufunc core dimensions should not broadcast

[jaimefrio]

No description provided.

[jaimefrio]

Some working code to substantiate the discussion. This currently disallows non-matching core dimensions, but prepends dimensions of size 1 as needed to complete the core dimensions. Some examples of the behavior...

First the obvious ones:

>>> inner1d([1,1,1], [1])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: inner1d: Operand 1 has a mismatch in its core dimension 0, with gufunc signature (i),(i)->() (size 1 is different from 3)
>>> matrix_multiply([[1], [1]], [[1], [1]])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: matrix_multiply: Operand 1 has a mismatch in its core dimension 0, with gufunc signature (m,n),(n,p)->(m,p) (size 2 is different from 1)

The prepending dimensions of size 1 thing at work:

>>> matrix_multiply([1,1], [[1], [1]])
array([[2]])
>>> matrix_multiply([[1], [1]], [1, 1, 1])
array([[1, 1, 1],
       [1, 1, 1]])

And a failure that needs a rewording of the error, because the non-matching dimension is a prepended 1, so the error is misleading:

>>> matrix_multiply([[1, 1], [1, 1]], [1, 1, 1])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: matrix_multiply: Operand 1 has a mismatch in its core dimension 0, with gufunc signature (m,n),(n,p)->(m,p) (size 1 is different from 2)

[charris]

I think that the core dimensions should not be filled out, i.e., if the signature is '(L,M), (M,N) -> (L,N)', then then each argument should have at least dimension 2. A trickier case is probably '(L,M), (M) -> (L)', where the arguments should be of at least dimension 2 and1 respectively. Does the broadcasting works as expected in this case. Those signatures won't work for scalars, but I think that is the right thing to do.

[njsmith]

-1 on the dimension prepending until the mailing list discussion is resolved. (I really don't understand why anyone thinks this is a good idea? Why would we want the equivalent of dot([[1, 2, 3]], [1, 2, 3]) to succeed, when it's obviously a transposition error and will give the wrong answer?)

[njsmith]

Oh, now it sounds like Chuck and I are in agreement after all :-)

[njsmith]

btw you have universal build failure

[jaimefrio]

Yes, the failing tests are the ones checking that gufuncs do what this PR no longer allows...

The don't-prepend-ones behavior is already in place for output args, so getting it for inputs is trivial.

[jaimefrio]

So now an exact match between the sizes of core dimensions is required, and all core dimensions must be present. I kind of like the result, but It feels like too big of a change to do so happily...

[charris]

If this goes in, and I'm in favor of that pending the discussion on the list, it needs to be documented in the release notes, and also in doc/source/reference/c-api.generalized-ufuncs.rst. The latter is inadequate even for the old behavior, as it doesn't completely document the new functions. It should look more like doc/source/reference/c-api.ufunc.rst.

[jaimefrio]

This is probably also a good opportunity to put together a more comprehensive set of tests that validate that the behavior is what we have been discussing.

[charris]

I have some functions that we can use for that if we want to eventually remove umath_tests.c.src.

[jaimefrio]

@charris Take a look at the new code when you have a chance, and let me know what you think. I have tried to minimize code repetition, so I may have departed a little bit from what we discussed on the list.

I haven't tested the new path (creating a ufunc with PyUFunc_FromFuncAndDataAndStrictSignature will set core_enabled to 2 instead of 1, and trigger the more strict rules when the gufunc gets called), but the old one seems to be working fine. Will start looking into testing this tomorrow.

[pv]

Probably best to replace the magic numbers "0", "1" and "2" by named
constants "CORE_DISABLED", "CORE_RELAXED", "CORE_STRICT". Probably also
best to leave room for considering it a bitmask in the future...

[njsmith]

So if I understand the plan correctly:

• the "broadcasting" here has zero known use cases
• we plan to deprecate the broadcasting version and remove it
• therefore we'll be leaving all existing ufuncs in broadcasting mode for
  the deprecation cycle

So, uh... Who is this strict mode for? It seems like no one will ever use
it? What am I missing?

I'm also curious why we aren't just treating this as a bug - the behaviour
seems pretty useless, inconsistent with how all "gufunc-like" operations
have always worked in numpy (e.g. dot), and I have no knowledge of any
gufuncs being used in the wild before 1.8. I mean, if everyone just wants
to cautious then I'm happy with that, better to err on the side of caution
and all. But if there's more of a motivation than that, then I'm curious
what it is (plus documenting it in the tracker seems useful)? Is there some
use case that we're worried about in particular? Did someone discover a
case where someone was actually relying on broadcasting core dimensions?

[seberg]

I tend to agree, we probably don't need the current mode (or at least for now we can assume we don't), and can opt instead for a normal deprecation.

[charris]

@njsmith The current function is in the API, so the choices are

1. It's a bug, just change the behavior and document it.
2. Add a new function. Not sure why the old needs to be deprecated in that case.

I would normally pick 2 out of caution, but Jaime is suggesting that there may be more changes/options in the future,, so it seems reasonable to have one constructor for all the potential options.

[njsmith]

As a general rule, I think we should try hard to minimize the number of
options we support. As a somewhat silly example, when we added .at() we
didn't also add a new ufunc registration function or flag or something that
had to be used to enable it - we just turned it on unconditionally for all
ufuncs. And ufuncs have a lot of complexity already. So I'm -1 on
officially continuing to support core dimension broadcasting without some
compelling use case.

I don't really understand the future proofing argument - the new API
function in this patch is no more flexible than the old one, and increases
the number of functions that we'll have to handle of when making future
changes...?

@njsmith https://github.com/njsmith The current function is in the API,
so the choices are

1. It's a bug, just change the behavior and document it.
2. Add a new function. Not sure why the old needs to be deprecated in
   that case.

I would normally pick 2 out of caution, but Jaime is suggesting that there
may be more changes/options in the future,, so it seems reasonable to have
one constructor for all the potential options.

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

[charris]

About the current behavior, I think the arrays must have at least as many dimensions as in the signature, but broadcasting within that smaller scope could be useful in some circumstances, so another option is to change the original to require the right sizes but continue the broadcasting, then add a new function that doesn't allow broadcasting.

[njsmith]

Can you give a concrete example of these "some circumstances"?
On 19 Sep 2014 09:31, "Charles Harris" notifications@github.com wrote:

About the current behavior, I think the arrays must have at least as many
dimensions as in the signature, but broadcasting within that smaller scope
could be useful in some circumstances, so another option is to change the
original to require the right sizes but continue the broadcasting, then add
a new function that doesn't allow broadcasting.

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

[charris]

@njsmith It's already parameterized internally, although that is an implementation detail, the question is whether of not it would be useful to expose that parameterization. That said, I'm sympathetic to the idea of just adding new functions to the API, designing for the (unknown) future sometimes creates a mess.

One example I can think of would be multiplication of a stack of matrices by a stack of diagonal matrices. An efficient way to do that, instead of the full matrix multiplication, would be to store the diagonals as (1,n) matrices and broadcast the * operator rather than do a dot.

[charris]

Another option is to extend the signature notation, allowing an explicit 1 or . for a dimension intended to be broadcast.

[charris]

For example, (M, N), (1, N) -> (M,N) would do what is currently done, adding a new dimension of 1 and broadcasting. A . might be better than a 1, as in some cases (stacked quaternions) it might be preferable to use numbers to specify a required size, 4 in the quaternion case. But that is something for the future...

[njsmith]

For the broadcasting matrix-times-diagonal matrix case, isn't that the way
that * already works? And the matmul gufunc doesn't work this way
regardless of whether core dimension broadcasting is enabled?

You seem very gung ho about figuring out how to expose this but I feel
like we skipped the question of whether we should expose this.
On 19 Sep 2014 09:59, "Charles Harris" notifications@github.com wrote:

Another option is to extend the signature notation, allowing an explicit 1
or . for a dimension intended to be broadcast.

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

[charris]

It was a suggestion on the list, there was a "maybe" attached. I was hoping for discussion, which we now have ;) The current * will do that if the dimension 1 is explicitly added to the stacked diagonal vectors.

I'm rather liking the idea of letting the signature specify that behavior.

[njsmith]

Maybe I am just being dense but I still have no idea how you are imagining
this working :-(. AFAICT enabling this functionality for the matmat gufunc
will not do anything useful - certainly it won't do anything like your
matrix @ diagonal matrix multiplication example. So that doesn't seem like
an argument for this functionality. But * (np.multiply) is not a gufunc, so
it can't be affected by this functionality either. What am I missing?
On 19 Sep 2014 10:18, "Charles Harris" notifications@github.com wrote:

It was a suggestion on the list, there was a "maybe" attached. I was
hoping for discussion, which we now have ;) The current * will do that if
the dimension 1 is explicitly added to the stacked diagonal vectors.

I'm rather liking the idea of letting the signature specify that behavior.

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

[charris]

Matmul is a separate issue, this is just about gufuncs in general.

[njsmith]

Yes, but I still haven't seen any example I understand of why any gufunc
would want this. You seemed to have an example - which might involve
matmul, I'm not sure - but I don't understand it. So I was hoping I'd I
elaborated on my confusion then you might elaborate on your example :-)
On 19 Sep 2014 10:29, "Charles Harris" notifications@github.com wrote:

Matmul is a separate issue, this is just about gufuncs in general.

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

[jaimefrio]

I think my reference for what I would not like to happen is @seberg's rewrite of indexing: everything looked great to everyone, all the new behaviors were sensible and better than the ones they were replacing... Then we put it out and seemed to break everyone's code, and had to go back and make some of the old behavior fit in again in a hurry.

If we extend the signature parsing code following Nathaniel's idea of "a trailing ? means optional dimension", we could also go with "a trailing means core dimension broadcasts." This will add the most complexity and may still create interesting backwards compatibility issues if anyone is relying on broadcast of core dimensions in their code, so I am -0.5 on it, and would rather go the "this is a bug" route.

If we add a new constructor to the API, I would also change all the existing linalg gufuncs to the new behavior, add deprecation warnings to the old constructor, and eventually, if no major issue is found, make both constructors an alias of each other and remove the old behavior.

But quite frankly, I don't care if the bike shed gets painted green or orange...

[charris]

The thought is that you can stack vectors with broadcasting explicitly as v[..., None, :] or have the signature tell the gufunc machinery to do so. The advantage is the signature would also distinguish "broadcast that dimension" from "it's an plain old dimension, don't broadcast it". I'm not thinking of this for matmul, the matmul functions work fine with Jaime's original. The disadvantage would be making a function for that special case when it is already pretty simple with the existing machinery.

The simplest options at this point, ignoring all these possibilites, are

1. Treat current behavior as a bug and change the current function.
2. Add a new function.

The second gives backwards compatibility in case there is a function in the wild depending on the old behavior, the first is a bit of a gamble. Probably safe, but I've been surprised before...

[charris]

So +1 for Jaime's approach. We can keep the other possibilities in mind for the future.

[charris]

My answer to this question is: no. I think the long term costs of supporting a behavior like this - which includes testing, maintenance, extra complexity in making future changes, documentation, cognitive load for users who have to figure out what broadcasting means in this case, etc. - will be substantially higher than even a very careful and extended deprecation cycle.

Here is where we differ. If I thought there was a great burden in supporting the current behavior, I would agree with you. But I don't see that. For the old behavior, we can leave support and testing as is and wait for complaints, if any. For the new behavior, tests need to be written in any case, so that part is a wash. Documentation can mention up front the distinction, and that should not be difficult to understand. And the code itself is not much more complicated; code for the new behavior is available, the old code is still there, and the choice between the two is determined by the value of core_enabled. So I don't see support as a burden and by keeping the old behavior as is we maintain backward compatibility and don't need to worry if some poor soul depends on it.

[charris]

I've searched github and Google Code for any independent uses of PyUFunc_FromFuncAndDataAndSignature and can't find any, so we might as well gamble on changing the behavior of the function if that gets it in.

[jaimefrio]

Are you thinking of going cold turkey? Or will we still want to support the old behavior with a deprecation warning for the next release cycle?

[charris]

I figure it's cold turkey and keep the current name. If we add a new function we might as well support both, the upkeep isn't much.

[njsmith]

Here's an example of the kind of "costs" I'm thinking of: suppose scipy or whoever decides to add some gufuncs (maybe to scipy.linalg). Obviously they want the strict behaviour where available, but when running against an older version of numpy they would rather have the broadcasty behaviour than nothing at all. If we have two functions, then how do they accomplish this? The "easy" approach would be to write a compile-time configure check for the new function, and use it only if available. But, this won't work -- it means that if you build against an older numpy and then run against a newer numpy, you will get the old broken behaviour even though the installed numpy supports the correct strict behaviour. But if you build against a newer numpy to get the correct behaviour, then this build can no longer run against an older numpy. Instead we need to somehow check for the API version at runtime, and somehow call a C function when it exists at runtime even if it didn't exist at compile time... it's almost certainly doable somehow, I can sorta see some possible ways to do it, but my point is that any time you or anyone else spends inventing/implementing/debugging a system like this would be a massive waste of precious developer time.

Numba apparently has an API for defining user-defined gufuncs. If both options are officially supported documented things, then presumably someone will have to waste some time adding a Numba-level wrapper for selecting which behaviour they want, documenting this, etc. Again this is wasted time, given that AFAWK everyone wants the strict behaviour.

[charris]

@njsmith So are you happy with cold turkey?

[njsmith]

I'm happy with cold turkey.

On Wed, Oct 8, 2014 at 3:10 AM, Charles Harris notifications@github.com
wrote:

@njsmith https://github.com/njsmith So are you happy with cold turkey?

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

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

[jaimefrio]

A gufunc user wannabe has come out of the closet today in StackOverflow:

http://stackoverflow.com/questions/26285595/generalized-universal-function-in-numpy

[njsmith]

I don't think the request on stackoverflow is actually coherent with how gufuncs work? I left a comment anyway.

@jaimefrio: It sounds like we're at least closer to consensus... are you up for modifying this PR to make it all-strict-all-the-time?

[jaimefrio]

Yes, I should find some time in the next couple of days: we moved last week, so things are still a little hectic around here. I think I will still see if a deprecation warning can be added, if only to remove it afterwards.

[jaimefrio]

OK, redid the whole thing: cold turkey it is. Much simpler indeed. Lets see what Travis thinks...

[njsmith]

We should probably also have a testcase that would previously invoked the prepend-1 behaviour. I guess inner1d(scalar, scalar) would work.

[jaimefrio]

That's what the next test (Extend core dimensions) does.

[njsmith]

otherwise LGTM

On Thu, Oct 16, 2014 at 5:46 PM, Jaime notifications@github.com wrote:

OK, redid the whole thing: cold turkey it is. Much simpler indeed. Lets
see what Travis thinks...

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

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

[charris]

prepended.

[charris]

And same for other uses.

[jaimefrio]

I have added your suggestions, plus fixed a couple more places in the docs, and added a note to the release notes.

[charris]

@jaimefrio Just a spelling nitpick. I think this is almost ready to go.

[jaimefrio]

@charris I think I fixed them all... There is a comment of yours about pre-pended vs prepended, but that's on a removed line.

[charris]

OK, merged. Thanks Jaime.