<div dir="ltr"><div class="gmail_quote"><div dir="ltr">On Mon, Jun 11, 2018 at 11:59 PM Eric Wieser <<a href="mailto:wieser.eric%2Bnumpy@gmail.com">wieser.eric+numpy@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="m_-6042199907354862978markdown-here-wrapper" style="font-size:1em;font-family:Helvetica,arial,freesans,clean,sans-serif;color:rgb(34,34,34);background-color:rgb(255,255,255);border:none;line-height:1.2"><blockquote style="margin:1em 0px;border-left:4px solid rgb(221,221,221);padding:0px 1em;color:rgb(119,119,119);quotes:none">
<p style="margin:1em 0px">I don’t understand your alternative here. If we overload np.matmul using <strong>array_function</strong>, then it would not use <em>ether</em> of these options for writing the operation in terms of other gufuncs. It would simply look for an <strong>array_function</strong> attribute, and call that method instead.  </p>
</blockquote>
<p style="margin:1em 0px">Let me explain that suggestion a little more clearly.</p>
<ol style="padding-left:2em;margin:1em 0px">
<li style="margin:1em 0px">There’d be a <code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">linalg.matmul2d</code> that performs the real matrix case, which would be easy to make as a ufunc right now.</li>
<li style="margin:1em 0px"><code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">__matmul__</code> and <code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">__rmatmul__</code> would just call <code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">np.matmul</code>, as they currently do (for consistency between <code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">np.matmul</code> and <code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">operator.matmul</code>, needed in python pre-<code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">@</code>-operator)</li>
<li style="margin:1em 0px"><code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">np.matmul</code> would be implemented as:<pre style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;font-size:1em;line-height:1.2em;overflow:auto;margin:1em 0px"><code class="m_-6042199907354862978hljs m_-6042199907354862978language-python" style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline;white-space:pre-wrap;border-radius:3px;border:1px solid rgb(204,204,204);padding:0.5em 0.7em;display:block;padding:0.5em;color:rgb(51,51,51);background:rgb(248,248,255)"><span class="m_-6042199907354862978hljs-decorator">@do_array_function_overrides</span>
<span class="m_-6042199907354862978hljs-function"><span class="m_-6042199907354862978hljs-keyword">def</span> <span class="m_-6042199907354862978hljs-title">matmul</span><span class="m_-6042199907354862978hljs-params">(a, b)</span>:</span>
    <span class="m_-6042199907354862978hljs-keyword">if</span> a.ndim != <span class="m_-6042199907354862978hljs-number">1</span> <span class="m_-6042199907354862978hljs-keyword">and</span> b.ndim != <span class="m_-6042199907354862978hljs-number">1</span>:
        <span class="m_-6042199907354862978hljs-keyword">return</span> matmul2d(a, b)
    <span class="m_-6042199907354862978hljs-keyword">elif</span> a.ndim != <span class="m_-6042199907354862978hljs-number">1</span>:
        <span class="m_-6042199907354862978hljs-keyword">return</span> matmul2d(a, b[:,<span class="m_-6042199907354862978hljs-keyword">None</span>])[...,<span class="m_-6042199907354862978hljs-number">0</span>]
    <span class="m_-6042199907354862978hljs-keyword">elif</span> b.ndim != <span class="m_-6042199907354862978hljs-number">1</span>:
        <span class="m_-6042199907354862978hljs-keyword">return</span> matmul2d(a[<span class="m_-6042199907354862978hljs-keyword">None</span>,:], b)
    <span class="m_-6042199907354862978hljs-keyword">else</span>:
        <span class="m_-6042199907354862978hljs-comment"># this one probably deserves its own ufunf</span>
        <span class="m_-6042199907354862978hljs-keyword">return</span> matmul2d(a[<span class="m_-6042199907354862978hljs-keyword">None</span>,:], b[:,<span class="m_-6042199907354862978hljs-keyword">None</span>])[<span class="m_-6042199907354862978hljs-number">0</span>,<span class="m_-6042199907354862978hljs-number">0</span>]
</code></pre>
</li>
<li style="margin:1em 0px"><code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">Quantity</code> can just override <code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">__array_ufunc__</code> as with any other ufunc</li>
<li style="margin:1em 0px"><code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">DataArray</code>, knowing the above doesn’t work, would implement something like<pre style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;font-size:1em;line-height:1.2em;overflow:auto;margin:1em 0px"><code class="m_-6042199907354862978hljs m_-6042199907354862978language-python" style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline;white-space:pre-wrap;border-radius:3px;border:1px solid rgb(204,204,204);padding:0.5em 0.7em;display:block;padding:0.5em;color:rgb(51,51,51);background:rgb(248,248,255)"><span class="m_-6042199907354862978hljs-decorator">@matmul.register_array_function(DataArray)</span>
<span class="m_-6042199907354862978hljs-function"><span class="m_-6042199907354862978hljs-keyword">def</span> <span class="m_-6042199907354862978hljs-title">__array_function__</span><span class="m_-6042199907354862978hljs-params">(a, b)</span>:</span>
    <span class="m_-6042199907354862978hljs-keyword">if</span> a.ndim != <span class="m_-6042199907354862978hljs-number">1</span> <span class="m_-6042199907354862978hljs-keyword">and</span> b.ndim != <span class="m_-6042199907354862978hljs-number">1</span>:
        <span class="m_-6042199907354862978hljs-keyword">return</span> matmul2d(a, b)
    <span class="m_-6042199907354862978hljs-keyword">else</span>:
        <span class="m_-6042199907354862978hljs-comment"># either:</span>
        <span class="m_-6042199907354862978hljs-comment"># - add/remove dummy dimensions in a dataarray-specific way</span>
        <span class="m_-6042199907354862978hljs-comment"># - downcast to ndarray and do the dimension juggling there</span>
