<p dir="auto">In continuing to try to get astropy's <code class="notranslate">Quantity</code> to behave with <code class="notranslate">__numpy_ufunc__</code>, I found what probably is another bug.  One of our tests checks that an integer can be used as an <code class="notranslate">__index__</code> to, e.g., multiply a list, as in:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="import numpy as np
np.array([2]) * ['a', 'b']"><pre class="notranslate"><code class="notranslate">import numpy as np
np.array([2]) * ['a', 'b']
</code></pre></div>
<p dir="auto">(which yields <code class="notranslate">['a', 'b', 'a', 'b']</code>).</p>
<p dir="auto">But if I make a trivial class that implements <code class="notranslate">__numpy_ufunc__</code>, the <code class="notranslate">NotImplemented</code> does not seem to be passed on correctly:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="class MyA(np.ndarray):
    def __numpy_ufunc__(self, ufunc, method, i, inputs, **kwargs):
        result = getattr(ufunc, method)(*((input.view(np.ndarray)
                                           if isinstance(input, np.ndarray)
                                           else input) for input in inputs),
                                        **kwargs)
        print(&quot;Got result={}&quot;.format(result))
        return result

np.array([2]).view(MyA) * ['a', 'b']"><pre class="notranslate"><code class="notranslate">class MyA(np.ndarray):
    def __numpy_ufunc__(self, ufunc, method, i, inputs, **kwargs):
        result = getattr(ufunc, method)(*((input.view(np.ndarray)
                                           if isinstance(input, np.ndarray)
                                           else input) for input in inputs),
                                        **kwargs)
        print("Got result={}".format(result))
        return result

np.array([2]).view(MyA) * ['a', 'b']
</code></pre></div>
<p dir="auto">yields</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="Got result=NotImplemented
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
&lt;ipython-input-21-b1e259840c66&gt; in &lt;module&gt;()
----&gt; 1 np.array([2]).view(MyA) * ['a', 'b']

TypeError: __numpy_ufunc__ not implemented for this type."><pre class="notranslate"><code class="notranslate">Got result=NotImplemented
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
&lt;ipython-input-21-b1e259840c66&gt; in &lt;module&gt;()
----&gt; 1 np.array([2]).view(MyA) * ['a', 'b']

