here is a snippet I extracted from a project with similar aims
(integrating the functionality of einsum and numexpr, actually)
Not much to it, but in case someone needs a reminder on how to use
striding tricks:
http://pastebin.com/kQNySjcj
On Fri, Aug 15, 2014 at 5:20 PM, Eelco Hoogendoorn
<hoogendoorn.ee...@gmail.com <mailto:hoogendoorn.ee...@gmail.com>> wrote:
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 <mailto: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
<mailto: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 <mailto: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
<mailto: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)
<http://en.wikipedia.org/wiki/Kronecker_delta%29>).
> >
> > (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
<mailto:NumPy-Discussion@scipy.org>
> >
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >
> >
> >
> > _______________________________________________
> > NumPy-Discussion mailing list
> > NumPy-Discussion@scipy.org
<mailto:NumPy-Discussion@scipy.org>
> >
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >
> >
> >
> > _______________________________________________
> > NumPy-Discussion mailing list
> > NumPy-Discussion@scipy.org
<mailto:NumPy-Discussion@scipy.org>
> >
http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org <mailto:NumPy-Discussion@scipy.org>
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org <mailto:NumPy-Discussion@scipy.org>
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org <mailto:NumPy-Discussion@scipy.org>
http://mail.scipy.org/mailman/listinfo/numpy-discussion
