Supercharged Python: Take Your Code to the Next Level [ PART II ]

integer (giving the length of a one-dimensional array) or a tuple, describing N dimensions. The dtype is one of the values in Table 12.2. The order is either 'C' (row-major order, as in the C language) or 'F' (column-major order, as in FORTRAN). Note The name of this function is tricky, because the English word “zeros” can also be spelled “zeroes.” Remember to use the shorter spelling, zeros, only. Ah, English spelling—never mastered even by native English speakers! Here is a simple example creating a 3 × 3 two-dimensional array using the default float type. >>> import numpy as np >>> a = np.zeros((3,3)) >>> a array([[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]]) Here’s another example, this time creating a 2 × 2 × 3 array of integers. Click here to view code image >>> a = np.zeros((2, 2, 3), dtype=np.int16) >>> a array([[[0, 0, 0], [0, 0, 0]], [[0, 0, 0], [0, 0, 0]]], dtype=int16) Finally, here’s a one-dimensional array of Booleans. Notice that all the zero values are realized as the Boolean value False. Click here to view code image >>> a = np.zeros(5, dtype=np.bool) >>> a array([False, False, False, False, False])

12.5.8 The “full” Function The full function creates a numpy array using the same arguments shown earlier for empty, ones, and zeros; however, full has one additional argument: a value to assign to each and every element. Click here to view code image numpy.full(shape, fill_value, dtype=None, order='C') Briefly, shape is either an integer (giving the length of a one- dimensional array) or a tuple, describing N dimensions. The dtype is one of the values in Table 12.2. The order is either 'C' (row-major order, as in the C language) or 'F' (column- major order, as in FORTRAN). If the dtype argument is either omitted or set to None, the function uses the data type of the fill_value, which is required for this function. Here’s a simple example creating a 2 × 2 array in which each element is set to 3.14. >>> import numpy as np >>> a = np.full((2, 2), 3.14) >>> a array([[3.14, 3.14], [3.14, 3.14]]) Here’s another example; this one creates an array of eight integers, each set to 100. Click here to view code image

>>> a = np.full(8, 100) >>> a array([100, 100, 100, 100, 100, 100, 100, 100]) This final example takes advantage of the fact that you can create a numpy array of strings—provided that all these strings observe a fixed maximum size. Click here to view code image >>> a = np.full(5,'ken') >>> a array(['ken', 'ken', 'ken', 'ken', 'ken'], dtype='<U3') After the array is created in this way, with strings of size 3, each such string has, in effect, a maximum size. You can assign a longer string to any of these array elements, but it will be truncated. Click here to view code image a[0] = 'tommy' # Legal, but only 'tom' is assigned. 12.5.9 The “copy” Function The numpy copy function copies all the elements of an existing array. Because data is stored contiguously rather than through references, deep copying versus shallow copying is generally not an issue for numpy arrays. An example should suffice. Suppose you already have an array, a_arr, and you want to make a full copy of it and call the copy b_arr. import numpy as np b_arr = np.copy(a_arr)

12.5.10 The “fromfunction” Function The numpy fromfunction function (yes, that’s a mouthful) is among the most powerful ways to create an array; we’ll use it in the next section to create a multiplication table. The fromfunction enables you to create and initialize an array by calling another function that works as if it were transforming indexes into arguments. Click here to view code image numpy.fromfunction(func, shape, dtype='float') The shape is an integer or tuple of integers, just as with some of the other functions. The length of this tuple determines the rank of the array (the number of dimensions); it also determines how many arguments that the func—a callable— must accept. There’s a twist here in that shape must be a tuple and not a scalar, so you may have to use a tuple expression such as (5,) in order to create one-dimensional data sets. Here’s a simple example: We want to create a one- dimensional array corresponding to the first five natural numbers. You can do this with arange, but fromfunction provides another way. It requires us, however, to provide a callable. Click here to view code image import numpy as np def simple(n): return n + 1 a = np.fromfunction(simple, (5,), dtype='int32')

The resulting array is displayed as follows in IDLE: Click here to view code image array([1, 2, 3, 4, 5], dtype=int32) This might be better expressed with a lambda function. (See Chapter 3, “Advanced Capabilities of Lists,” for more information on lambdas.) Click here to view code image a = np.fromfunction(lambda n:n+1, (5,), dtype='int32') But higher-dimensional arrays are common. Here’s an example that creates a two-dimensional array in which each element is equal to the total of its two indexes. Click here to view code image def add_it(r, c): return r + c a = np.fromfunction(add_it, (3, 3), dtype='int32') This code could also be written with a lambda. Click here to view code image a = np.fromfunction(lambda r,c:r+c, (3, 3), dtype='int') In either case, if the resulting array is displayed in IDLE, it has the following representation: array([[0, 1, 2], [1, 2, 3], [2, 3, 4]])

At the beginning of this section, we stated that fromfunction works as if the function were being called for each element, with the arguments being the index or indexes at that position. What fromfunction actually does is create an array or arrays, in which one of the dimensions (or axis) is set to a series of whole numbers. For a one-dimensional array of size 6, this is the numbers 0 through 5. [0 1 2 3 4 5] This is an identity array in that each element is equal to its index. For the two-dimensional 3 × 3 array used in the previous example, fromfunction creates two arrays: one for each of the two axes. [[0 0 0], [1 1 1], [2 2 2]] [[0 1 2], [0 1 2], [0 1 2]] These are identity arrays along specific axes. In the first array, each element is equal to its row index; in the second array, each element is equal to its column index. The implementation of fromfunction operates on the arrays. As a result, the callable argument (the other function being called) is executed only once! But it is executed on one or more arrays—one for each dimension—enabling batch processing. If you use fromfunction the way it was designed to be used, this underlying implementation works. But if you do unorthodox things, strange results are possible. Consider the following code, which should produce a 3 × 3 array.

