<div dir="ltr"><div>Hi Nathaniel,<br><br></div>Thanks for the detailed thoughts.<br><div><div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Jul 3, 2018 at 4:27 AM, Nathaniel Smith <span dir="ltr"><<a href="mailto:njs@pobox.com" target="_blank">njs@pobox.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><span class="gmail-">On Sat, Jun 30, 2018 at 6:51 AM, Marten van Kerkwijk<br>
<<a href="mailto:m.h.vankerkwijk@gmail.com">m.h.vankerkwijk@gmail.com</a>> wrote:<br>
> Hi All,<br>
><br>
> In case it was missed because people have tuned out of the thread: Matti and<br>
> I proposed last Tuesday to accept NEP 20 (on coming Tuesday, as per NEP 0),<br>
> which introduces notation for generalized ufuncs allowing fixed, flexible<br>
> and broadcastable core dimensions. For one thing, this will allow Matti to<br>
> finish his work on making matmul a gufunc.<br>
><br>
> See <a href="http://www.numpy.org/neps/nep-0020-gufunc-signature-enhancement.html" rel="noreferrer" target="_blank">http://www.numpy.org/neps/nep-<wbr>0020-gufunc-signature-<wbr>enhancement.html</a><br>
<br>
</span>So I still have some of the same concerns as before...<br>
<br>
For the possibly missing dimensions: matmul is really important, and<br>
making it a gufunc solves the problem of making it overridable by duck<br>
arrays (via __array_ufunc__). Also, it will help later when we rework<br>
dtypes: new dtypes will be able to implement matmul by the normal<br>
ufunc loop registration mechanism, which is much nicer than the<br>
current system where every dtype has a special-case method just for<br>
handling matmul.</blockquote><div><br></div><div>Indeed, I only became recently aware of the ->dot member of the dtype struct. Pretty ugly!<br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"> The ? proposal isn't the most elegant idea ever, but<br>
we've been tossing around ideas for solving these problems for a<br>
while, and so far this seems to be the least-bad one, so... sure,<br>
let's do it. <br></blockquote><div><br></div><div>Yes, I think this ugliness is the price we pay for matmul not just meaning matrix multiplication but also no-broadcast vector-matrix, matrix-vector, and vector-vector. (Of course, my after-the-fact annoyance with that design does make me more sensitive to the argument that one should not try to shoehorn gufuncs into cases they are not meant for.) <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<br>
For the fixed-size dimensions: this makes me nervous. It's aimed at a<br>
real use case, which is a major point in it's favor. But a few things<br>
make me wary. For input dimensions, it's sugar – the gufunc loop can<br>
already raise an error if it doesn't like the size it gets.</blockquote><div><br></div><div>Can it? I would never have thought that the inner loop would need to do any checking; it is certainly not obvious from the code or from the documentation. Does the iterator itself check for errors every iteration? If so, that might be one place for a quick speed-up for ufunc.at...<br> <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"> For output<br>
dimensions, it does solve a real problem. But... only part of it. It's<br>
awkward that right now you only have a few limited ways to choose<br>
output dimensions, but this just extends the list of special cases,<br>
rather than solving the underlying problem. For example,<br>
'np.linalg.qr' needs a much more generic mechanism to choose output<br>
shape, and parametrized dtypes will need a much more generic mechanism<br>
to choose output dtype, so we're definitely going to end up with some<br>
phase where arbitrary code gets to describe the output array.</blockquote><div><br>I think this is a much rarer case, which will indeed need some type of hook no matter what.<br>(I have been thinking about solving that by moving the call to the type resolver earlier. That one gets passed `op`, so it can create an output array if it doesn't exist; it doesn't need any of the sizes.)<br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"> Are we<br>
going to look back on fixed-size dimensions as a quirky, redundant<br>
thing?<br></blockquote><div><br></div><div></div><div>For this particular case, I find the signature so much clearer that that in itself is a reason to do it (readability counts and all that).<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<br>
Also, as currently proposed, it seems to rule out the possibility of<br>
using name-based axis specification in the future, right? (See<br>
<a href="https://github.com/numpy/numpy/pull/8819#issuecomment-366329325" rel="noreferrer" target="_blank">https://github.com/numpy/<wbr>numpy/pull/8819#issuecomment-<wbr>366329325</a>) Are<br>
we sure we want to do that?<br></blockquote><div><br></div><div>I think it only excludes having the choice of keying the dict with the default axis number *or* with its dimension name, which I think would in fact be a good idea: if one uses a dict to describe axes entries, the keys should be the names of the axes, where name can be one of the fixed numbers.<br></div><div>(Actually, strictly, since the default axis number is always negative, one can even have both. But that would be awful.)<br></div><div><br></div><div>Should add that I don't think the scheme will work all that well anyway - there are quite a few cases where one duplicates a name (say, a square matrix (n,n)), for which keying by "n" would be less than useful.<br></div><div> <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
If everyone else is comfortable with all these things then I won't<br>
block it though.<br>
<br>
For broadcasting: I'm sorry, but I think I'm -1 on this. I feel like<br>
it falls into a classic anti-pattern in numpy, where someone sees a<br>
cool thing they could do and then goes looking for problems to justify<br>
it. (A red flag for me is that "it's easy to implement" keeps being<br>
mentioned as justification for doing it.) The all_equal and<br>
weighted_mean examples both feel pretty artificial -- traditionally<br>
we've always implemented these kinds of functions as regular functions<br>
that use (g)ufuncs internally, and it's worked fine (cf. np.allclose,<br>
ndarray.mean). In fact in some sense the whole point of numpy is to<br>
help people implement functions like this, without having to write<br>
their own gufuncs. Is there some reason these need to be gufuncs? And<br>
if there is, are these the only things that need to be gufuncs, or is<br>
there a broader class we're missing? The design just doesn't feel<br>
well-justified to me.<br></blockquote><div><br></div><div>Those are all fair points. For all_equal one cannot really write a separate function easily given the absence of short-circuiting. But possibly that just argues for the need to short-circuit...<br> <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
And in the past, when we've implemented things like this, where the<br>
use cases are thin but hey why not it's easy to do, it's ended up<br>
causing two problems: first people start trying to force it into cases<br>
where it doesn't quite work, which makes everyone unhappy... and then<br>
when we eventually do try to solve the problem properly, we end up<br>
having to do elaborate workarounds to keep the old not-quite-working<br>
use cases from breaking.<br></blockquote><div><br></div><div>Given my sentiments about the multiple meanings of `@`, I'm sensitive to this argument.<br><br>But as Matti pointed out, we do not have the accept the whole NEP. Indeed, the broadcasting is in a separate PR.<br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<br>
I'm pretty sure we're going to end up rewriting most of the ufunc code<br>
over the next few years as we ramp up duck array and user dtype<br>
support, and it's already going to be very difficult, both to design<br>
in the first place and then to implement while carefully keeping shims<br>
to keep all the old stuff working. Adding features has a very real<br>
cost, because it adds extra constraints that all this future work will<br>
have to work around. I don't think this meets the bar.<br></blockquote><div><br></div><div>I think by now it is clear that moving incrementally is the way to go; the ufunc code is in fact being rewritten, if slowly.<br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
By the way, I also think we're getting well past the point where we<br>
should be switching from a string-based DSL to a more structured<br>
representation. (This is another trap that numpy tends to fall into...<br>
the dtype "language" is also a major offender.) This isn't really a<br>
commentary on any part of this in particular, but just something that<br>
I've been noticing and wanted to mention :-).<span class="gmail-im gmail-HOEnZb"><br></span></blockquote><div><br></div><div>Yet, dtype('(3,4)f8') is really clear, unlike all the other forms... It is similar to the string formatting mini-language. Anyway, that is a bit off topic.<br><br></div><div>Overall, would one way to move forward be to merge the first PR (flexible and frozen) and defer the broadcastable dimensions?<br><br></div><div>All the best,<br><br></div><div>Marten<br></div><div><br></div><div>p.s.  I'm amused that the broadcastable dimensions were in fact the default originally. At some point, I should try to find out why that default was changed.<br></div></div></div></div></div></div>