As for performance, it's hard to get consistent results, but this adds about 5μs to 10μs, or 4% to 5% decrease in speed (tested with %timeit and multiply and add.)

Note that we want similar logic also for dot --- there are two versions of that, one in dotblas.c and the other one is array_matrixproduct in multiarraymodule.c. The argument parsing function would be best to put into a separate header file.

In the ufunc code, you might want to run the override check/parsing code only if get_ufunc_arguments sees an object with an __array_priority__. Although 5% increase is not so bad.

The question whether the dict should be keyed by ufunc objects rather than ufunc names should also be addressed. I think one cannot assume ufunc names are unique.

Considering this, should I change the name to something other than __ufunc_override__?

Also, do tests for this go in numpy/numpy/core/tests/test_umath.py? Where would I add docs?

It think __ufunc_override__ is fine even if it's also dot. Yes, tests go there. Docs go into the same place as for __array_priority__, ie., arrays.classes.rst.

I'd perhaps not expose that function in the public API, as it's not really something a user of Numpy C-API wants to see. It would be better put in a private/something.h and include it where necessary.

I'm having a problem with checking for the numpy.dot key in __ufunc_override__ since there is no numpy.dot type in the C-API for me to compare with numpy.dot. Ideas? I'm trying to check it with function pointers, but It's going kind of slow.

@cowlicks: You can either grab the Python function object via PyImport_ImportModule + PyObject_GetAttrString and store the result in a static variable inside the functions that need it (i.e., _dotblas.c:dotblas_matrixproduct and multiarraymodule.c:array_matrixproduct), or, you can grab the function object to a static variable in the module's init method.

This is really cool. I think the API's still not 100% there, though. Some thoughts:

• We should definitely be referring to ufuncs by object, not by name.
• We should be able to override ufunc methods like reduce and accumulate, not just __call__, but this doesn't seem to be supported by the proposed machinery?
• I'd also worry a bit about the (somewhat overcomplicated) ufunc calling conventions -- e.g. an out argument might be positional, or might be a kwarg, depending on the call. IMHO we should have some system for normalizing these before delegating to an override function. Anything else will create weird and subtle bugs as we force everyone to try and reimplement these details.
• The dict lookup is needlessly limiting. np.ma could support arbitrary (non-generalized) ufuncs without caring at all about overriding individual ones -- pull out the masks from the arguments, combine them to create a mask for the operation, then call the ufunc on the underlying
• I don't like the use of __array_priority__. This mechanism has always seemed to be based on a kind of wishful thinking... it would be nice if we could just look at each possible array-like object in the world and sort them into some linear order, in which every "higher" object knows how to appropriately handle every "lower" object, but it's just not true. And it's not necessary, either. I think a good goal would be for this new mechansim to be good enough to handle everything the current __array_priority__/__array_prepare__/__array_wrap__ stuff does, so that we can eventually deprecate that cumbersome and limited API.

Python's standard approach to multimethods seems like a good fit here -- I'm thinking of things like __add__, which are in fact multimethods: to perform a + b, Python first asks the a object whether it can perform the operation, and if it can't (a.__add__(b) returns NotImplemented) then it asks the b object whether it can perform it (by calling b.__radd__(a)), and if both return NotImplemented then the operation fails with a TypeError. (This is more limited than classic multimethod implementations, but IMHO superior -- basically it imposes the restriction that one of the types involved has to know how to handle the whole operation, as opposed to allowing arbitrary third party libraries to manipulate how arbitrary other types get handled.)

Of course, ufuncs can have arbitrarily many arguments, so we don't want to mess about with "rufunc" or whatever. Here's my suggestion:

• First, we normalize the ufunc arguments into a tuple of input arrays, and a dict of kwargs. So if someone calls np.add(a, b, out=c) then we have input arrays (a, b) and kwargs {"out": c}. And if someone calls np.add(a, b, c), then we again have input arrays (a, b) and kwargs {"out": c}.
• We scan the input arrays from left to right. If inputs[i] has a __numpy_ufunc__ attribute, then we call:

inputs[i].__numpy_ufunc__(ufunc_object, ufunc_method, i, inputs, kwargs)

For example, in our example above, if b had a __numpy_ufunc__ method, it would be called as:

b.__numpy_ufunc__(np.add, np.add.__call__, 1, (a, b), {"out": c})

If this returns NotImplemented, then we continue scanning the input arguments. If it returns any other value, we return that. If it raises an error, we propagate it.

• If we finish scanning the input arrays, then there are two possiblities. If we found at least one __numpy_ufunc__ attribute, then the fact that we've reached the end means that they've all returned NotImplemented. In this case, we raise TypeError. If we found no __numpy_ufunc__ attributes, then we fall back on the current ufunc dispatch behaviour.

Note 1: This scheme doesn't allow dispatch on out= or where= arguments. (And in the case of reduceat, it doesn't dispatch on the indices argument.) I think this is the right behaviour, though it can be extended without breaking compatibility in the future if I'm wrong. These other arguments all have to be ndarray-compatible objects for their semantics to make sense -- e.g., where has to be a homogenous boolean array. __array_interface__ and friends are already flexible enough to support any object which is supportable here. And besides dispatching on output type is just weird and unPythonic; where you store the result of a function shouldn't be allowed to affect what the function does.

Note 2: If a class defines __numpy_ufunc__, then it opts out of __array_wrap__ etc. entirely; they will never be called. This makes it easier to eventually transition, because no-one will be using both system together.

Note 3: I think this system is strictly more powerful than the existing __array_wrap__ stuff. To reimplement __array_priority__ and __array_wrap__, just do:

class NDArrayWithSpice(object):
  def __numpy_ufunc__(self, ufunc_object, ufunc_method, i, inputs, kwargs):
      # If someone else has higher priority, defer to them.
      this_priority = self.__array_priority__
      for arr in inputs[i + 1:]:
          if getattr(arr, "__array_priority__", this_priority) > this_priority:
              return NotImplemented
      # Delegate to the numpy implementation
      result = ufunc_method(*[np.asarray(input) for input in inputs], **kwargs)
      # Wrap the result and return it.
      return self.__array_wrap__(result)

(This leaves out __array_prepare__, because it's mostly undocumented so I'm not sure how it works. And I'm not sure anyone uses it. But it should be doable too.)

Oh right, I was going to also say:

Note 4: For a generic implementation (like the one np.ma might want), where you might need to know whether the current operation was e.g. a reduce, even if you don't care about exactly which ufunc is being reduced, you can check ufunc_method.__name__, which will be "__call__", "reduce", "accumulate", "outer", or "reduceat".

An alternative would be to have a different override method for each ufunc method. __numpy_ufunc_call__, __numpy_ufunc_reduce__, etc. I don't have much of an opinion either way on this detail. Maybe the multiple-method version is slightly cleaner? And they could have slightly different and more appropriate signatures?

Regarding np.dot: I think the solution is to turn np.dot into a gufunc? Like most of linalg is now?

@njsmith: Doesn't the weird way dot chooses the operation axes cause problems for gufunc core loop selection? I may be wrong, though. OTOH, just allowing dot use the same mechanism as ufuncs is not a big sin in my opinion.

Regarding the points:

• argument normalization: agreed
• NotImplemented instead of __array_priority__ --- agreed, I like the idea of using the same mechanism as in Python's operator overrides. Maybe checking also the override method in the output argument should be done (which would also cover __array_prepare__)?
• Exclusive __numpy_ufunc__ over __array_wrap__ / __array_prepare__: agreed --- that mechanism frankly seems overly complicated.
• I would prefer having separate reduce etc. methods. Their implementation would for sparse matrices likely be quite different from __call__, so it would be simplest to make them separate.

I thought 'dot' used the same rule as gufuncs -- iterate over the first n
axes, pass the last k to the inner operation loop? (This is sort of
annoying, actually, because 'dot' and 'ufunc.reduce' are both gufuncs in
principle, but ufunc.reduce defaults to operating on axis=0 and dot on
axis=(-2, -1). I assumed that the reason gufuncs work the way they do was
for compatibility with dot.) Or are you thinking of something else?

Do you have any use case for dispatching on the output argument? It seems
really problematic conceptually and I couldn't think of any reason why
anyone would want it. (Isn't array_prepare specifically designed for
cases where the ufunc machinery creates the output array object and it
needs some initialization? That by definition isn't happening if someone's
passed an out= argument.)

On Wed, Jul 17, 2013 at 3:43 PM, Pauli Virtanen notifications@github.comwrote:

@njsmith https://github.com/njsmith: Doesn't the weird way dot chooses
the operation axes cause problems for gufunc core loop selection? I may be
wrong, though. OTOH, just allowing dot use the same mechanism as ufuncs
is not a big sin in my opinion.

Regarding the points:

• argument normalization: agreed
• NotImplemented instead of array_priority --- agreed, I like the
  idea of using the same mechanism as in Python's operator overrides. Maybe
  checking also the override method in the output argument should be done
  (which would also cover array_prepare)?
• Exclusive numpy_ufunc over array_wrap / array_prepare:
  agreed --- that mechanism frankly seems overly complicated.
• I would prefer having separate reduce etc. methods. Their
  implementation would for sparse matrices likely be quite different from
  call, so it would be simplest to make them separate.

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/3524#issuecomment-21117601
.

@njsmith: dot broadcasts the leading axes:

>>> np.dot(np.random.rand(2,3,4), np.random.rand(2,4,5)).shape
(2, 3, 2, 5)

For gufunc, this would be (2, 3, 5). There's an (unexposed) implementation for a matrix_product ufunc working like this in numpy.linalg._gufunc_linalg, waiting for polish, though.

I don't really have an use case for broadcasting on output, apart from probably misremembering how __array_prepare__. The docs seem to say it is called also on output, though?. Right now I don't again understand how it was supposed to work :)

I agree that standardizing the input argument format is good. But why standardize on dict? I think this could slow thing down. Why not use a tuple like (inputs..., outputs...)? Or 2 tuples, one for inputs and one for outputs?

I don't know about the other arguments you talked about.

Ufuncs can take tons of kwargs, can't avoid that. For call the current
list is: out, where, casting, dtype, sig, extobj, subok. These have to be
passed on somehow. And usually we receive all of these together packaged by
python into a dictionary. So dictionaries are pretty handy: in many cases
(when the user used the kwarg form, add(a, b, out=c) instead of add(a, b,
c)) they're actually the fastest option; they can easily represent missing
arguments (add(a, b)); they can easily be extended as new arguments are
added.

On Wed, Jul 17, 2013 at 4:21 PM, Frédéric Bastien
notifications@github.comwrote:

I agree that standardizing the input argument format is good. But why
standardize on dict? I think this could slow thing down. Why not use a
tuple like (inputs..., outputs...)? Or 2 tuples, one for inputs and one for
outputs?

I don't know about the other arguments you talked about.

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/3524#issuecomment-21120577
.

Okay, I'm rewriting the NEP to describe the implementation y'all are suggesting. I think __numpy_ufunc__ sounds better too. Thanks for the help.

@pv @njsmith It looks like ufunc methods (reduce, etc.) except for outer, take a different code path than the standard ufunc_generic_call. So I think this API can be told what the ufunc method being used is. Since I don't think there is a way to tell just being given the ufunc, args, and kwargs.

...Except handling outer might be a little different.

So the new function to check for an override can be like:

PyObject *
PyUFunc_CheckOverride(PyUFuncObject *ufunc, PyObject *ufunc_method, 
                      PyObject *args, PyObject *kwds)
{
...
}

Multi-methods are definitely the way to go here and should replace the __array_wrap__, __array_prepare__, __array_priority__ mechanism. Nathaniel's proposal to implement type-based multi-methods using the __numpy_ufunc__ attribute on the objects is thought-provoking and might be a good way to go, but it seems more complicated than other schemes and requires the introduction of another special method that doesn't seem really necessary. It seems more straight-forward to add the functions-to-be-called on the ufunc itself (via a decorator) and use a traditional covariant type hierarchy to resolve calls which are not handled exactly.

I recently wrote a blog-post where I linked to several common multi-method implementations that are worth looking at before going down this route of adding yet another special method to the objects. Here is the blog-post: http://technicaldiscovery.blogspot.com/2013/07/thoughts-after-scipy-2013-and-specific.html

I think looking into doing an actual multi-method implementation would be extremely beneficial for NumPy and would strongly encourage that whatever is done is at least done with awareness of the mulit-method literature such as this paper: http://dl.acm.org/citation.cfm?id=203096&dl=ACM&coll=

Traditional multi-methods are basically a simple look-up table (i.e. dictionary) whose keys are the types of the arguments (here actual Python types would work well). Obviously, if there is an exact match, you know what to do --- call the function and propagate errors. If the function called returns NotImplemented or if there is no such function registered for the types, then what? A simple approach is to follow-the method resolution order of the types (starting left-to-right on the input arguments only seems appropriate) to see if a function is registered for the super-types of the argument types.

