[SciPy-user] += Operator and Slicing of Arrays
Josh Lawrence
josh.k.lawrence@gmail....
Fri Jun 12 15:29:31 CDT 2009
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
>> > http://www.scipy.org/Tentative_NumPy_Tutorial#head-0dffc419afa7d77d51062d40d2d84143db8216c2
>> > 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
>
--
Josh Lawrence
Ph.D. Student
Clemson University
More information about the SciPy-user
mailing list