Click here to view code image a = np.fromfunction(lambda r, c: 1, (3, 3), dtype='int') You’d probably expect the result to be a 3 × 3 array in which each element is set to 1. Yet this function call (you can easily try it out) returns 1 as a scalar value! 12.6 EXAMPLE: CREATING A MULTIPLICATION TABLE Suppose you want to create the classic multiplication table for numbers from 1 to 10. There is more than one way to do that with numpy. You could create an empty array, for example, and assign values to the elements. With numpy, you could also do something similar. You could create an array initialized to all-zero values, for example, and then write a nested loop to assign R * C to each element— actually, it’s (R+1)*(C+1). By far the most efficient approach would be to use fromfunction to create an array that called a function to generate the values, without writing any loops at all. This is the numpy philosophy: As much as possible, let the package do all the work by using batch operations. You should be writing relatively few loops. Here’s the solution: Click here to view code image import numpy as np def multy(r, c): return (r + 1) * (c+ 1) a = np.fromfunction(multy, (10, 10), dtype=np.int16)

You could write this more compactly by using a lambda function. (Lambdas are explained in Chapter 3, “Advanced Capabilities of Lists.”) Click here to view code image a = np.fromfunction(lambda r,c: (r+1) * (c+1), (10, 10), dtype=np.int16) Printing the result, a, produces a nice-looking multiplication table. Click here to view code image >>> print(a) [[ 1 2 3 4 5 6 7 8 9 10] [ 2 4 6 8 10 12 14 16 18 20] [ 3 6 9 12 15 18 21 24 27 30] [ 4 8 12 16 20 24 28 32 36 40] [ 5 10 15 20 25 30 35 40 45 50] [ 6 12 18 24 30 36 42 48 54 60] [ 7 14 21 28 35 42 49 56 63 70] [ 8 16 24 32 40 48 56 64 72 80] [ 9 18 27 36 45 54 63 72 81 90] [ 10 20 30 40 50 60 70 80 90 100]] You can improve the appearance by getting rid of the brackets in the display. That’s relatively easy to do if we first convert to a string and then use the str class replace method. s = str(a) s = s.replace('[', '') s = s.replace(']', '') s=''+s As mentioned in Chapter 4, replacing a character with the empty string is a convenient way to purge all instances of a character. This example calls the replace method to get rid of both kinds of brackets. Finally, a space is inserted at the front of the string to make up for the loss of two open brackets.

Now, printing this string produces a pleasing display. Click here to view code image >>> print(s) 1 2 3 4 5 6 7 8 9 10 2 4 6 8 10 12 14 16 18 20 3 6 9 12 15 18 21 24 27 30 4 8 12 16 20 24 28 32 36 40 5 10 15 20 25 30 35 40 45 50 6 12 18 24 30 36 42 48 54 60 7 14 21 28 35 42 49 56 63 70 8 16 24 32 40 48 56 64 72 80 9 18 27 36 45 54 63 72 81 90 10 20 30 40 50 60 70 80 90 100 12.7 BATCH OPERATIONS ON “NUMPY” ARRAYS The real power and speed of a numpy array become clear when you start to do large-scale, or batch, operations on the array— either on the whole array or on selected pieces created by slicing. This lets you operate on selected rows and columns, or even intersections. Once you’ve created a numpy array, you can do any number of arithmetic operations on it by combining it with a scalar value. Table 12.3 lists some of the things you can do. This is far from an exhaustive list. In this table, A is a numpy array, and n is a scalar value, such as a single integer or floating-point number. Table 12.3. Some Scalar Operations on “numpy” Arrays Operation Produces A + Array with n added to each element of A.

n A - Array with n subtracted from each element of A. n A * Array with n multiplied with each element of A. n n A number raised to the power of each element of A, producing ** another array with the results. A A Each element in A is raised to the power n. ** n A / Array with n dividing into each element of A. n A Array with n dividing into each element of A but using ground // division. n Each of these operations has an assignment operator associated with it, as do ordinary Python operations. For example, to double each member of a numpy array named my_array, use the following: Click here to view code image my_array *= 2 # Double each element of my_array

Another simple, and very powerful, version of numpy batch operations is to use arithmetic on two numpy arrays of the same shape—which implies that the number of dimensions match, as well as the size of each dimension. Table 12.4 shows some of these operations. Table 12.4. Some Array-to-Array Operations Operation Produces A Array generated by adding each element of A to the corresponding + element of B. B A Array generated by subtracting each element of B from the - corresponding element of A. B A Array generated by multiplying each element of A with the * corresponding element of B. B A Array generated by raising each element in A to the exponent ** taken from the corresponding element in B. B A Array generated by dividing each element of B into the / corresponding element of A. B A Array generated by dividing each element of B into the // corresponding element of A, but using ground division. B

For example, let’s start with a simple 4 × 4 array. Click here to view code image import numpy as np A = np.array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11], [12, 13, 14, 15]]) print(A) Printing A gives a nice result. [[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11] [12 13 14 15]] This is a familiar pattern. Is there a way to produce such an array without having to enter all the data directly? Yes, there are at least two! The simplest technique is to generate the numbers from 1 to 15 to form a simple array, and then use the numpy reshape function to rearrange it into a 4 × 4 array with the same elements. Click here to view code image A = np.arange(16).reshape((4,4)) The other way (which is slightly more verbose) is to use fromfunction to do the job. In either case, you can apply this pattern to much, much larger arrays, with a shape such as 200 × 100 or even 1,000 × 3,000. In the case of creating a 4 × 4 array, the function call would be Click here to view code image A = np.fromfunction(lambda r, c: r*4 + c, (4, 4))

Suppose, also, that we have an array named B of matching shape and size. B = np.eye(4, dtype='int16') print(B) These statements print the following: [[ 1 0 0 0] [ 0 1 0 0] [ 0 0 1 0] [ 0 0 0 1]] Now the fun begins. We can produce a new array, for example, by multiplying every element in A by 10. From within IDLE, here’s what that looks like: >>> C = A * 10 >>> print(C) [[ 0. 10. 20. 30.] [ 40. 50. 60. 70.] [ 80. 90. 100. 110.] [120. 130. 140. 150.]] To produce the array now referred to by the variable C, each and every element of A has been multiplied by 10. We might also produce an array that contains the squares of all the elements in A. We can do this by multiplying A by itself, which does a member-by-member multiplication of each element. >>> C = A * A >>> print(C) [[ 0. 1. 4. 9.] [ 16. 25. 36. 49.] [ 64. 81. 100. 121.] [144. 169. 196. 225.]] Keep in mind there’s no requirement that numpy arrays have a perfect square or perfect cubic shape—only that they’re