TypeError: __numpy_ufunc__ not implemented for this type.
</code></pre></div>
<p dir="auto">p.s. This may be related to <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="36045125" data-permission-text="Title is private" data-url="https://github.com/numpy/numpy/issues/4815" data-hovercard-type="pull_request" data-hovercard-url="/numpy/numpy/pull/4815/hovercard" href="https://github.com/numpy/numpy/pull/4815">#4815</a></p>
<p dir="auto">There is a complex set of questions around how to handle method resolution in the presence of <code class="notranslate">__numpy_ufunc__</code>. Currently in master is an extremely complicated set of rules that isn't documented and that I don't actually understand (see <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="66389458" data-permission-text="Title is private" data-url="https://github.com/numpy/numpy/issues/5748" data-hovercard-type="pull_request" data-hovercard-url="/numpy/numpy/pull/5748/hovercard" href="https://github.com/numpy/numpy/pull/5748">#5748</a> for the latest set of changes to this), so it's kinda hard to know whether they are correct, but I suspect not. And this is a blocker for 1.10, b/c whatever we release in 1.10 will be set in stone forever.</p>
<p dir="auto">I strongly feel that we cannot include <code class="notranslate">__numpy_ufunc__</code> in a release without at least having a document somewhere describing what the actual dispatch rules are. I hope that doesn't mean we have to defer <code class="notranslate">__numpy_ufunc__</code> for another release, but if it does then it does.</p>
<p dir="auto">AFAICT this is how <code class="notranslate">a op b</code> dispatch works for ndarrays, BEFORE <code class="notranslate">__numpy_ufunc__</code> (i.e., this is how 1.9 works):</p>
<ul dir="auto">
<li>First Python uses the subclass rule to decide whether to invoke <code class="notranslate">a.__op__(b)</code> or <code class="notranslate">b.__rop__(a)</code>. So in the case where one of these objects is a proper subclass of the other, that object always gets to do absolutely anything, so that's fine. The interesting cases are the ones where neither is a proper subclass of the other (either because it's like, <code class="notranslate">matrix + masked array</code>, or because it's like <code class="notranslate">ndarray + scipy.sparse</code>). So without loss of generality, let's focus on the case where Python calls <code class="notranslate">a.__op__(b)</code>, and <code class="notranslate">a</code> is either an instance of <code class="notranslate">ndarray</code> or else an instance of a subclass of <code class="notranslate">ndarray</code> which has not overridden <code class="notranslate">__op__</code>, i.e. we're getting <code class="notranslate">ndarray.__op__(a, b)</code>.</li>
<li><code class="notranslate">ndarray.__op__</code> has the following logic (see <code class="notranslate">PyArray_GenericBinaryFunction</code> in <code class="notranslate">number.c</code>):
<ul dir="auto">
<li>If <code class="notranslate">b</code> is not an <code class="notranslate">ndarray</code> at all (even a subclass), and <code class="notranslate">b</code> has a higher <code class="notranslate">__array_priority__</code> than <code class="notranslate">a</code>, then we return <code class="notranslate">NotImplemented</code> and let control pass to <code class="notranslate">b.__rop__(a)</code>.</li>
<li>Otherwise, we call <code class="notranslate">np.op(a, b)</code> and let the ufunc machinery take over.</li>
</ul>
</li>
<li><code class="notranslate">np.op(a, b)</code> does the following (see <code class="notranslate">PyUFunc_GenericFunction</code>, <code class="notranslate">PyUFunc_GeneralizedFunction</code>,  in <code class="notranslate">ufunc_object.c</code>, and also <code class="notranslate">ufunc_generic_call</code> which converts <code class="notranslate">-2</code> return values from the previous into <code class="notranslate">NotImplemented</code> so you have to audit their whole call stack):
<ul dir="auto">
<li>If <code class="notranslate">b</code> is not an <code class="notranslate">ndarray</code>, and calling <code class="notranslate">np.array(b)</code> returns an object array (presumably because coercion failed... though I guess this could also be hit if <code class="notranslate">b.__array__()</code> return an object array or something), AND <code class="notranslate">b</code> has a higher <code class="notranslate">__array_priority__</code> than <code class="notranslate">a</code>, and <code class="notranslate">b</code> has an <code class="notranslate">__rop__</code> method, then return <code class="notranslate">NotImplemented</code>.</li>
<li>If any of our arrays contain structured dtypes or strings, and there are no special struct ufunc loops registered, but not if any of our arrays contain objects, then return <code class="notranslate">NotImplemented</code>. (This is buried in <code class="notranslate">get_ufunc_arguments</code>, search for <code class="notranslate">return -2</code>.)</li>
<li>Otherwise we return the actual ufunc result.</li>
</ul>
</li>
</ul>
<p dir="auto">Now, my suggestion is that the way we would EVENTUALLY like this to look is:</p>
<ul dir="auto">
<li>First, Python uses the subclass rule to decide whether to invoke <code class="notranslate">a.__op__(b)</code> or <code class="notranslate">b.__rop__(a)</code>. As above, let's assume that it invokes <code class="notranslate">ndarray.__op__(a, b)</code>.</li>
<li><code class="notranslate">ndarray.__op__(a, b)</code> calls <code class="notranslate">np.op(a, b)</code> (which in turn invokes all the standard ufunc stuff, including <code class="notranslate">__numpy_ufunc__</code> resolution).</li>
<li>There is no step 3.</li>
</ul>
<p dir="auto">I submit that it is obvious that IF we can make this work, then it is obviously the ideal outcome, because it is the simplest possible solution. But is it too simple? To determine this we have to answer two questions: (1) Will it adequately address all the relevant use cases? (2) Can we get there from here?</p>
<p dir="auto">So let's compare the current rules to my dream rules.</p>
<p dir="auto">First, we observe that everything that currently happens inside the ufunc machinery looks like it's totally wrong. The first check can only be triggered if <code class="notranslate">b</code> is a non-<code class="notranslate">ndarray</code> that has a higher <code class="notranslate">__array_priority__</code> (among other things), but if we look above, we see that those conditions are sufficient to trigger the check in <code class="notranslate">ndarray.__op__</code>, so checking again at the ufunc level is redundant at best. And the second check is just incoherent nonsense AFAICT. The only reason to return <code class="notranslate">NotImplemented</code> is b/c you want to pass control to another <code class="notranslate">__(r)op__</code> method, and there's no reason arrays containing structured dtypes in particular should somehow magically have different <code class="notranslate">__(r)op__</code> methods available than other arrays. So we can just get rid of all the ufunc stuff immediately, great.</p>
<p dir="auto">That leaves the <code class="notranslate">__array_priority__</code> stuff. We have two problems here: we can't just drop this immediately b/c of backcompat issues, and we need to have some way to continue to support all the use cases that this currently supports. The first problem is just a matter of having a deprecation period. For the second, observe that a class which defines a <code class="notranslate">__numpy_ufunc__</code> method gets complete control over what any ufunc call does, so it has <em>almost</em> as much power as a class that currently sets <code class="notranslate">__array_priority__</code>. The only additional power that <code class="notranslate">__array_priority__</code> currently gives you is that it lets you distinguish between e.g. a call to <code class="notranslate">ndarray.__add(a, b)</code> versus a call to <code class="notranslate">np.add(a, b)</code>. So the only code that really loses out from my proposed change is code which wants <code class="notranslate">a + b</code> and <code class="notranslate">add(a, b)</code> to do different things.</p>
<p dir="auto">AFAIK in the entire history of numpy there is only one situation where this power has been used on purpose: the definition of matrix classes where <code class="notranslate">a * b</code> is matmul, but <code class="notranslate">np.multiply(a, b)</code> is elmul. And we've all agreed that such classes should be deprecated and eventually phased out (cite: PEP 465).</p>
<p dir="auto">So, I conclude that EVENTUALLY my dream rules should work great. The only problem is that we need some temporary compromises to get us from here to there. Therefore, I propose we use the following dispatch rules in numpy 1.10, with the goal of moving to my "dream rules" in some future version:</p>
<ul dir="auto">
<li>First, Python uses the subclass rule to decide whether to invoke <code class="notranslate">a.__op__(b)</code> or <code class="notranslate">b.__rop__(a)</code>. As above, let's assume that it invokes <code class="notranslate">ndarray.__op__(a, b)</code>.</li>
<li><code class="notranslate">ndarray.__op__(a, b)</code> does the following:
<ul dir="auto">
<li>If <code class="notranslate">b</code> <strong>does not define <code class="notranslate">__numpy_ufunc__</code> and</strong> is not an <code class="notranslate">ndarray</code> at all (even a subclass), and <code class="notranslate">b</code> has a higher <code class="notranslate">__array_priority__</code> than <code class="notranslate">a</code>, then we <strong>issue a deprecation warning and</strong> return <code class="notranslate">NotImplemented</code> and let control pass to <code class="notranslate">b.__rop__(a)</code>. (bolded parts are changes compared to the current behaviour)</li>
<li>If <code class="notranslate">__op__</code> is <code class="notranslate">__mul__</code> and <code class="notranslate">b-&gt;tp_class-&gt;tp_name.startswith("scipy.sparse.")</code>, then return <code class="notranslate">NotImplemented</code>. (This rule is necessary in addition to the above, because <code class="notranslate">scipy.sparse</code> has already made a release containing <code class="notranslate">__numpy_ufunc__</code> methods, so the exception above doesn't apply.)</li>
<li>Otherwise, we call <code class="notranslate">np.op(a, b)</code> and let the ufunc machinery take over.</li>
</ul>
</li>
</ul>
<p dir="auto">I believe that this is adequate to covers all practical use cases for the current dispatch machinery, and gives us a clean path to better dispatch machinery in the future.</p>
<p dir="auto">The main alternative proposal is Pauli's, which involves a very complicated check (I won't try to summarize here, see <a href="https://github.com/numpy/numpy/blob/eecb2e3c07f29c0ac991d364a846a2f8293a432a/numpy/core/src/multiarray/number.c#L103">this comment and following code</a>). The goal of that approach is to continue supporting classes where <code class="notranslate">a + b</code> and <code class="notranslate">add(a, b)</code> do different things. I don't think that keeping substantial additional complexity around indefinitely is worth it in order to support functionality that no-one has ever found a use for except in one very specific case (overriding <code class="notranslate">__mul__</code>), and where we generally agree that that one specific case should be phased out as possible.</p>
<p dir="auto">I would very much appreciate feedback from scipy.sparse and astropy in particular on whether the above covers all their concerns.</p>
<p dir="auto">(Partial) History: <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="36045125" data-permission-text="Title is private" data-url="https://github.com/numpy/numpy/issues/4815" data-hovercard-type="pull_request" data-hovercard-url="/numpy/numpy/pull/4815/hovercard" href="https://github.com/numpy/numpy/pull/4815">#4815</a>, <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="66389458" data-permission-text="Title is private" data-url="https://github.com/numpy/numpy/issues/5748" data-hovercard-type="pull_request" data-hovercard-url="/numpy/numpy/pull/5748/hovercard" href="https://github.com/numpy/numpy/pull/5748">#5748</a><br>
CC: <a class="user-mention notranslate" data-hovercard-type="user" data-hovercard-url="/users/pv/hovercard" data-octo-click="hovercard-link-click" data-octo-dimensions="link_type:self" href="https://github.com/pv">@pv</a>, <a class="user-mention notranslate" data-hovercard-type="user" data-hovercard-url="/users/cowlicks/hovercard" data-octo-click="hovercard-link-click" data-octo-dimensions="link_type:self" href="https://github.com/cowlicks">@cowlicks</a>, <a class="user-mention notranslate" data-hovercard-type="user" data-hovercard-url="/users/mhvk/hovercard" data-octo-click="hovercard-link-click" data-octo-dimensions="link_type:self" href="https://github.com/mhvk">@mhvk</a></p>