# [SciPy-user] += Operator and Slicing of Arrays

Emmanuelle Gouillart emmanuelle.gouillart@normalesup....
Fri Jun 12 17:26:11 CDT 2009

Hi Josh,

I forgot to tell you before, but this discussion should rather be on the
numpy-disussion mailing-list: if you have similar questions about numpy
arrays and indexing, try to post on the numpy list.

I didn't understand that you want to sum multidimensional instead of
scalars. I think it is impossible to use np.histogram for your case
(unless the arrays to add have indeed a (3,1) shape, in which case I
would just use three times np.histogram inside a for loop... Sometimes
you should just leave well enough alone instead of removing *ALL* for
loops).

I had a look at the source code of np.histogram to understand how the
function works, and it is actually fairly simple to use the same
algorithm for you case with multidimensional arrays as weights. Here is
an example below, you can try to adapt it to your data.

>>> a = np.array([ 3,  9,  9,  7,  5,  1,  9,  3,  5, 10])#indices
>>> w = np.ones((10, 2))#weights
>>> w[::2] = 2#unequal weights
>>> # The code below is inspired by the code of np.histogram
>>> bins = np.arange(12)
>>> sorting_index = np.argsort(a)
>>> sa = a[sorting_index]
>>> sw = w[sorting_index]
>>> cw = np.concatenate((np.array([0,0]).reshape((1,2)),
>>> sw.cumsum(axis=0)), axis=0)
>>> bin_index = np.r_[sa.searchsorted(bins[:-1], 'left'), \
sa.searchsorted(bins[-1], 'right')]
>>> n = cw[bin_index]
>>> n = np.diff(n, axis=0)
>>> a
array([ 3,  9,  9,  7,  5,  1,  9,  3,  5, 10])
>>> w
array([[ 2.,  2.],
[ 1.,  1.],
[ 2.,  2.],
[ 1.,  1.],
[ 2.,  2.],
[ 1.,  1.],
[ 2.,  2.],
[ 1.,  1.],
[ 2.,  2.],
[ 1.,  1.]])
>>> n
array([[ 0.,  0.],
[ 1.,  1.],
[ 0.,  0.],
[ 3.,  3.],
[ 0.,  0.],
[ 4.,  4.],
[ 0.,  0.],
[ 1.,  1.],
[ 0.,  0.],
[ 5.,  5.],
[ 1.,  1.]])

Cheers,

Emmanuelle

On Fri, Jun 12, 2009 at 04:29:31PM -0400, Josh Lawrence wrote:
> Emmanuelle,

> In the example I gave to begin with, tri_idx_plus has shape (8,) (if
> tri_idx_plus were a numpy array), edge_unknown has shape (8,..), and
> basis_p has shape (8,...). In practice, the shape for edge_unknown
> would be (8,1,n) and basis_p would have shape (8,3,1). Thus, weight =
> edge_unknown*basis_p would have a shape (8,3,n) in the example. Below
> is a better example.

> import numpy as np

> tri_idx_plus = np.array([0, 3, 6, 3, 2, 4, 1, 4])
> edge_unknown = np.random.rand(8,1,1) + 1.0j*np.random.rand(8,1,1)
> basis_p = np.random.rand(8,3,1)

> weights = edge_unknown * basis_p
> tri_quantity = np.zeros((20,3,1)) + 0j

> for i in range(tri_idx_plus.size):
>     tri_quantity[tri_idx_plus[i],:,:] += weights[i,:,:]

> I believe that does exactly what I want. Now, my desire is to get rid
> of the for loop. That is why I at first tried to do

> tri_quantity[tri_idx_plus,:,:] += weights

> but to no avail as only the last reference in tri_idx_plus to 3 and 4
> would be summed to tri_quantity.

> I hope this is more clear.

> Cheers,

> Josh Lawrence

> On Fri, Jun 12, 2009 at 4:04 PM, Emmanuelle
> Gouillart<emmanuelle.gouillart@normalesup.org> wrote:

> >> It tells me that the weights argument needs to be the same shape as
> >> tri_idx_plus. Also, edge_unknown is complex valued. So I have +- real
> >> and imaginary components. I do not know much about numpy.histogram...
> >> Does the complex dtype preclude it from use?

> >        Yes, the weight array must have the same shape as the values
> > array ("a" in the help of np.histogram). This is the case in the example
> > I provided. Can you check that the shapes of the two arrays are the same
> > in the code you execute? And if you can't find the origin of the error,
> > can you post a minimal example of code that reproduces the error?

> >        Also, as Anne mentioned (thanks for your answer, Anne!), the
> > weights keyword argument can be used with negative numbers, but also with
> > complex numbers (did you check it before posting?).
> >>>> a = np.arange(10)
> >>>> b = 1j*np.ones_like(a)
> >>>> np.histogram(a, a, weights=b)
> > (array([ 0.+1.j,  0.+1.j,  0.+1.j,  0.+1.j,  0.+1.j,  0.+1.j,  0.+1.j,
> >        0.+1.j,  0.+2.j]), array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]))
> > So you don't need to separate real and imaginary parts, you can use
> > directly your complex array as weights.

> >        Cheers,

> >        Emmanuelle

> >> Cheers,

> >> Josh Lawrence
> >> Ph.D. Student
> >> Clemson University

> >> On Jun 12, 2009, at 2:28 PM, Emmanuelle Gouillart wrote:

> >> > Hi Josh,

> >> > what kind of problem do you have exactly ? Do you have trouble
> >> > implementing the computation you describe, or do you get unexpected
> >> > results?

> >> > If I understood well what you want to do, you cannot use directly
> >> > use +=
> >> > with fancy indexing (quantity[tri_idx_plus,...]) because the repeated
> >> > elements will be incremented just once (see
> >> > for more details).

> >> > However, I think you can solve your problem by using a weighted
> >> > histogram. Using your notations

> >> > weights = edge_unknown[:len(tri_idx_plus),...] *
> >> >     basis_p[:len(tri_idx_plus),...]
> >> > histogram_values = np.histogram(tri_idx_plus,
> >> >     np.arange(tri_idx_plus.max() +2), weights=weights)
> >> > unique_plus = np.unique1d(tri_idx_plus)
> >> > quantity[unique_plus,...] = histogram_values[0][unique_plus]

> >> > The weighted histogram allows you to make the sums corresponding to
> >> > each
> >> > triangle.

> >> > Here is an example

> >> >>>> quantity = np.zeros(10)
> >> >>>> tri_idx_plus = np.array([0, 3, 6, 3, 2, 4, 1, 4])
> >> >>>> weights = 2*tri_idx_plus + 1
> >> >>>> histogram_values = np.histogram(tri_idx_plus,
> >> >        np.arange(tri_idx_plus.max() +2), weights=weights)
> >> >>>> unique_plus = np.unique1d(tri_idx_plus)
> >> >>>> quantity[unique_plus,...] = histogram_values[0][unique_plus]

> >> > Actually, this may only work with positive values of weights (not
> >> > checked)...

> >> > Please tell us if this meets your needs or not.

> >> > Cheers,

> >> > Emmanuelle

> >> > On Fri, Jun 12, 2009 at 01:14:16PM -0400, Josh Lawrence wrote:
> >> >> Greetings,

> >> >> I am in need of some help. I am trying to use the += operator to sum
> >> >> over edge elements on a triangular mesh. Each edge has an unknown
> >> >> associated with it. After solving for the unknowns, I am trying to
> >> >> compute a quantity at the centroid of all triangles in the mesh. On a
> >> >> flat surface, each edge will be connected to either one or two
> >> >> triangles. The orientation of the edge and the normals of the
> >> >> triangles determines whether each triangle attached to the edge is a
> >> >> "plus" or "minus" triangle for that edge. It is possible for one
> >> >> triangle to be referenced three times as a "plus" triangle, three
> >> >> times as a "minus" triangle or any combination of "plus" and
> >> >> "minus" (1 and 2 or 2 and 1, respectively).

> >> >> I have a variable tri_idx which relates the edges to the "plus" and
> >> >> "minus" triangles. I then compute the quantity at the centroid for
> >> >> the
> >> >> "plus" triangle and "minus" triangle attached to each edge. An
> >> >> example
> >> >> follows:

> >> >> tri_idx_plus = [0 3 6 3 2 4 1 4]
> >> >> tri_idx_minus = [1 2 5 3 6 0 1 4]

> >> >> quantity[tri_idx_plus,...] += edge_unknown[:len(tri_idx_plus),...] *
> >> >> basis_p[:len(tri_idx_plus),...]
> >> >> quantity[tri_idx_minus,...] += edge_unknown[:len(tri_idx_minus)...] *
> >> >> basis_m[:len(tri_idx_minus),...]

> >> >> where basis_p and basis_m are basis functions that expand the unknown
> >> >> of each edge into a surface function over the "plus" or "minus"
> >> >> triangle.

> >> >> I am pretty sure the problem I am encountering is that tri_idx_plus
> >> >> mentions indices 3 and 4 twice and tri_idx_minus contains index 1
> >> >> twice. Is there a way of doing this operation without reverting to
> >> >> looping over each edge (read: not doing this the slow way).

> >> >> Thanks in advance!

> >> >> Josh Lawrence
> >> >> Ph.D. Student
> >> >> Clemson University

> >> >> _______________________________________________
> >> >> SciPy-user mailing list
> >> >> SciPy-user@scipy.org
> >> >> http://mail.scipy.org/mailman/listinfo/scipy-user
> >> > _______________________________________________
> >> > SciPy-user mailing list
> >> > SciPy-user@scipy.org
> >> > http://mail.scipy.org/mailman/listinfo/scipy-user

> >> _______________________________________________
> >> SciPy-user mailing list
> >> SciPy-user@scipy.org
> >> http://mail.scipy.org/mailman/listinfo/scipy-user
> > _______________________________________________
> > SciPy-user mailing list
> > SciPy-user@scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user