Numpy array slicing vs. for-loop results

Question

Consider a sphere, composed of shells of varying density.
I have two arrays, one for the outer radius of each shell (rad[]) and one for the density of each shell (den[]). I want to calculate the mass, out to a given radius, called mass[].

The following for-loop approach achieves the desired result by first finding the mass of the innermost shell (the inner-radius is zero, so it's a sphere), then the mass of each subsequent shell - added to the previous (summed) mass:

mass = numpy.zeros(len(rad))                                   # create array
mass[0] = den[0]**(rad[0]**3)                                  # find inner sphere mass
for i in range(1,len(mass)):
    mass[i] = mass[i-1] + den[i]*(rad[i]**3 - rad[i-1]**3)     # Find mass out to shell i

Note: I only need the scalings, so I'm not worried about factors of pi.

Can anyone explain why the following slicing result does not achieve the same result?

mass = numpy.zeros(len(rad))
mass[0]  = den[0]*(rad[0]**3)
mass[1:] = mass[0:-1] + den[1:]*(rad[1:]**3-rad[0:-1]**3)

@IonutHulub: My answer was incorrect (or at least incomplete) so I deleted it, at least temporarily, until I can improve it. — David Robinson
– David Robinson, Commented Mar 10, 2013 at 0:41
It was correct. I was about to post the same thing when yours appeared. In the first case he modifies an element in mass and then uses that element to compute the next. in the second case the old value of the element is used to compute the new one, because the value of mass will be set only after the expression after the equals sign is evaluated. simple as that. — Ionut Hulub
– Ionut Hulub, Commented Mar 10, 2013 at 0:43
@IonutHulub: I wanted to explain why the example test = linspace(1,10,num=10); test[1:] += test[0:-1] works differently. (I am almost done with my revision to that effect). — David Robinson
– David Robinson, Commented Mar 10, 2013 at 0:46

David Robinson · Accepted Answer · 2013-03-10 00:51:26Z

3

The reason is that all the elements in the numpy array will be computed simultaneously. The array mass[0:-1] in your second line will be evaluated as den[0]*(rad[0]**3) followed by nothing but zeros. (The fact that mass[1] will no longer be zero once the line is calculated is irrelevant- by then it is too late).

You noted that the example:

test = numpy.linspace(1,10,num=10)
test[1:] += test[0:-1]
# [  1.   3.   6.  10.  15.  21.  28.  36.  45.  55.]

works differently, as though the addition does happen iteratively. The difference in your example is the addition of a value to the right side- that addition makes it a new array in memory (x + y is not the same array as x), such that numpy no longer treats it as adding to itself. See this example

test = numpy.linspace(1,10,num=10)
test[1:] += test[0:-1] + 0
# [  1.   3.   5.   7.   9.  11.  13.  15.  17.  19.]

If you want to do a vectorized version of your for loop, you can do:

mass[1:] += den[1:]*(rad[1:]**3-rad[0:-1]**3)
mass[1:] += mass[0:-1]

edited Mar 10, 2013 at 0:51

answered Mar 10, 2013 at 0:28

David Robinson

78.8k16 gold badges172 silver badges189 bronze badges

Sign up to request clarification or add additional context in comments.

9 Comments

DilithiumMatrix Over a year ago

That's what I guessed at first, but then I tried: test = linspace(1,10,num=10), test[1:] += test[0:-1] - and that correctly sums all the way to the given index

David Robinson Over a year ago

@zhermes: try changing your original code to mass[1:] += mass[0:-1] + den[1:]*(rad[1:]**3-rad[0:-1]**3).

David Robinson Over a year ago

@zhermes: Interesting. I'm looking into it.

seberg Over a year ago

Be careful with these in-place things, they are dangerous, and sometimes may appear to work when they don't. This one is probably fine, but just use cumsum or add.accumulate...

seberg Over a year ago

@zhermes I mean if you do in place operations on ufuncs np.add(a, b, out=c) or that a += c and a and c share memory. In this simple case it works, but numpy may buffer them, so there are some cases where there is no way to predict what happens.

|

Collectives™ on Stack Overflow

Numpy array slicing vs. for-loop results

1 Answer 1

9 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

9 Comments

Your Answer

Sign up or log in

Post as a guest

Related