rectangular. You can always reshape an array. For example, the 4 × 4 array just shown can be reshaped into a 2 × 8 array. Click here to view code image >>> print(C.reshape((2,8))) [[ 1. 4. 9. 16. 25. 36. 49. 64.] [ 81. 100. 121. 144. 169. 196. 225. 256.]] If we want to change A in place, that’s entirely doable with an assignment operator, *=. Arrays are mutable. >>> A *= A Finally, let’s assume that the square numbers are assigned to A itself using the statement just shown (A *= A). The next operation multiplies A by B. Because B is an eye (or identity) array, what do you think the result is? >>> C = A * B >>> print(C) [[ 0. 0. 0. 0.] [ 0. 25. 0. 0.] [ 0. 0. 100. 0.] [ 0. 0. 0. 225.]] Again, the result is a member-by-member multiplication. 12.8 ORDERING A SLICE OF “NUMPY” You can take slices of a one-dimensional numpy array just as you can with a Python list; the next section deals with higher- dimensional arrays. Given a numpy array, you can print a slice of it, just as you can with a Python list. Here’s an example:

>>> A = np.arange(1, 11) 9 10] >>> print(A) [1 2 3 4 5 6 7 8 >>> print(A[2:5]) [3 4 5] One of the interesting things you can do with a numpy slice is to assign a scalar value to it. The result is to assign this same value to each position in the slice. >>> A[2:5] = 0 >>> print(A) [ 1 2 0 0 0 6 7 8 9 10] And you can do more. You can operate on a slice of an array just as you can with a full array. If you use assignment, such operations happen in place. For example, you could add 100 to each of these three elements, rather than setting them to zero. >>> A[2:5] += 100 >>> print(A) [ 1 2 103 104 105 6 7 8 9 10] Remember that when you combine arrays together through standard operations, the size of the two arrays must match. This applies to slices as well. For example, the following is valid, because the shapes match. A[2:5] *= [100, 200, 300] The effect is to multiply the third, fourth, and fifth elements of A by 100, 200, and 300, respectively. That operation produces the following array (assuming we apply it to the original value of A): [ 1 2 300 800 1500 6 7 8 9 10]

Now, how can we use some of these features to solve practical problems? One of the classic benchmarks is the algorithm known as the Sieve of Eratosthenes, which is an efficient way of producing a large group of prime numbers. Let’s start with the numbers between 0 and 50, inclusive. The procedure (which we’ll generalize later) is to eliminate all the numbers that are not prime and then print the ones left. First, here’s the array we start with. Click here to view code image >>> A = np.arange(51) >>> print(A) [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50] We want to zero out all the numbers that are not prime. Zero out element A[1], because 1 is not a prime number. Zero out all multiples of 2, starting with 2 squared. Zero out all multiples of 3, starting with 3 squared. Repeat the procedure for 5 and 7. The following code carries out the first two steps. Click here to view code image >>> A[1] = 0 >>> A[2 * 2::2] = 0 >>> print(A) [ 0 0 2 3 0 5 0 7 0 9 0 11 0 13 0 15 0 17 0 19 0 21 0 23 0 25 0 27 0 29 0 31 0 33 0 35 0 37 0 39 0 41 0 43 0 45 0 47 0 49 0] The meaning of A[2 * 2::2] is to take a slice beginning with the index number 2 squared, which is 4, then go to the end of the array (because the middle argument is blank), and then

move through the array two elements at a time. Everything in this slice is set to 0. Notice that each index position, in this particular example, corresponds to the value of a number from 0 to 50 in this case. So to zero out the number 8, for example, we set A[8] to zero. This keeps the programming simple. The results show that A[1] has been zeroed out, as well as all the even numbers other than 2 itself. We can do the same thing for multiples of 3. Click here to view code image >>> A[3 * 3::3] = 0 0 11 0 13 0 0 0 [0 0 2 3 0 5 0 7 0 0 0 0 29 0 31 0 0 17 0 19 0 0 0 23 0 25 0 0 0 0 47 0 49 0 35 0 37 0 0 0 41 0 43 0] After repeating the procedure for multiples of 5 and 7, we finally get an array in which all the values are either 0 or are prime numbers. Click here to view code image [ 0 0 2 3 0 5 0 7 0 0 0 11 0 13 0 0 0 17 0 19 0 0 0 23 0 0 0 0 0 29 0 31 0 0 0 0 0 37 0 0 0 41 0 43 0 0 0 47 0 0 0] Now, how do we print all the nonzero values? You could write a loop, of course, that goes through the array and either prints a value if it is nonzero, or adds it to a list. Here’s an example: Click here to view code image my_prime_list = [i for i in A if i > 0] That’s not bad, but numpy provides a way that’s even more efficient and more compact! You can create a Boolean array just by specifying a condition.

A>0 The Boolean array that gets produced, which we’ll look at more closely in Section 12.10, can be applied to array A itself as a mask—just by using the indexing operation. The effect in this case is to produce a new array from A, satisfying the condition that the element is greater than 0. We previously zeroed out all nonprime numbers in A; therefore, by taking the nonzero values remaining in A, we get only prime numbers as a result. Click here to view code image >>> P = A[A > 0] >>> print(P) [ 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47] So there you have it: all the primes not larger than 50. 12.9 MULTIDIMENSIONAL SLICING numpy arrays provide an even more powerful slicing ability: getting slices of any number of dimensions supported by the source array. We can start by seeing the effect of an (apparently) one-dimensional slice of a two-dimensional array. Let’s start with a familiar 4 × 4 array. Click here to view code image >>> A = np.arange(1,17).reshape((4,4)) >>> print(A) [[ 1 2 3 4] [ 5 6 7 8] [ 9 10 11 12] [13 14 15 16]]

What would happen if we took the two middle elements—1 and 2? >>> print(A[1:3]) [[ 5 6 7 8] [ 9 10 11 12]] The result, clearly, is to produce the middle two rows. Now, how can we get the middle two columns? Actually, that turns out to be almost as easy. >>> A[:, 1:3] array([[ 2, 3], [ 6, 7], [10, 11], [14, 15]]) Take another look at that array expression. A[:, 1:3] The colon before the comma says, “Select everything in this dimension”—in this case, the row dimension. The expression 1:3 selects all the columns beginning with index 1 (the second column) up to but not including index 3 (the fourth column). Therefore, the expression says, “Select all rows, with the intersection of columns 1 up to but not including column 3— that is, the second and third columns.” The general syntax for indexing and slicing an array of N dimensions is shown here. array_name[ i1, i2, i3,... iN ]

In this syntax, each of the arguments i1 through iN may be either a scalar value—which must be an index in range—or a slice. You can use at most N such arguments, where N is the number of dimensions (the rank) of the array. For each scalar used, the dimensionality of the result is reduced by one. Therefore, slicing a two-dimensional array as A[2, 1:4] produces a one-dimensional array as a result. Slicing it as A[2:3, 1:4] would get the same elements but would be two- dimensional, even though it would have only one row that was not empty. (This issue matters, because most operations on arrays must match in size and number of dimensions.) Any of the i values in this syntax may be omitted; in each such case, its value is assumed to be the colon (:), which says, “Select all the elements in this dimension.” If there are fewer than N arguments, then the first M dimensions (where M is the number of arguments) get the values assigned, and the last N– M dimensions assume the colon as their default. Table 12.5 lists a series of examples. In this table, A is a two- dimensional array and A3D is a three-dimensional array. Table 12.5. Examples of “numpy” Indexing and Slicing Example Description A[3] The entire fourth row, returned as a one-dimensional array. A[3,: Same as above. ] A[3,] Same as above. A[:,2 The entire third column, returned as a one-dimensional array. ]

A[::2 Get every other row from the third column. ,2] A[1:3 The intersection of the second and third columns with the ,1:3] second and third rows, returned as a two-dimensional array. A3D[2 The third row of the third plane (a plane being a level in the ,2] third dimension), returned as a one-dimensional array. Takes all columns in that row. A3D[2 Same as above. , 2, :] A3D[: A two-dimensional array containing the third column from , the second and third rows, intersecting with all planes. 1:3, 2] A3D[: Same as above, but get every other plane (as opposed to every :2, plane), starting with the first. 1:3, 2] A3D[0 A single element, taken from second column of the first row of , 0, the first plane. 1] Let’s look at a more practical example. Suppose you’re writing a computer simulation called the Game of Life, and you have the following grid, realized as a numpy array that we’ll call G, for “grid.” The 1’s are in bold for clarity.

[[0 0 0 0 0 0] [0 0 1 0 0 0] [0 0 1 0 0 0] [0 0 1 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0]] A 1 represents a living cell; a 0 represents a dead cell. You’d like to get a count of all the neighbors around a specific cell, say, G[2, 2]. A fast way to do that is to get the two-dimensional slice that includes the columns before, after, and including 2, intersecting with the rows before, after, and including 2. >>> print(G[1:4, 1:4]) [[0 1 0] [0 1 0] [0 1 0]] Remember that the index numbers used are 1 and 4, and not 1 and 3, because the slice expressions always mean up to but not including the endpoint. This gives us a nice cross-section of all the cells neighboring G[2,2], not including the cell itself. To get the neighbor count, therefore, it’s only necessary to sum this cross-section and then subtract the value of the cell itself. Click here to view code image neighbor_count = np.sum(G[1:4, 1:4]) - G[2, 2] The result is 2. In the Game of Life, that would indicate that the cell in the middle is “stable”: In the next generation, it will experience neither a birth nor a death event. 12.10 BOOLEAN ARRAYS: MASK OUT THAT “NUMPY”!

We’ve already shown a use for Boolean arrays used as masks. Section 12.7, “Batch Operations on ‘numpy’ Arrays,” used the following expression at the end: A>0 Assuming that A is a numpy array, this expression means “For each element of A, produce an element that is True if the element in A is greater than 0, and False if the element is not greater than 0.” The resulting array will have the same shape as A. For example, start with an array named B: B = np.arange(1,10).reshape(3,3) B, when printed, looks like this: [[1 2 3] [4 5 6] [7 8 9]] Now let’s apply the condition B > 4. B1 = B > 4 The result is to produce a Boolean array, B1, with the following contents: [[False False False] [False True True] [ True True True]] B1 has the shame shape as B, but each element is either True or False rather than being an integer. The general rule is as follows.

Whenever a comparison operator (such as ==, <, or >) is applied to a numpy array, the result is a Boolean array of the same shape. One way to use this array is to zero out all elements that don’t meet the condition of being greater than 4, by multiplying the two arrays—B and (B>4)—together. >>> print(B * (B > 4)) [[0 0 0] [0 5 6] [7 8 9]] When working with Boolean arrays, you should note that the use of parentheses is often critical, because comparison operators have low precedence. But an even better use of a Boolean array is to use it as a mask—in which case it selects in elements with a corresponding True value in the mask, and selects out elements with a corresponding False value. Using a Boolean array as a mask produces a one-dimensional array, regardless of the shape of the array operands. Click here to view code image array_name[bool_array] # Use bool_array as a mask. For example, we could use a mask to get all elements greater than 7. The result is a one-dimensional array containing 8 and 9.

>>> print(B[B > 7]) [8 9] Here’s a more sophisticated use: Get all elements whose remainder, when divided by 3, is 1. There are three elements of B that meet this condition: 1, 4, and 7. >>> print(B[B % 3 == 1]) [1 4 7] A difficulty arises when you try to introduce complex conditions involving and and or keywords, even though these operations should work with Booleans. A good solution is to apply the bitwise operators (&, |) to the Boolean masks. The & symbol performs bitwise AND, whereas the | symbol performs bitwise OR. You can also use multiplication (*) and addition (+) to get the same results. For example, to create a Boolean array in which the test for each element of B is “greater than 2 and less than 7,” you could use the following: Click here to view code image B2 = (B > 2) & (B < 7) # \"AND\" operation Let’s break this down semantically. B is a two-dimensional array of integers. B > 2 is a Boolean array of the same shape as B. B < 7 is another Boolean array, again of the same shape as B. The expression (B > 2) & (B < 7) uses binary AND (&) to achieve an “and” effect between these two Boolean arrays. The resulting Boolean array is assigned to the variable B2. This array will contain True and False values, in effect produced by Boolean operations on the two arrays which serve as operands.

You can then apply the mask to B itself to get a one- dimensional array of results in which each element is greater than 2 and less than 7. >>> print(B[B2]) [3 4 5 6] In this next example, bitwise OR is used to create a Boolean array from an “or” operation. That resulting Boolean array is then applied as a mask to B, and the final result selects all elements of B that are either equal to 1 or greater than 6. Click here to view code image >>> print(B[ (B == 1) | (B > 6)]) # \"OR\" operation [1 7 8 9] 12.11 “NUMPY” AND THE SIEVE OF ERATOSTHENES Let’s return to the Sieve of Eratosthenes example and benchmark how numpy does compared to a standard Python list solution. The goal of the algorithm is to produce all the prime numbers up to N, where N is any number you choose in advance. Here’s the algorithm. Create a one-dimensional Boolean array indexed from 0 to N. Set all elements to True, except for the first two elements, which are set to False.

For I running from 2 to N: If array[I] is True, For J running from I*I to N, by steps of I: Set array[J] to False The result of these steps is a Boolean array. For each index number higher than 2, corresponding to a True element, add that index number to results. Here is an obvious way to implement this algorithm as a Python function: Click here to view code image def sieve(n): b_list = [True] * (n + 1) for i in range(2, n+1): if b_list[i]: for j in range(i*i, n+1, i): b_list[j] = False primes = [i for i in range(2, n+1) if b_list[i]] return primes Can we do better with numpy? Yes. We can improve performance by taking advantage of slicing and Boolean masking. In keeping with the general flavor of the algorithm, we use an array of Booleans, indexed from 2 to N–1. Click here to view code image import numpy as np def np_sieve(n): # Create B, setting all elements to True. B = np.ones(n + 1, dtype=np.bool) B[0:2] = False for i in range(2, n + 1): if B[i]: B[i*i: n+1: i] = False return np.arange(n + 1)[B]

So where does this implementation of the algorithm manage to do better? The function still has to loop through the array, one member at a time, looking for each element with the value True. This indicates that the index number is prime, because its corresponding element in the Boolean array has not been eliminated yet. But the inner loop is replaced by a slice operation that sets every element in the slice to False. Assuming there are many elements, we can perform all these operations more efficiently with a batch operation rather than a loop. B[i*i: n+1: i] = False The other advanced technology used here is a Boolean mask to produce the final results: a numpy array from 0 to n, inclusive, which after the masking operation will contain only the prime numbers in that range. return np.arange(n + 1)[B] Now, we’d like to know about the performance of this operation. By using the time package, you can perform benchmark tests designed to show which approach is faster. The following code adds lines that report the number of milliseconds taken. The added lines are in bold. The return statement is omitted, because you probably don’t need to print all the primes up to 1 million, for example, if you’re only interested in speed. Click here to view code image import numpy as np import time def np_sieve(n): t1 = time.time() * 1000 B = np.ones(n + 1, dtype=np.bool) B[0:2] = False

for i in range(2, n + 1): if B[i]: B[i*i: n+1: i] = False P = np.arange(n + 1)[B] t2 = time.time() * 1000 print('np_sieve took', t2-t1, 'milliseconds.') You can put in similar lines of timing code to benchmark the non-numpy version. What the benchmarks show is that for relatively small numbers, the numpy version takes more time, and not less, than the other version. But for large N, especially greater than 1,000, np_sieve starts pulling ahead. Once N gets greater than 10,000 or so, the numpy version takes half the time the other version odes. That may not be the spectacular results we were looking for, but it’s an increase in speed of 100 percent. Not bad. Is this section playing fair? Yes. It’s admittedly possible to implement the non-numpy version, sieve, by using more lists and more list comprehension. However, we’ve found that such attempts at code enhancement actually make the function run more slowly. Therefore, for large N, the numpy version remains the high-speed champ. 12.12 GETTING “NUMPY” STATS (STANDARD DEVIATION) One of the areas in which numpy excels is getting statistics on large data sets. Although you could get any of this information yourself by working with standard Python lists, this is where numpy arrays are many times faster. Table 12.6 lists the statistical-analysis functions for numpy arrays. Each of these works by calling the corresponding method for the ndarray class; so you can use either the function or the method version.

These functions have a series of important arguments, which we’ll cover later. Table 12.6. Statistical Functions for “numpy” Arrays Function Returns m The lowest element in the data set. Will return the element lowest i along each dimension if the axis argument is specified; the same n applies to each function listed here. ( A ) m The highest element. a x ( A ) m The arithmetic mean, which is the sum of the elements divided by e the number of elements. When applied to an individual axis (as a you’ll see in the next section), it will sum and divide along rows or n columns as appropriate. ( A ) m The median element, which is the element in the group having an e equal number of elements higher and lower. d i a n (

A ) s The number of elements. i z e ( A ) s Standard deviation, a classic measure of variance. t d ( A ) s The sum of all elements in the data set, or the sum of all elements in u the specified subset. m ( A ) Let’s start by looking at how these functions apply to simple one-dimensional arrays. As the next section shows, you have more choices when applying them to higher-dimensional arrays. The performance speeds of the statistical functions are impressive, as you’ll see. First, let’s generate an array to operate on. We can use a subpackage of numpy to generate random numbers—specifically, the numpy.random function, rand, which takes an array shape as input and generates an array of

that shape. Each element is a random floating-point number between 0.0 and 1.0. import numpy as np import numpy.random as ran A = ran.rand(10) Printing A produces an array full of random floating-point values. Click here to view code image [0.49353738 0.88551261 0.69064065 0.93626092 0.17818198 0.16637079 0.55144244 0.16262533 0.36946706 0.61859074] The numpy package can handle much bigger data sets and handle them efficiently, such as the following large array. But note: Don’t print this unless you want Python to be busy for a long, long time! A = ran.rand(100000) This statement creates an array of 100,000 elements, each of which is a random floating-point value, and it does it in a fraction of a second. Even more astonishing is the speed with which numpy statistical functions process this array. The following IDLE session demonstrates how quickly you can get stats on this large data set. >>> import numpy as np >>> import numpy.random as ran >>> A = ran.random(100000) >>> np.mean(A) 0.49940282901121 >>> np.sum(A) 49940.282901121005 >>> np.median(A) 0.5005147698475437

>>> np.std(A) 0.2889516828729947 If you try this session yourself, you should experience the response times as instantaneous, even the standard deviation. Most of these stats are straightforward in meaning. Because the probability distribution in this case is a uniform distribution from 0.0 to 1.0, you’d reasonably expect the mean to be close to 0.5, which it is: approximately 0.4994. The sum is exactly 100,000 times that, or about 49,940. The median is not the same as the mean, although you’d expect it to be close to the center of values, which it is: just over 0.50. The standard deviation is what statisticians would predict for a uniform distribution like this: just under 0.29. So roughly 60 percent of the values fall within one standard deviation (plus or minus) of the mean. Using numpy saves you from having to do this calculation yourself, but it’s useful to review how standard deviation is calculated and what it means. Assume A and A2 represent arrays, and i refers to elements: A2 = (i – mean(A)) ^ 2, for all i in A. std(A) = sqrt(mean(A2)) In plainer language, these equations mean the following: Figure out the average of the elements in array A. This is also called the mean. Subtract the mean from each element in A to create a new array full of “deviations.” In this array of deviations, square each member, and call the resulting array A2. Find the average of all elements in A2, take the square root of the result, and voila! You have produced the standard deviation of the array you started with.

Although numpy provides the standard-deviation function for free, it’s useful to see what it would take to produce the result through the standard batch operations. First, getting A2 would be easy enough: Subtracting the mean of A (a scalar) from A itself gives us an array filled with deviations. All these are then squared. A2 = (A – mean(A)) ** 2 Having obtained this new array, we need only get the square root of the mean of the deviations. result = (mean(A2)) ** 0.5 Or, in terms of Python code, it requires the np qualifier to call the mean function: >>> A2 = (A - np.mean(A)) ** 2 >>> result = (np.mean(A2)) ** 0.5 >>> result 0.2889516828729947 >>> np.std(A) 0.2889516828729947 The results, as you can see, are precisely the same in both cases—calculating standard deviation “the hard way” and getting it from np.std—which is good evidence that the numpy routines are following the same algorithm. It’s instructive to run the standard deviation function on an even larger array—say, an array of 1 million random numbers— and the equivalent code in Python, using standard lists. Now comes the interesting part: If you benchmark this technique with a list of 1 million elements in size, against a numpy version with an array containing the same data, the numpy version—getting the standard deviation directly—beats out the non-numpy version by a factor of more than 100 to 1! Here’s the full benchmark code:

import numpy as np import time import numpy.random as ran def get_std1(ls): t1 = time.time() m = sum(ls)/len(ls) ls2 = [(i - m) ** 2 for i in ls] sd = (sum(ls2)/len(ls2)) ** .5 t2 = time.time() print('Python took', t2-t1) def get_std2(A): t1 = time.time() A2 = (A - np.mean(A)) ** 2 result = (np.mean(A2)) ** .5 t2 = time.time() print('Numpy took', t2-t1) def get_std3(A): t1 = time.time() result = np.std(A) t2 = time.time() print('np.std took', t2-t1) A = ran.rand(1000000) get_std1(A) get_std2(A) get_std3(A) Running all three gets the following results, expressed in parts of a second. Remember this is the time taken to get standard deviation for 1 million elements. Python took 0.6885709762573242 Numpy took 0.0189220905303955 np.std took 0.0059509277343750 You can see how enormous the gains in performance are, as we go from Python lists, to numpy arrays, to finally using numpy to get standard deviation directly, with a single function call.

The increase in speed from not using numpy at all, compared to using np.std (the numpy standard deviation function) is well over 100 to 1. Now that’s greased lightning! 12.13 GETTING DATA ON “NUMPY” ROWS AND COLUMNS Section 12.12, “Getting ‘numpy’ Stats (Standard Deviation),” assumed the simplest possible case: getting statistics on a floating-point array of one dimension. However, all those functions accept other arguments. You can look these up in online documentation or in IDLE, through the help command. The most important argument, other than the array itself, is the axis argument, which is used with higher dimensions— that is, dimensions greater than one. Let’s start with another array of random numbers—this time, integers. To produce such an array, you can use the randint function of the numpy.random package. Here’s an example: import numpy as np import numpy.random as ran A = ran.randint(1, 20, (3, 4)) print(A) Here’s some sample output. Your results, of course, will vary. [[ 4 13 11 8] [ 7 14 16 1] [4 1 5 9]] The numpy.random package has its own randint function, just as the random package does. This is another reason that using namespace qualifiers is so important. In this case, using the numpy random package, the function takes begin and

end arguments, as you’d expect, but it also takes an additional argument: a tuple giving the shape of the array. Another thing to note is that with ran.randint, the begin and end arguments include numbers starting with the begin argument, up to but not including the end argument. So this example produces numbers up to 19. Finally, the shape argument—which comes after the begin and end arguments—is (3, 3), which causes the random integers to be placed throughout a 3 × 3 array. Here it is again. Your mileage will vary. [[ 4 13 11 8] [ 7 14 16 1] [4 1 5 9]] As you learned in the previous section, the numpy statistical functions can be used to study this array as one large source of data. If np.mean is applied directly to the array, for example, it gets the mean of all 20 elements. >>> np.mean(A) 7.75 Likewise, we can sum the data or get the standard deviation. Click here to view code image >>> np.sum(A) 93 >>> np.std(A) # standard deviation 4.780603169754489 The fun part comes when we collect statistics along an axis: either the row or the column dimension. Such operations enable the treatment of a numpy array as if it were a spreadsheet, containing totals for each row and column.

However, it’s easy to get the axis confused. Table 12.7 should help clarify. Table 12.7. Use of the “axis” Argument Setting Description axi Create a row collecting data for each column. Size of the s= resulting one-dimensional array is the number of columns. 0 axi Create a column collecting data for each row. Size of the s= resulting one-dimensional array is the number of rows. 1 For even higher dimensional arrays, the axis settings can run higher. The axis settings can even be tuples. Although it may be confusing at first, the way to approach the word “axis” is to think of it like a Cartesian coordinate system, as the name suggests. Look at A again. [[ 4 13 11 8] [ 7 14 16 1] [4 1 5 9]] The argument setting axis=0 refers to the first axis, which means rows (because row-major order is assumed). Therefore, to sum along axis=0 is to sum along the traditional X axis. Summing as it goes, the function sums each column in turn, starting with the lowest-numbered column and moving right. The result is [15 28 32 18]

The argument setting axis=1 refers to the second axis, which means columns. Therefore, to sum along axis=1 is to sum along the traditional Y axis. In this case, the sums start with the lowest-numbered row and move downward. The result is [36 38 19] When summation is done along the X axis, the numpy package collects data on the other dimensions. So, although axis=0 refers to rows, columns are being summed. Figure 12.3 illustrates how this works. Figure 12.3. How “axis=0” and “axis=1” work Let’s take another example; this one is easier to see how its effects work. Let’s start with an array in which each element is equal to its column number. Click here to view code image B = np.fromfunction(lambda r,c: c, (4, 5), dtype=np.int32) Printing this array produces [[0 1 2 3 4] [0 1 2 3 4]

[0 1 2 3 4] [0 1 2 3 4]] Summing along axis 0 (producing column totals) should give us a multiple of 4 each time. Summing along axis 1 (producing row totals) should give us 10 each time. And, in fact, that’s what we get. Click here to view code image >>> np.sum(B, axis = 0) # row, totaling cols. array([ 0, 4, 8, 12, 16]) # col, totaling rows. >>> np.sum(B, axis = 1) array([10, 10, 10, 10]) This is admittedly confusing because axis=0, which should refer to rows, actually sums all the dimensions except rows (in this case, columns). And axis=1 actually sums all the dimensions except columns (in this case, rows). Can we use this data to produce something like a spreadsheet? What we’d like to do is to total all the rows, for example, and then concatenate the results onto the array, using the results as an additional column. Let’s start by, once again, getting both the starting array and the row-bytwo totals. Click here to view code image B = np.fromfunction(lambda r,c:c, (4, 5), dtype=np.int32) B_rows = np.sum(B, axis = 1) Now, can we “glue on” the one-dimensional B_rows to the two-dimensional array, B? Yes, the solution is to use the c_ operator, as follows: B1 = np.c_[B, B_rows]

The array B1 is similar to B, the array we started with, but B1 has an additional column, this one made up of the totals of each row. When printed, it’s displayed as [[ 0 1 2 3 4 10] [ 0 1 2 3 4 10] [ 0 1 2 3 4 10] [ 0 1 2 3 4 10]] This is part of a “spreadsheet” display, with the last column representing sums of rows. With a few more lines of code, we can produce a more complete spreadsheet, in which the bottom row contains sums of the columns. To do this, we get the sums of all the columns of B1. The setting axis=0 moves along the X axis, getting totals as it moves; therefore, it creates a row containing column totals. B_cols = np.sum(B1, axis = 0) The following statement then glues on the row, along the bottom of B1. B2 = np.r_[B1, [B_cols]] Printing B2 prints the following results: [[ 0 1 2 3 4 10] [ 0 1 2 3 4 10] [ 0 1 2 3 4 10] [ 0 1 2 3 4 10] [ 0 4 8 12 16 40]] So there we have it: transformation of an ordinary array into a spreadsheet display that includes totals of all rows and columns along the bottom and the right. The whole procedure can be placed in a function that will operate on any two-dimensional array.

def spreadsheet(A): AC = np.sum(A, axis = 1) A2 = np.c_[A, AC] AR = np.sum(A2, axis = 0) return np.r_[A2, [AR] ] For example, suppose you have the following array: Click here to view code image >>> arr = np.arange(15).reshape(3, 5) >>> print(arr) [[ 0 1 2 3 4] [ 5 6 7 8 9] [10 11 12 13 14]] Here’s what happens if you use the spreadsheet function and print the results: >>> print(spreadsheet(arr)) [[ 0 1 2 3 4 10] [ 5 6 7 8 9 35] [ 10 11 12 13 14 60] [ 15 18 21 24 27 105]] The spreadsheet function can be altered to print summary statistics for other operations, such as mean, median, standard deviation (std), and so on. CHAPTER 12 SUMMARY The numpy package supports manipulation and statistical analysis of large data sets, with abilities that go far beyond those of standard Python arrays. But this chapter, long though it is, has only begun to explore those abilities. One simple test of performance speed is to add up a large set of numbers. In the test of adding up all the numbers from 1 to 1 million, the numpy version of the program beats the ordinary

version by a factor of 10 to 1. But when the benchmark is run on the manipulation of data and not on array creation, the contrast is much greater still. The numpy package provides many ways to create a standard numpy array, called an ndarray, or an “N-dimensional array.” The type is distinguished by the ease with which you can create arrays of more than one dimension. This numpy type has built-in support for statistical analysis, including addition, mean, median, and standard deviation. You can perform these operations on rows, columns, and slices. Much of the power of the numpy ndarray type stems from the ability to take slices of these arrays, either one-dimensional or higher-dimensional, and then perform sophisticated batch operations on them—that is, doing many calculations at once. The slicing ability extends smoothly to any number of dimensions. This numpy type has built-in support for statistical analysis, including addition, mean, median, and standard deviation. You can perform these operations on rows, columns, and slices. In the next chapter, Chapter 13, we’ll explore more advanced capabilities that are built on top of numpy standard types (ndarray), particularly the ability to plot mathematical equations. CHAPTER 12 REVIEW QUESTIONS 1 What are the advantages, if any, of the built-in array package? 2 What are some limitations of the array package? 3 State some major differences between the array package and numpy.

4 Describe the differences between the empty, ones, and zeros functions. 5 What is the role of the callable argument in the fromfunction function used to create new arrays? 6 What happens when a numpy array is combined with a single-value operand (a scalar, such as an int or a floating- point value) through addition, as in the expression A + n? 7 Can combined operation-assign operators (such as += or *=) be used in array-to-scalar operation? What is the effect? 8 Can fixed-length strings be included in a numpy array? What happens when a string of longer length is assigned to such an array? 9 What happens when two numpy arrays are combined through an operation such as addition (+) or multiplication (*)? What requirements must be met to combine two numpy arrays? 10 How do you use a Boolean array as a mask for another array? 11 What are three distinct ways, using both standard Python and its packages, to get the standard deviation of a large set of data? Rank the three of them according to execution speed. 12 What is the dimensionality of an array produced by a Boolean mask? CHAPTER 12 SUGGESTED PROBLEMS 1 Revise the benchmarks used in Section 12.4, “Introduction to ‘numpy’: Sum 1 to 1 Million,” so that they measure the

creation and summation of the data sets separately—that is, report the relative speeds of Python list creation versus numpy array creation, and then report the speeds of adding the numbers together. 2 Use the numpy.random package to generate an array of 1,000 ×1,000 random floating-point numbers. Measure the speed of creating this array. Measure the mean and standard deviation of the numbers in this array. 3 Generate an array of random integers between 0 and 99, inclusive. Then, using Boolean arrays, mask out all integers except those meeting any of the following three conditions: N == 1, N is a multiple of 3, or N > 75. Print the results. 4 Generate a 10 × 10 array full of 1’s. Then zero out the middle 8 × 8 portion, leaving only the outer regions of the array set to 1. That would include the four corners and the edges. Print the results. (Hint: There is more than one way to do this, but slicing is particularly efficient.) 5 Perform a similar operation for a 5 × 5 × 5 “cube,” leaving all visible portions set to 1, while the inner 3 × 3 × 3 cube is zeroed out. Then print the five horizontal planes constituting this cube.

13. Advanced Uses of “numpy” Your introduction to the world of numpy has just begun. One of the most interesting things you can do is plot charts, building on top of numpy data types and functions. You’ll need to download and import the matplotlib package as well as numpy. There are many things you can set, such as color, range, extent, and other factors, but getting started is simple. After covering the plotting capabilities, this chapter examines other advanced uses of the numpy package: Financial applications Linear algebra: dot products and outer products Operations on fixed-length records of all kinds Reading and writing large amounts of data 13.1 ADVANCED MATH OPERATIONS WITH “NUMPY” Before you start plotting curves, you’ll need to know about the mathematical functions operating on the standard numpy ndarray type. Table 13.1 lists the most common of the functions. These are extremely useful in plotting graphs. Most of them take an array as argument and return an array of the same shape. Table 13.1. High-Level “numpy” Math Operations Operation Description

nump Return the cosine of each element in A, as explained in Chapter y.c 11, “The Random and Math Packages.” The input to this os( function is assumed to be in radians and not degrees. That’s A) assumed for the other trig functions as well. nump Return the sine of each element in A. As with cos, the results y.s are returned in an array of matching shape. in( A) nump Return the tangent of each element. y.t an( A) nump Return e raised to the power of each element in A. y.e xp( A) nump Raise each element in X by its corresponding value in Y. These y.p two arrays must have the same shape, or one or both must be owe scalar. r(X , Y) nump Convert degrees to radians; the argument may be an array of the y.r same shape or may be scalar. adi ans (A | x) nump Take the absolute value of each element in A.

y.a bs( A) nump Take the natural logarithm of each element in A. y.l og( A) nump Take the logarithm base 10 of each element in A. y.l og1 0(A ) nump Take the logarithm base 2 of each element in A. y.l og2 (A) nump Take the square root of each element in A. y.s qrt (A) nump Inverse cosine function. y.a rcc os( A) nump Inverse sine function. y.a rcs in( A)

nump Inverse tangent function. y.a rct an( A) nump Hyperbolic cosine. y.h cos (A) nump Hyperbolic sine. y.h sin (A) nump Hyperbolic tangent. y.h tan (A) nump Create a new array by appending the contents of array B onto y.a the end of array A. ppe nd( A, B) nump Get the value of pi. y.p i nump Get the value of e. y.e

These functions have additional arguments. A common argument is the out argument, which names an output array in which to place the results. Such an array must match the source array in size and shape. Here’s an example: import numpy as np A = np.linspace(0, np.pi, 10) B = np.empty(10) np.sin(A, out=B) This last line puts the results into B. The following statement would have the same effect, except that every time you call it, you allocate a new array. In some cases, it’s more efficient to use an existing array to store the results. Click here to view code image import numpy as np A = np.linspace(0, np.pi, 10) B = np.sin(A) To see all the arguments available to be called for each of these functions, use the np.info command. Click here to view code image numpy.info(numpy.function_name) For example, the following commands, given from within IDLE, provide manageable chunks of information on specific numpy functions: