<div dir="ltr"><div>Well, there is the numpy-API C level, and then there is the arcane macro C level. The two might as well be a completely different language.<br><br>Indeed, it should be doing something similar for the inputs. Actually, I think I wrote a wrapper around einsum/numexpr once that performed this generalized indexing once... ill see if I can dig that up.</div>
</div><div class="gmail_extra"><br><br><div class="gmail_quote">On Fri, Aug 15, 2014 at 5:01 PM, Sebastian Berg <span dir="ltr"><<a href="mailto:sebastian@sipsolutions.net" target="_blank">sebastian@sipsolutions.net</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div>On Fr, 2014-08-15 at 16:42 +0200, Eelco Hoogendoorn wrote:<br>
> Agreed; this addition occurred to me as well. Note that the<br>
> implemenatation should be straightforward: just allocate an enlarged<br>
> array, use some striding logic to construct the relevant view, and let<br>
> einsums internals act on the view. hopefully, you wont even have to<br>
> touch the guts of einsum at the C level, because id say that isn't for<br>
> the faint of heart...<br>
><br>
<br>
</div>I am not sure if einsum isn't pure C :). But even if, it should be doing<br>
something identical already for duplicate indices on the inputs...<br>
<span class="HOEnZb"><font color="#888888"><br>
- Sebastian<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
><br>
> On Fri, Aug 15, 2014 at 3:53 PM, Sebastian Berg<br>
> <<a href="mailto:sebastian@sipsolutions.net">sebastian@sipsolutions.net</a>> wrote:<br>
>         On Do, 2014-08-14 at 12:42 -0700, Stephan Hoyer wrote:<br>
>         > I think this would be very nice addition.<br>
>         ><br>
>         ><br>
>         > On Thu, Aug 14, 2014 at 12:21 PM, Benjamin Root<br>
>         <<a href="mailto:ben.root@ou.edu">ben.root@ou.edu</a>><br>
>         > wrote:<br>
>         >         You had me at Kronecker delta... :-)  +1<br>
>         ><br>
><br>
><br>
>         Sounds good to me. I don't see a reason for not relaxing the<br>
>         restriction, unless there is some technical issue, but I doubt<br>
>         that.<br>
><br>
>         - Sebastian<br>
><br>
>         ><br>
>         ><br>
>         >         On Thu, Aug 14, 2014 at 3:07 PM, Pierre-Andre Noel<br>
>         >         <<a href="mailto:noel.pierre.andre@gmail.com">noel.pierre.andre@gmail.com</a>> wrote:<br>
>         >                 (I created issue 4965 earlier today on this<br>
>         topic, and<br>
>         >                 I have been<br>
>         >                 advised to email to this mailing list to<br>
>         discuss<br>
>         >                 whether it is a good<br>
>         >                 idea or not. I include my original post<br>
>         as-is,<br>
>         >                 followed by additional<br>
>         >                 comments.)<br>
>         ><br>
>         >                 I think that the following new feature would<br>
>         make<br>
>         >                 `numpy.einsum` even<br>
>         >                 more powerful/useful/awesome than it already<br>
>         is.<br>
>         >                 Moreover, the change<br>
>         >                 should not interfere with existing code, it<br>
>         would<br>
>         >                 preserve the<br>
>         >                 "minimalistic" spirit of `numpy.einsum`, and<br>
>         the new<br>
>         >                 functionality would<br>
>         >                 integrate in a seamless/intuitive manner for<br>
>         the<br>
>         >                 users.<br>
>         ><br>
>         >                 In short, the new feature would allow for<br>
>         repeated<br>
>         >                 subscripts to appear<br>
>         >                 in the "output" part of the `subscripts`<br>
>         parameter<br>
>         >                 (i.e., on the<br>
>         >                 right-hand side of `->`). The corresponding<br>
>         dimensions<br>
>         >                 in the resulting<br>
>         >                 `ndarray` would only be filled along their<br>
>         diagonal,<br>
>         >                 leaving the off<br>
>         >                 diagonal entries to the default value for<br>
>         this `dtype`<br>
>         >                 (typically zero).<br>
>         >                 Note that the current behavior is to raise<br>
>         an<br>
>         >                 exception when repeated<br>
>         >                 output subscripts are being used.<br>
>         ><br>
>         >                 This is simplest to describe with an example<br>
>         involving<br>
>         >                 the dual behavior<br>
>         >                 of `numpy.diag`.<br>
>         ><br>
>         >                 ```python<br>
>         >                 # Extracting the diagonal of a 2-D array.<br>
>         >                 A = arange(16).reshape(4,4)<br>
>         >                 print(diag(A)) # Output: [ 0 5 10 15 ]<br>
>         >                 print(einsum('ii->i', A)) # Same as previous<br>
>         line<br>
>         >                 (current behavior).<br>
>         ><br>
>         >                 # Constructing a diagonal 2-D array.<br>
>         >                 v = arange(4)<br>
>         >                 print(diag(v)) # Output: [[0 0 0 0] [0 1 0<br>
>         0] [0 0 2<br>
>         >                 0] [0 0 0 3]]<br>
>         >                 print(einsum('i->ii', v)) # New behavior<br>
>         would be same<br>
>         >                 as previous line.<br>
>         >                 # The current behavior of the previous line<br>
>         is to<br>
>         >                 raise an exception.<br>
>         >                 ```<br>
>         ><br>
>         >                 By opposition to `numpy.diag`, the approach<br>
>         >                 generalizes to higher<br>
>         >                 dimensions: `einsum('iii->i', A)` extracts<br>
>         the<br>
>         >                 diagonal of a 3-D array,<br>
>         >                 and `einsum('i->iii', v)` would build a<br>
>         diagonal 3-D<br>
>         >                 array.<br>
>         ><br>
>         >                 The proposed behavior really starts to shine<br>
>         in more<br>
>         >                 intricate cases.<br>
>         ><br>
>         >                 ```python<br>
>         >                 # Dummy values, these should be<br>
>         probabilities to make<br>
>         >                 sense below.<br>
>         >                 P_w_ab = arange(24).reshape(3,2,4)<br>
>         >                 P_y_wxab = arange(144).reshape(3,3,2,2,4)<br>
>         ><br>
>         >                 # With the proposed behavior, the following<br>
>         two lines<br>
>         >                 should be equivalent.<br>
>         >                 P_xyz_ab = einsum('wab,xa,ywxab,zy->xyzab',<br>
>         P_w_ab,<br>
>         >                 eye(2), P_y_wxab,<br>
>         >                 eye(3))<br>
>         >                 also_P_xyz_ab = einsum('wab,ywaab->ayyab',<br>
>         P_w_ab,<br>
>         >                 P_y_wxab)<br>
>         >                 ```<br>
>         ><br>
>         >                 If this is not convincing enough, replace<br>
>         `eye(2)` by<br>
>         >                 `eye(P_w_ab.shape[1])` and replace `eye(3)`<br>
>         by<br>
>         >                 `eye(P_y_wxab.shape[0])`,<br>
>         >                 then imagine more dimensions and repeated<br>
>         indices...<br>
>         >                 The new notation<br>
>         >                 would allow for crisper codes and reduce the<br>
>         >                 opportunities for dumb<br>
>         >                 mistakes.<br>
>         ><br>
>         >                 For those who wonder, the above computation<br>
>         amounts to<br>
>         >                 $P(X=x,Y=y,Z=z|A=a,B=b) = \sum_w P(W=w|<br>
>         A=a,B=b) P(X=x|<br>
>         >                 A=a)<br>
>         >                 P(Y=y|W=w,X=x,A=a,B=b) P(Z=z|Y=y)$ with<br>
>         $P(X=x|A=a)=<br>
>         >                 \delta_{xa}$ and<br>
>         >                 $P(Z=z|Y=y)=\delta_{zy}$ (using LaTeX<br>
>         notation, and<br>
>         >                 $\delta_{ij}$ is<br>
>         >                 [Kronecker's<br>
>         ><br>
>          delta](<a href="http://en.wikipedia.org/wiki/Kronecker_delta)" target="_blank">http://en.wikipedia.org/wiki/Kronecker_delta)</a>).<br>
>         ><br>
>         >                 (End of original post.)<br>
>         ><br>
>         >                 I have been told by @jaimefrio that "The<br>
>         best way of<br>
>         >                 getting a new<br>
>         >                 feature into numpy is putting it in<br>
>         yourself." Hence,<br>
>         >                 if discussions<br>
>         >                 here do reveal that this is a good idea,<br>
>         then I may<br>
>         >                 give a try at coding<br>
>         >                 it myself. However, I currently know nothing<br>
>         of the<br>
>         >                 inner workings of<br>
>         >                 numpy/ndarray/einsum, and I have higher<br>
>         priorities<br>
>         >                 right now. This means<br>
>         >                 that it could take a long while before I<br>
>         contribute<br>
>         >                 any code, if I ever<br>
>         >                 do. Hence, if anyone feels like doing it,<br>
>         feel free to<br>
>         >                 do so!<br>
>         ><br>
>         >                 Also, I am aware that storing a lot of zeros<br>
>         in an<br>
>         >                 `ndarray` may not, a<br>
>         >                 priori, be a desirable avenue. However,<br>
>         there are<br>
>         >                 times where you have<br>
>         >                 to do it: think of `numpy.eye` as an<br>
>         example. In my<br>
>         >                 case of application,<br>
>         >                 I use such diagonal structures in the<br>
>         initialization<br>
>         >                 of an `ndarray`<br>
>         >                 which is later updated through an iterative<br>
>         process.<br>
>         >                 After these<br>
>         >                 iterations, most of the zeros will be gone.<br>
>         Do other<br>
>         >                 people see a use<br>
>         >                 for such capabilities?<br>
>         ><br>
>         >                 Thank you for your time and have a nice day.<br>
>         ><br>
>         >                 Sincerely,<br>
>         ><br>
>         >                 Pierre-André Noël<br>
>         ><br>
>          _______________________________________________<br>
>         >                 NumPy-Discussion mailing list<br>
>         >                 <a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
>         ><br>
>          <a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
>         ><br>
>         ><br>
>         ><br>
>         >         _______________________________________________<br>
>         >         NumPy-Discussion mailing list<br>
>         >         <a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
>         ><br>
>          <a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
>         ><br>
>         ><br>
>         ><br>
>         > _______________________________________________<br>
>         > NumPy-Discussion mailing list<br>
>         > <a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
>         > <a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
><br>
><br>
>         _______________________________________________<br>
>         NumPy-Discussion mailing list<br>
>         <a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
>         <a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
><br>
><br>
><br>
> _______________________________________________<br>
> NumPy-Discussion mailing list<br>
> <a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
> <a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
<br>
<br>
_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
<a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
</div></div></blockquote></div><br></div>