CODEWITHSUNDEEP
pythonarraysnumpyslice
DilithiumMatrix
DilithiumMatrix

Reputation: 18677

Numpy array slicing vs. for-loop results

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)

Upvotes: 3

Views: 1563

Answers (1)

David Robinson
David Robinson

Reputation: 78630

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]

Upvotes: 3

Related Questions