</code></pre>
</li>
</ol>
<p style="margin:1em 0px">Advantages of this approach:</p>
<ul style="padding-left:2em;margin:1em 0px">
<li style="margin:1em 0px"><p style="margin:1em 0px">Neither the ufunc machinery, nor <code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">__array_ufunc__</code>, nor the inner loop, need to know about optional dimensions.</p>
</li>
<li style="margin:1em 0px"><p style="margin:1em 0px">We get a <code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">matmul2d</code> ufunc, that all subclasses support out of the box if they support <code style="font-size:0.85em;font-family:Consolas,Inconsolata,Courier,monospace;margin:0px 0.15em;padding:0px 0.3em;white-space:nowrap;border:1px solid rgb(234,234,234);background-color:rgb(248,248,248);border-radius:3px;display:inline">matmul</code></p>
</li>
</ul>
<p style="margin:1em 0px">Eric</p></div></div></blockquote><div>OK, this sounds pretty reasonable to me -- assuming we manage to figure out the __array_function__ proposal!</div><div><br></div><div>There's one additional ingredient we would need to make this work well: some way to guarantee that "ndim" and indexing operations are available without casting to a base numpy array.</div><div><br></div><div>For now, np.asanyarray() would probably suffice, but that isn't quite right (e.g., this would fail for np.matrix).</div><div> </div><div>In the long term, I think we need a new coercion protocol for "duck" arrays. Nathaniel Smith and I started writing a NEP on this, but it isn't quite ready yet. </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="m_-6042199907354862978markdown-here-wrapper" style="font-size:1em;font-family:Helvetica,arial,freesans,clean,sans-serif;color:rgb(34,34,34);background-color:rgb(255,255,255);border:none;line-height:1.2"><div title="MDH:Jmd0OyBJIGRvbid0IHVuZGVyc3RhbmQgeW91ciBhbHRlcm5hdGl2ZSBoZXJlLiBJZiB3ZSBvdmVy
bG9hZCBucC5tYXRtdWwgdXNpbmcgX19hcnJheV9mdW5jdGlvbl9fLCB0aGVuIGl0IHdvdWxkIG5v
dCB1c2UgKmV0aGVyKiBvZiB0aGVzZSBvcHRpb25zIGZvciB3cml0aW5nIHRoZSBvcGVyYXRpb24g
aW4gdGVybXMgb2Ygb3RoZXIgZ3VmdW5jcy4gSXQgd291bGQgc2ltcGx5IGxvb2sgZm9yIGFuIF9f
YXJyYXlfZnVuY3Rpb25fXyBhdHRyaWJ1dGUsIGFuZCBjYWxsIHRoYXQgbWV0aG9kIGluc3RlYWQu
wqDCoDxicj48YnI+TGV0IG1lIGV4cGxhaW4gdGhhdCBzdWdnZXN0aW9uIGEgbGl0dGxlIG1vcmUg
Y2xlYXJseS48ZGl2Pjxicj48L2Rpdj48ZGl2PjEuIFRoZXJlJ2QgYmUgYSBgbGluYWxnLm1hdG11
bDJkYCB0aGF0IHBlcmZvcm1zIHRoZSByZWFsIG1hdHJpeCBjYXNlLCB3aGljaCB3b3VsZCBiZSBl
YXN5IHRvIG1ha2UgYXMgYSB1ZnVuYyByaWdodCBub3cuPC9kaXY+PGRpdj4yLiBgX19tYXRtdWxf
X2AgYW5kIGBfX3JtYXRtdWxfX2Agd291bGQganVzdCBjYWxsIGBucC5tYXRtdWxgLCBhcyB0aGV5
IGN1cnJlbnRseSBkbyAoZm9yIGNvbnNpc3RlbmN5IGJldHdlZW4gYG5wLm1hdG11bGAgYW5kIGBv
cGVyYXRvci5tYXRtdWxgLCBuZWVkZWQgaW4gcHl0aG9uIHByZS1gQGAtb3BlcmF0b3IpPC9kaXY+
PGRpdj4zLiBgbnAubWF0bXVsYCB3b3VsZCBiZSBpbXBsZW1lbnRlZCBhczo8YnI+YGBgcHl0aG9u
PC9kaXY+PGRpdj5AZG9fYXJyYXlfZnVuY3Rpb25fb3ZlcnJpZGVzPGJyPmRlZiBtYXRtdWwoYSwg
Yik6PC9kaXY+PGRpdj4mbmJzcDsgJm5ic3A7IGlmIGEubmRpbSAhPSAxIGFuZCBiLm5kaW0gIT0g
MTo8L2Rpdj48ZGl2PiZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyByZXR1cm4gbWF0bXVsMmQo
YSwgYik8L2Rpdj48ZGl2PiZuYnNwOyAmbmJzcDsgZWxpZiBhLm5kaW0gIT0gMTo8L2Rpdj48ZGl2
PiZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyByZXR1cm4gbWF0bXVsMmQoYSwgYls6LE5vbmVd
KVsuLi4sMF08L2Rpdj48ZGl2PiZuYnNwOyAmbmJzcDsgZWxpZiBiLm5kaW0gIT0gMTo8L2Rpdj7C
oCDCoCDCoCDCoCByZXR1cm4gbWF0bXVsMmQoYVtOb25lLDpdLCBiKTxkaXY+Jm5ic3A7ICZuYnNw
OyBlbHNlOjwvZGl2PjxkaXY+Jm5ic3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7ICMgdGhpcyBvbmUg
cHJvYmFibHkgZGVzZXJ2ZXMgaXRzIG93biB1ZnVuZjwvZGl2PjxkaXY+Jm5ic3A7ICZuYnNwOyAm
bmJzcDsgJm5ic3A7IHJldHVybiBtYXRtdWwyZChhW05vbmUsOl0sIGJbOixOb25lXSlbMCwwXTwv
ZGl2PjxkaXY+YGBgPC9kaXY+PGRpdj40LiBgUXVhbnRpdHlgIGNhbiBqdXN0IG92ZXJyaWRlIGBf
X2FycmF5X3VmdW5jX19gIGFzIHdpdGggYW55IG90aGVyIHVmdW5jPC9kaXY+PGRpdj41LiBgRGF0
YUFycmF5YCwga25vd2luZyB0aGUgYWJvdmUgZG9lc24ndCB3b3JrLCB3b3VsZCBpbXBsZW1lbnQg
c29tZXRoaW5nIGxpa2U8L2Rpdj48ZGl2PmBgYHB5dGhvbjxicj5AbWF0bXVsLnJlZ2lzdGVyX2Fy
cmF5X2Z1bmN0aW9uKERhdGFBcnJheSk8YnI+ZGVmIF9fYXJyYXlfZnVuY3Rpb25fXyhhLCBiKTo8
L2Rpdj7CoCDCoCBpZiBhLm5kaW0gIT0gMSBhbmQgYi5uZGltICE9IDE6PGRpdj4mbmJzcDsgJm5i
c3A7ICZuYnNwOyAmbmJzcDsgcmV0dXJuIG1hdG11bDJkKGEsIGIpPC9kaXY+PGRpdj4mbmJzcDsg
Jm5ic3A7IGVsc2U6PC9kaXY+PGRpdj4mbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsgIyBlaXRo
ZXI6PC9kaXY+PGRpdj4mbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsgIyAtIGFkZC9yZW1vdmUg
ZHVtbXkgZGltZW5zaW9ucyBpbiBhIGRhdGFhcnJheS1zcGVjaWZpYyB3YXk8L2Rpdj48ZGl2PiZu
YnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAjIC0gZG93bmNhc3QgdG8gbmRhcnJheSBhbmQgZG8g
dGhlIGRpbWVuc2lvbiBqdWdnbGluZyB0aGVyZTwvZGl2PjxkaXY+PGRpdj5gYGA8L2Rpdj48ZGl2
Pjxicj48L2Rpdj48ZGl2PkFkdmFudGFnZXMgb2YgdGhpcyBhcHByb2FjaDo8L2Rpdj48ZGl2Piom
bmJzcDsgTmVpdGhlciB0aGUgdWZ1bmMgbWFjaGluZXJ5LCBub3IgYF9fYXJyYXlfdWZ1bmNfX2As
IG5vciB0aGUgaW5uZXIgbG9vcCwgbmVlZCB0byBrbm93IGFib3V0IG9wdGlvbmFsIGRpbWVuc2lv
bnMuPC9kaXY+PC9kaXY+PGRpdj4qIFdlIGdldCBhIGBtYXRtdWwyZGAgdWZ1bmMsIHRoYXQgYWxs
IHN1YmNsYXNzZXMgc3VwcG9ydCBvdXQgb2YgdGhlIGJveCBpZiB0aGV5IHN1cHBvcnQgYG1hdG11
bGA8L2Rpdj48ZGl2Pjxicj48L2Rpdj48ZGl2PkVyaWM8L2Rpdj48ZGl2Pjxicj48L2Rpdj48ZGl2
Pjxicj48L2Rpdj4=" style="height:0;width:0;max-height:0;max-width:0;overflow:hidden;font-size:0em;padding:0;margin:0">​</div></div></div>
</blockquote></div></div>