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