One nice thing about the suggested __numpy_ufunc__ attribute, thought, is that it allows objects to "opt-out" of the __array_wrap__ business and provides a path-way to deprecating those attributes. This could perhaps also be accomplished by setting __array_wrap__ to a known, defined singleton object.

It is really great to see activity and progress towards building better multi-methods in NumPy.

@teoliphant: I think the __numpy_ufunc__ or __add__ approach is a limited form of multimethod dispatch, and these "Python-style multimethods" are actually superior to "classic multimethods" in a number of ways.

(1) They preserve Python's overall dispatch model. The classic multimethod literature all assumes that you're using a strongly-typed language with type-based dispatch. Python doesn't do that at all -- Python uses "duck typing", in which dispatch always goes directly through attributes of the object, without caring about the actual value of ->tp_class. (Of course many objects use the type hierarchy as a convenient way to modularize attribute lookup, but this is an implementation detail. Python "types" aren't types at all in the programming language sense, they're just a convenient way to create lots of similar-looking objects. Python isn't quite a prototype-based language like Self or Javascript, but you can definitely see the influence.) In Python IMO one should always feel a little dirty when calling isinstance explicitly, and classic multimethods are all about isinstance. Python-style multimethods, OTOH, are multimethods that dispatch on duck type.

(2) They preserve locality: if I see an expression a + b, and I want to trace down the code that's being executed, then barring monkeypatching, there are exactly two places I have to look -- a's class definition and b's class definition. With classic multimethods, then the search space is potentially unbounded -- any arbitrary module anywhere in the code could have registered a handler for this operation. Basically in classic multimethods, monkeypatching is considered to be a core feature instead of a nasty sin.

(3) They play better with having a fallback "other type" handler: For a numpy ufunc, there are lots of objects that we don't want to handle via any kind of multimethod dispatch, but instead want to fall back on using asarray to coerce them to ndarrays. It would be impossible to register multimethods on every new type implementing the buffer interface, for example. But what if some joker comes along and registers a handle for add(list, list)? (I'm sure they would have some good reason why this was useful in the context of their particular module! But the registration would of course be global.) With python-style multimethods, every type either does or doesn't have a handler for a particular multimethod, and this is immutable. So we know that builtins like lists will always hit the fallback, and so will any third-party buffer-interface object that didn't explicitly plan to interoperate with numpy, etc.

(4) They're familiar: While I love playing with cool designs for sexy programming language technology, I don't really think numpy should be in the business of trying to define a new dispatch system from scratch. But that's what we're doing as soon as we reach the hand-wavy bit of your design where we have to start walking multiple MRO hierarchies in parallel. (What happens when you have op(a, b), and there's no method registered for (type(a), type(b)), but there are methods for both (parent(type(a)), type(b)) and (type(a), parent(type(b)))? How do you decide which to use? Or can you disallow this ahead of time? Can you even detect such things ahead of time? These details matter, and the reason I'm using Python in the first place is because the language designers knew what they were doing.)

@njsmith A complication, what it an object defines __numpy_ufunc__ but wants default ufunc behavior for certain ufuncs?

I think it should probably return None. Then checking inputs would continue and at the end if only None has been returned, continue default ufunc behavior.

@teoliphant I think @njsmith's arguments are persuasive. However I'd like to hear what you think.

@cowlicks: if you want default ufunc behavior, you can cast your data to an ndarray and call the ufunc again. If you are a subclass and __numpy_ufunc__ disables array_wrap/prepare, the default ufunc behavior would very likely just produces something crazy (exhibit A: sparse matrices). So I'd rather not allow falling back to it.

@cowlicks: yeah, I had that kind of fallback behaviour in my initial sketch but then took it out before posting, because of what @pv says. With __numpy_ufunc__ It's still very easy to get the core ndarray ufunc behaviour if you want it. It's not so easy to get __array_wrap__ itself -- but that's a good thing. And it's still easy to do anything that __array_wrap__ would have let you do.

I'm concerned about the use of static variables in some of the functions. They are like globals and best avoided if possible.

Ping travisbot. Looks like pip is failing.

@cowlicks universal failure of build on account of pip.

Looks good modulo a few nitpicks.

Merged, thanks.