A Python Integer Is More Than Just an Integer

The standard Python implementation is written in C. This means that every Python object is simply a cleverly disguised C structure, which contains not only its value, but other information as well. For example, when we define an integer in Python, such as x = 10000 , x is not just a “raw” integer. It’s actually a pointer to a compound C structure, which contains several values. Looking through the Python 3.4 source code, we find that the integer (long) type definition effectively looks like this (once the C macros are expanded):

 struct    _longobject    {        long    ob_refcnt   ;        PyTypeObject    *   ob_type   ;        size_t    ob_size   ;        long    ob_digit   [   1   ];    };   

A single integer in Python 3.4 actually contains four pieces:

• ob_refcnt , a reference count that helps Python silently handle memory allocation and deallocation
• ob_type , which encodes the type of the variable
• ob_size , which specifies the size of the following data members
• ob_digit , which contains the actual integer value that we expect the Python variable to represent

This means that there is some overhead in storing an integer in Python as compared to an integer in a compiled language like C, as illustrated in .

Figure 2-1. The difference between C and Python integers

Here PyObject_HEAD is the part of the structure containing the reference count, type code, and other pieces mentioned before.

Notice the difference here: a C integer is essentially a label for a position in memory whose bytes encode an integer value. A Python integer is a pointer to a position in memory containing all the Python object information, including the bytes that contain the integer value. This extra information in the Python integer structure is what allows Python to be coded so freely and dynamically. All this additional information in Python types comes at a cost, however, which becomes especially apparent in structures that combine many of these objects.

Fixed-Type Arrays in Python

Python offers several different options for storing data in efficient, fixed-type data buffers. The built-in array module (available since Python 3.3) can be used to create dense arrays of a uniform type:

 In   [   6   ]:    import    array           L    =    list   (   range   (   10   ))           A    =    array   .   array   (   'i'   ,    L   )           A   

Out[6]: array('i', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

Here 'i' is a type code indicating the contents are integers.

Much more useful, however, is the ndarray object of the NumPy package. While Python’s array object provides efficient storage of array-based data, NumPy adds to this efficient operations on that data. We will explore these operations in later sections; here we’ll demonstrate several ways of creating a NumPy array.

We’ll start with the standard NumPy import, under the alias np :

 In   [   7   ]:    import    numpy    as    np   

Creating Arrays from Python Lists

First, we can use np.array to create arrays from Python lists:

 In   [   8   ]:    # integer array:           np   .   array   ([   1   ,    4   ,    2   ,    5   ,    3   ])   

Out[8]: array([1, 4, 2, 5, 3])

Remember that unlike Python lists, NumPy is constrained to arrays that all contain the same type. If types do not match, NumPy will upcast if possible (here, integers are upcast to floating point):

 In   [   9   ]:    np   .   array   ([   3.14   ,    4   ,    2   ,    3   ])   

Out[9]: array([ 3.14,  4.  ,  2.  ,  3.  ])

If we want to explicitly set the data type of the resulting array, we can use the dtype keyword:

 In   [   10   ]:    np   .   array   ([   1   ,    2   ,    3   ,    4   ],    dtype   =   'float32'   )   

Out[10]: array([ 1.,  2.,  3.,  4.], dtype=float32)

Finally, unlike Python lists, NumPy arrays can explicitly be multidimensional; here’s one way of initializing a multidimensional array using a list of lists:

 In   [   11   ]:    # nested lists result in multidimensional arrays            np   .   array   ([   range   (   i   ,    i    +    3   )    for    i    in    [   2   ,    4   ,    6   ]])   

Out[11]: array([[2, 3, 4],                 [4, 5, 6],                 [6, 7, 8]])

The inner lists are treated as rows of the resulting two-dimensional array.

Creating Arrays from Scratch

Especially for larger arrays, it is more efficient to create arrays from scratch using routines built into NumPy. Here are several examples:

 In   [   12   ]:    # Create a length-10 integer array filled with zeros            np   .   zeros   (   10   ,    dtype   =   int   )   

Out[12]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

 In   [   13   ]:    # Create a 3x5 floating-point array filled with 1s            np   .   ones   ((   3   ,    5   ),    dtype   =   float   )   

Out[13]: array([[ 1.,  1.,  1.,  1.,  1.],                 [ 1.,  1.,  1.,  1.,  1.],                 [ 1.,  1.,  1.,  1.,  1.]])

 In   [   14   ]:    # Create a 3x5 array filled with 3.14            np   .   full   ((   3   ,    5   ),    3.14   )   

Out[14]: array([[ 3.14,  3.14,  3.14,  3.14,  3.14],                 [ 3.14,  3.14,  3.14,  3.14,  3.14],                 [ 3.14,  3.14,  3.14,  3.14,  3.14]])

 In   [   15   ]:    # Create an array filled with a linear sequence            # Starting at 0, ending at 20, stepping by 2            # (this is similar to the built-in range() function)            np   .   arange   (   0   ,    20   ,    2   )   

Out[15]: array([ 0,  2,  4,  6,  8, 10, 12, 14, 16, 18])

 In   [   16   ]:    # Create an array of five values evenly spaced between 0 and 1            np   .   linspace   (   0   ,    1   ,    5   )   

Out[16]: array([ 0.  ,  0.25,  0.5 ,  0.75,  1.  ])

 In   [   17   ]:    # Create a 3x3 array of uniformly distributed            # random values between 0 and 1            np   .   random   .   random   ((   3   ,    3   ))   

Out[17]: array([[ 0.99844933,  0.52183819,  0.22421193],                 [ 0.08007488,  0.45429293,  0.20941444],                 [ 0.14360941,  0.96910973,  0.946117  ]])

 In   [   18   ]:    # Create a 3x3 array of normally distributed random values            # with mean 0 and standard deviation 1            np   .   random   .   normal   (   0   ,    1   ,    (   3   ,    3   ))   

Out[18]: array([[ 1.51772646,  0.39614948, -0.10634696],                 [ 0.25671348,  0.00732722,  0.37783601],                 [ 0.68446945,  0.15926039, -0.70744073]])

 In   [   19   ]:    # Create a 3x3 array of random integers in the interval [0, 10)            np   .   random   .   randint   (   0   ,    10   ,    (   3   ,    3   ))   

Out[19]: array([[2, 3, 4],                 [5, 7, 8],                 [0, 5, 0]])

 In   [   20   ]:    # Create a 3x3 identity matrix            np   .   eye   (   3   )   

Out[20]: array([[ 1.,  0.,  0.],                 [ 0.,  1.,  0.],                 [ 0.,  0.,  1.]])

 In   [   21   ]:    # Create an uninitialized array of three integers            # The values will be whatever happens to already exist at that            # memory location            np   .   empty   (   3   )   

Out[21]: array([ 1.,  1.,  1.])

NumPy Standard Data Types

NumPy arrays contain values of a single type, so it is important to have detailed knowledge of those types and their limitations. Because NumPy is built in C, the types will be familiar to users of C, Fortran, and other related languages.

The standard NumPy data types are listed in . Note that when constructing an array, you can specify them using a string:

 np   .   zeros   (   10   ,    dtype   =   'int16'   )   

Or using the associated NumPy object:

 np   .   zeros   (   10   ,    dtype   =   np   .   int16   )   

More advanced type specification is possible, such as specifying big or little endian numbers; for more information, refer to the . NumPy also supports compound data types, which will be covered in .

NumPy Array Attributes

First let’s discuss some useful array attributes. We’ll start by defining three random arrays: a one-dimensional, two-dimensional, and three-dimensional array. We’ll use NumPy’s random number generator, which we will seed with a set value in order to ensure that the same random arrays are generated each time this code is run:

 In   [   1   ]:    import    numpy    as    np           np   .   random   .   seed   (   0   )     # seed for reproducibility            x1    =    np   .   random   .   randint   (   10   ,    size   =   6   )     # One-dimensional array           x2    =    np   .   random   .   randint   (   10   ,    size   =   (   3   ,    4   ))     # Two-dimensional array           x3    =    np   .   random   .   randint   (   10   ,    size   =   (   3   ,    4   ,    5   ))     # Three-dimensional array   

Each array has attributes ndim (the number of dimensions), shape (the size of each dimension), and size (the total size of the array):

 In   [   2   ]:    print   (   "x3 ndim: "   ,    x3   .   ndim   )           print   (   "x3 shape:"   ,    x3   .   shape   )           print   (   "x3 size: "   ,    x3   .   size   )   

x3 ndim:  3 x3 shape: (3, 4, 5) x3 size:  60

Another useful attribute is the dtype , the data type of the array (which we discussed previously in ):

 In   [   3   ]:    print   (   "dtype:"   ,    x3   .   dtype   )   

dtype: int64

Other attributes include itemsize , which lists the size (in bytes) of each array element, and nbytes , which lists the total size (in bytes) of the array:

 In   [   4   ]:    print   (   "itemsize:"   ,    x3   .   itemsize   ,    "bytes"   )           print   (   "nbytes:"   ,    x3   .   nbytes   ,    "bytes"   )   

itemsize: 8 bytes nbytes: 480 bytes

In general, we expect that nbytes is equal to itemsize times size .

Array Indexing: Accessing Single Elements

If you are familiar with Python’s standard list indexing, indexing in NumPy will feel quite familiar. In a one-dimensional array, you can access the i ^th value (counting from zero) by specifying the desired index in square brackets, just as with Python lists:

 In   [   5   ]:    x1   

Out[5]: array([5, 0, 3, 3, 7, 9])

 In   [   6   ]:    x1   [   0   ]   

Out[6]: 5

 In   [   7   ]:    x1   [   4   ]   

Out[7]: 7

To index from the end of the array, you can use negative indices:

 In   [   8   ]:    x1   [   -   1   ]   

Out[8]: 9

 In   [   9   ]:    x1   [   -   2   ]   

Out[9]: 7

In a multidimensional array, you access items using a comma-separated tuple of indices :

 In   [   10   ]:    x2   

Out[10]: array([[3, 5, 2, 4],                 [7, 6, 8, 8],                 [1, 6, 7, 7]])

 In   [   11   ]:    x2   [   0   ,    0   ]   

Out[11]: 3

 In   [   12   ]:    x2   [   2   ,    0   ]   

Out[12]: 1

 In   [   13   ]:    x2   [   2   ,    -   1   ]   

Out[13]: 7

You can also modify values using any of the above index notation:

 In   [   14   ]:    x2   [   0   ,    0   ]    =    12            x2   

Out[14]: array([[12,  5,  2,  4],                 [ 7,  6,  8,  8],                 [ 1,  6,  7,  7]])

Keep in mind that, unlike Python lists, NumPy arrays have a fixed type. This means, for example, that if you attempt to insert a floating-point value to an integer array, the value will be silently truncated. Don’t be caught unaware by this behavior!

 In   [   15   ]:    x1   [   0   ]    =    3.14159     # this will be truncated!            x1   

Out[15]: array([3, 0, 3, 3, 7, 9])

Array Slicing: Accessing Subarrays

Just as we can use square brackets to access individual array elements, we can also use them to access subarrays with the slice notation, marked by the colon (: ) character. The NumPy slicing syntax follows that of the standard Python list; to access a slice of an array x , use this:

 x   [   start   :   stop   :   step   ]   

If any of these are unspecified, they default to the values start=0 , stop=size of dimension , step=1 . We’ll take a look at accessing subarrays in one dimension and in multiple dimensions.

 In   [   16   ]:    x    =    np   .   arange   (   10   )            x   

Out[16]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

 In   [   17   ]:    x   [:   5   ]     # first five elements   

Out[17]: array([0, 1, 2, 3, 4])

 In   [   18   ]:    x   [   5   :]     # elements after index 5   

Out[18]: array([5, 6, 7, 8, 9])

 In   [   19   ]:    x   [   4   :   7   ]     # middle subarray   

Out[19]: array([4, 5, 6])

 In   [   20   ]:    x   [::   2   ]     # every other element   

Out[20]: array([0, 2, 4, 6, 8])

 In   [   21   ]:    x   [   1   ::   2   ]     # every other element, starting at index 1   

Out[21]: array([1, 3, 5, 7, 9])

 In   [   22   ]:    x   [::   -   1   ]     # all elements, reversed   

Out[22]: array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])

 In   [   23   ]:    x   [   5   ::   -   2   ]     # reversed every other from index 5   

Out[23]: array([5, 3, 1])

 In   [   24   ]:    x2   

Out[24]: array([[12,  5,  2,  4],                 [ 7,  6,  8,  8],                 [ 1,  6,  7,  7]])

 In   [   25   ]:    x2   [:   2   ,    :   3   ]     # two rows, three columns   

Out[25]: array([[12,  5,  2],                 [ 7,  6,  8]])

 In   [   26   ]:    x2   [:   3   ,    ::   2   ]     # all rows, every other column   

Out[26]: array([[12,  2],                 [ 7,  8],                 [ 1,  7]])

 In   [   27   ]:    x2   [::   -   1   ,    ::   -   1   ]   

Out[27]: array([[ 7,  7,  6,  1],                 [ 8,  8,  6,  7],                 [ 4,  2,  5, 12]])

 In   [   28   ]:    print   (   x2   [:,    0   ])     # first column of x2   

[12  7  1]

 In   [   29   ]:    print   (   x2   [   0   ,    :])     # first row of x2   

[12  5  2  4]

 In   [   30   ]:    print   (   x2   [   0   ])     # equivalent to x2[0, :]   

[12  5  2  4]

 In   [   31   ]:    print   (   x2   )   

[[12  5  2  4]  [ 7  6  8  8]  [ 1  6  7  7]]

 In   [   32   ]:    x2_sub    =    x2   [:   2   ,    :   2   ]            print   (   x2_sub   )   

[[12  5]  [ 7  6]]

 In   [   33   ]:    x2_sub   [   0   ,    0   ]    =    99            print   (   x2_sub   )   

[[99  5]  [ 7  6]]

 In   [   34   ]:    print   (   x2   )   

[[99  5  2  4]  [ 7  6  8  8]  [ 1  6  7  7]]

 In   [   35   ]:    x2_sub_copy    =    x2   [:   2   ,    :   2   ]   .   copy   ()            print   (   x2_sub_copy   )   

[[99  5]  [ 7  6]]

 In   [   36   ]:    x2_sub_copy   [   0   ,    0   ]    =    42            print   (   x2_sub_copy   )   

[[42  5]  [ 7  6]]

 In   [   37   ]:    print   (   x2   )   

[[99  5  2  4]  [ 7  6  8  8]  [ 1  6  7  7]]

Reshaping of Arrays

Another useful type of operation is reshaping of arrays. The most flexible way of doing this is with the reshape() method. For example, if you want to put the numbers 1 through 9 in a 3×3 grid, you can do the following:

 In   [   38   ]:    grid    =    np   .   arange   (   1   ,    10   )   .   reshape   ((   3   ,    3   ))            print   (   grid   )   

[[1 2 3]  [4 5 6]  [7 8 9]]

Note that for this to work, the size of the initial array must match the size of the reshaped array. Where possible, the reshape method will use a no-copy view of the initial array, but with noncontiguous memory buffers this is not always the case.

Another common reshaping pattern is the conversion of a one-dimensional array into a two-dimensional row or column matrix. You can do this with the reshape method, or more easily by making use of the newaxis keyword within a slice operation:

 In   [   39   ]:    x    =    np   .   array   ([   1   ,    2   ,    3   ])             # row vector via reshape            x   .   reshape   ((   1   ,    3   ))   

Out[39]: array([[1, 2, 3]])

 In   [   40   ]:    # row vector via newaxis            x   [   np   .   newaxis   ,    :]   

Out[40]: array([[1, 2, 3]])

 In   [   41   ]:    # column vector via reshape            x   .   reshape   ((   3   ,    1   ))   

Out[41]: array([[1],                 [2],                 [3]])

 In   [   42   ]:    # column vector via newaxis            x   [:,    np   .   newaxis   ]   

Out[42]: array([[1],                 [2],                 [3]])

We will see this type of transformation often throughout the remainder of the book.

Array Concatenation and Splitting

All of the preceding routines worked on single arrays. It’s also possible to combine multiple arrays into one, and to conversely split a single array into multiple arrays. We’ll take a look at those operations here.

 In   [   43   ]:    x    =    np   .   array   ([   1   ,    2   ,    3   ])            y    =    np   .   array   ([   3   ,    2   ,    1   ])            np   .   concatenate   ([   x   ,    y   ])   

Out[43]: array([1, 2, 3, 3, 2, 1])

 In   [   44   ]:    z    =    [   99   ,    99   ,    99   ]            print   (   np   .   concatenate   ([   x   ,    y   ,    z   ]))   

[ 1  2  3  3  2  1 99 99 99]

 In   [   45   ]:    grid    =    np   .   array   ([[   1   ,    2   ,    3   ],                             [   4   ,    5   ,    6   ]])   

 In   [   46   ]:    # concatenate along the first axis            np   .   concatenate   ([   grid   ,    grid   ])   

Out[46]: array([[1, 2, 3],                 [4, 5, 6],                 [1, 2, 3],                 [4, 5, 6]])

 In   [   47   ]:    # concatenate along the second axis (zero-indexed)            np   .   concatenate   ([   grid   ,    grid   ],    axis   =   1   )   

Out[47]: array([[1, 2, 3, 1, 2, 3],                 [4, 5, 6, 4, 5, 6]])

 In   [   48   ]:    x    =    np   .   array   ([   1   ,    2   ,    3   ])            grid    =    np   .   array   ([[   9   ,    8   ,    7   ],                             [   6   ,    5   ,    4   ]])             # vertically stack the arrays            np   .   vstack   ([   x   ,    grid   ])   

Out[48]: array([[1, 2, 3],                 [9, 8, 7],                 [6, 5, 4]])

 In   [   49   ]:    # horizontally stack the arrays            y    =    np   .   array   ([[   99   ],                          [   99   ]])            np   .   hstack   ([   grid   ,    y   ])   

Out[49]: array([[ 9,  8,  7, 99],                 [ 6,  5,  4, 99]])

 In   [   50   ]:    x    =    [   1   ,    2   ,    3   ,    99   ,    99   ,    3   ,    2   ,    1   ]            x1   ,    x2   ,    x3    =    np   .   split   (   x   ,    [   3   ,    5   ])            print   (   x1   ,    x2   ,    x3   )   

[1 2 3] [99 99] [3 2 1]

 In   [   51   ]:    grid    =    np   .   arange   (   16   )   .   reshape   ((   4   ,    4   ))            grid   

Out[51]: array([[ 0,  1,  2,  3],                 [ 4,  5,  6,  7],                 [ 8,  9, 10, 11],                 [12, 13, 14, 15]])

 In   [   52   ]:    upper   ,    lower    =    np   .   vsplit   (   grid   ,    [   2   ])            print   (   upper   )            print   (   lower   )   

[[0 1 2 3]  [4 5 6 7]] [[ 8  9 10 11]  [12 13 14 15]]

 In   [   53   ]:    left   ,    right    =    np   .   hsplit   (   grid   ,    [   2   ])            print   (   left   )            print   (   right   )   

[[ 0  1]  [ 4  5]  [ 8  9]  [12 13]] [[ 2  3]  [ 6  7]  [10 11]  [14 15]]

The Slowness of Loops

Python’s default implementation (known as CPython) does some operations very slowly. This is in part due to the dynamic, interpreted nature of the language: the fact that types are flexible, so that sequences of operations cannot be compiled down to efficient machine code as in languages like C and Fortran. Recently there have been various attempts to address this weakness: well-known examples are the , a just-in-time compiled implementation of Python; the , which converts Python code to compilable C code; and the , which converts snippets of Python code to fast LLVM bytecode. Each of these has its strengths and weaknesses, but it is safe to say that none of the three approaches has yet surpassed the reach and popularity of the standard CPython engine.

The relative sluggishness of Python generally manifests itself in situations where many small operations are being repeated — for instance, looping over arrays to operate on each element. For example, imagine we have an array of values and we’d like to compute the reciprocal of each. A straightforward approach might look like this:

 In   [   1   ]:    import    numpy    as    np           np   .   random   .   seed   (   0   )            def    compute_reciprocals   (   values   ):               output    =    np   .   empty   (   len   (   values   ))               for    i    in    range   (   len   (   values   )):                   output   [   i   ]    =    1.0    /    values   [   i   ]               return    output            values    =    np   .   random   .   randint   (   1   ,    10   ,    size   =   5   )           compute_reciprocals   (   values   )   

Out[1]: array([ 0.16666667,  1.        ,  0.25      ,  0.25      ,  0.125     ])

This implementation probably feels fairly natural to someone from, say, a C or Java background. But if we measure the execution time of this code for a large input, we see that this operation is very slow, perhaps surprisingly so! We’ll benchmark this with IPython’s %timeit magic (discussed in ):

 In   [   2   ]:    big_array    =    np   .   random   .   randint   (   1   ,    100   ,    size   =   1000000   )           %   timeit    compute_reciprocals   (   big_array   )   

1 loop, best of 3: 2.91 s per loop

It takes several seconds to compute these million operations and to store the result! When even cell phones have processing speeds measured in Giga-FLOPS (i.e., billions of numerical operations per second), this seems almost absurdly slow. It turns out that the bottleneck here is not the operations themselves, but the type-checking and function dispatches that CPython must do at each cycle of the loop. Each time the reciprocal is computed, Python first examines the object’s type and does a dynamic lookup of the correct function to use for that type. If we were working in compiled code instead, this type specification would be known before the code executes and the result could be computed much more efficiently.

Introducing UFuncs

For many types of operations, NumPy provides a convenient interface into just this kind of statically typed, compiled routine. This is known as a vectorized operation. You can accomplish this by simply performing an operation on the array, which will then be applied to each element. This vectorized approach is designed to push the loop into the compiled layer that underlies NumPy, leading to much faster execution.

Compare the results of the following two:

 In   [   3   ]:    print   (   compute_reciprocals   (   values   ))           print   (   1.0    /    values   )   

[ 0.16666667  1.          0.25        0.25        0.125     ] [ 0.16666667  1.          0.25        0.25        0.125     ]

Looking at the execution time for our big array, we see that it completes orders of magnitude faster than the Python loop:

 In   [   4   ]:    %   timeit    (   1.0    /    big_array   )   

100 loops, best of 3: 4.6 ms per loop

Vectorized operations in NumPy are implemented via ufuncs , whose main purpose is to quickly execute repeated operations on values in NumPy arrays. Ufuncs are extremely flexible — before we saw an operation between a scalar and an array, but we can also operate between two arrays:

 In   [   5   ]:    np   .   arange   (   5   )    /    np   .   arange   (   1   ,    6   )   

Out[5]: array([ 0.        ,  0.5       ,  0.66666667,  0.75      ,  0.8       ])

And ufunc operations are not limited to one-dimensional arrays — they can act on multidimensional arrays as well:

 In   [   6   ]:    x    =    np   .   arange   (   9   )   .   reshape   ((   3   ,    3   ))           2    **    x   

Out[6]: array([[  1,   2,   4],                [  8,  16,  32],                [ 64, 128, 256]])

Computations using vectorization through ufuncs are nearly always more efficient than their counterpart implemented through Python loops, especially as the arrays grow in size. Any time you see such a loop in a Python script, you should consider whether it can be replaced with a vectorized expression.

Exploring NumPy’s UFuncs

Ufuncs exist in two flavors: unary ufuncs , which operate on a single input, and binary ufuncs , which operate on two inputs. We’ll see examples of both these types of functions here.

 In   [   7   ]:    x    =    np   .   arange   (   4   )           print   (   "x     ="   ,    x   )           print   (   "x + 5 ="   ,    x    +    5   )           print   (   "x - 5 ="   ,    x    -    5   )           print   (   "x * 2 ="   ,    x    *    2   )           print   (   "x / 2 ="   ,    x    /    2   )           print   (   "x // 2 ="   ,    x    //    2   )     # floor division   

x     = [0 1 2 3] x + 5 = [5 6 7 8] x - 5 = [-5 -4 -3 -2] x * 2 = [0 2 4 6] x / 2 = [ 0.   0.5  1.   1.5] x // 2 = [0 0 1 1]

 In   [   8   ]:    print   (   "-x     = "   ,    -   x   )           print   (   "x ** 2 = "   ,    x    **    2   )           print   (   "x % 2  = "   ,    x    %    2   )   

-x     =  [ 0 -1 -2 -3] x ** 2 =  [0 1 4 9] x % 2  =  [0 1 0 1]

 In   [   9   ]:    -   (   0.5   *   x    +    1   )    **    2   

Out[9]: array([-1.  , -2.25, -4.  , -6.25])

 In   [   10   ]:    np   .   add   (   x   ,    2   )   

Out[10]: array([2, 3, 4, 5])

 In   [   11   ]:    x    =    np   .   array   ([   -   2   ,    -   1   ,    0   ,    1   ,    2   ])            abs   (   x   )   

Out[11]: array([2, 1, 0, 1, 2])

 In   [   12   ]:    np   .   absolute   (   x   )   

Out[12]: array([2, 1, 0, 1, 2])

 In   [   13   ]:    np   .   abs   (   x   )   

Out[13]: array([2, 1, 0, 1, 2])

 In   [   14   ]:    x    =    np   .   array   ([   3    -    4j   ,    4    -    3j   ,    2    +    0j   ,    0    +    1j   ])            np   .   abs   (   x   )   

Out[14]: array([ 5.,  5.,  2.,  1.])

 In   [   15   ]:    theta    =    np   .   linspace   (   0   ,    np   .   pi   ,    3   )   

 In   [   16   ]:    print   (   "theta      = "   ,    theta   )            print   (   "sin(theta) = "   ,    np   .   sin   (   theta   ))            print   (   "cos(theta) = "   ,    np   .   cos   (   theta   ))            print   (   "tan(theta) = "   ,    np   .   tan   (   theta   ))   

theta      =  [ 0.          1.57079633  3.14159265] sin(theta) =  [  0.00000000e+00   1.00000000e+00   1.22464680e-16] cos(theta) =  [  1.00000000e+00   6.12323400e-17  -1.00000000e+00] tan(theta) =  [  0.00000000e+00   1.63312394e+16  -1.22464680e-16]

 In   [   17   ]:    x    =    [   -   1   ,    0   ,    1   ]            print   (   "x         = "   ,    x   )            print   (   "arcsin(x) = "   ,    np   .   arcsin   (   x   ))            print   (   "arccos(x) = "   ,    np   .   arccos   (   x   ))            print   (   "arctan(x) = "   ,    np   .   arctan   (   x   ))   

x         =  [-1, 0, 1] arcsin(x) =  [-1.57079633  0.          1.57079633] arccos(x) =  [ 3.14159265  1.57079633  0.        ] arctan(x) =  [-0.78539816  0.          0.78539816]

 In   [   18   ]:    x    =    [   1   ,    2   ,    3   ]            print   (   "x     ="   ,    x   )            print   (   "e^x   ="   ,    np   .   exp   (   x   ))            print   (   "2^x   ="   ,    np   .   exp2   (   x   ))            print   (   "3^x   ="   ,    np   .   power   (   3   ,    x   ))   

x     = [1, 2, 3] e^x   = [  2.71828183   7.3890561   20.08553692] 2^x   = [ 2.  4.  8.] 3^x   = [ 3  9 27]

 In   [   19   ]:    x    =    [   1   ,    2   ,    4   ,    10   ]            print   (   "x        ="   ,    x   )            print   (   "ln(x)    ="   ,    np   .   log   (   x   ))            print   (   "log2(x)  ="   ,    np   .   log2   (   x   ))            print   (   "log10(x) ="   ,    np   .   log10   (   x   ))   

x        = [1, 2, 4, 10] ln(x)    = [ 0.          0.69314718  1.38629436  2.30258509] log2(x)  = [ 0.          1.          2.          3.32192809] log10(x) = [ 0.          0.30103     0.60205999  1.        ]

 In   [   20   ]:    x    =    [   0   ,    0.001   ,    0.01   ,    0.1   ]            print   (   "exp(x) - 1 ="   ,    np   .   expm1   (   x   ))            print   (   "log(1 + x) ="   ,    np   .   log1p   (   x   ))   

exp(x) - 1 = [ 0.          0.0010005   0.01005017  0.10517092] log(1 + x) = [ 0.          0.0009995   0.00995033  0.09531018]

 In   [   21   ]:    from    scipy    import    special   

 In   [   22   ]:    # Gamma functions (generalized factorials) and related functions            x    =    [   1   ,    5   ,    10   ]            print   (   "gamma(x)     ="   ,    special   .   gamma   (   x   ))            print   (   "ln|gamma(x)| ="   ,    special   .   gammaln   (   x   ))            print   (   "beta(x, 2)   ="   ,    special   .   beta   (   x   ,    2   ))   

gamma(x)     = [  1.00000000e+00   2.40000000e+01   3.62880000e+05] ln|gamma(x)| = [  0.           3.17805383  12.80182748] beta(x, 2)   = [ 0.5         0.03333333  0.00909091]

 In   [   23   ]:    # Error function (integral of Gaussian)            # its complement, and its inverse            x    =    np   .   array   ([   0   ,    0.3   ,    0.7   ,    1.0   ])            print   (   "erf(x)  ="   ,    special   .   erf   (   x   ))            print   (   "erfc(x) ="   ,    special   .   erfc   (   x   ))            print   (   "erfinv(x) ="   ,    special   .   erfinv   (   x   ))   

erf(x)  = [ 0.          0.32862676  0.67780119  0.84270079] erfc(x) = [ 1.          0.67137324  0.32219881  0.15729921] erfinv(x) = [ 0.          0.27246271  0.73286908         inf]

Advanced Ufunc Features

Many NumPy users make use of ufuncs without ever learning their full set of features. We’ll outline a few specialized features of ufuncs here.

 In   [   24   ]:    x    =    np   .   arange   (   5   )            y    =    np   .   empty   (   5   )            np   .   multiply   (   x   ,    10   ,    out   =   y   )            print   (   y   )   

[  0.  10.  20.  30.  40.]

 In   [   25   ]:    y    =    np   .   zeros   (   10   )            np   .   power   (   2   ,    x   ,    out   =   y   [::   2   ])            print   (   y   )   

[  1.   0.   2.   0.   4.   0.   8.   0.  16.   0.]

 In   [   26   ]:    x    =    np   .   arange   (   1   ,    6   )            np   .   add   .   reduce   (   x   )   

Out[26]: 15

 In   [   27   ]:    np   .   multiply   .   reduce   (   x   )   

Out[27]: 120

 In   [   28   ]:    np   .   add   .   accumulate   (   x   )   

Out[28]: array([ 1,  3,  6, 10, 15])

 In   [   29   ]:    np   .   multiply   .   accumulate   (   x   )   

Out[29]: array([  1,   2,   6,  24, 120])

 In   [   30   ]:    x    =    np   .   arange   (   1   ,    6   )            np   .   multiply   .   outer   (   x   ,    x   )   

Out[30]: array([[ 1,  2,  3,  4,  5],                 [ 2,  4,  6,  8, 10],                 [ 3,  6,  9, 12, 15],                 [ 4,  8, 12, 16, 20],                 [ 5, 10, 15, 20, 25]])

Ufuncs: Learning More

More information on universal functions (including the full list of available functions) can be found on the and documentation websites.

Recall that you can also access information directly from within IPython by importing the packages and using IPython’s tab-completion and help (? ) functionality, as described in .

Summing the Values in an Array

As a quick example, consider computing the sum of all values in an array. Python itself can do this using the built-in sum function:

 In   [   1   ]:    import    numpy    as    np   

 In   [   2   ]:    L    =    np   .   random   .   random   (   100   )           sum   (   L   )   

Out[2]: 55.61209116604941

The syntax is quite similar to that of NumPy’s sum function, and the result is the same in the simplest case:

 In   [   3   ]:    np   .   sum   (   L   )   

Out[3]: 55.612091166049424

However, because it executes the operation in compiled code, NumPy’s version of the operation is computed much more quickly:

 In   [   4   ]:    big_array    =    np   .   random   .   rand   (   1000000   )           %   timeit    sum   (   big_array   )           %   timeit    np   .   sum   (   big_array   )   

10 loops, best of 3: 104 ms per loop 1000 loops, best of 3: 442 µs per loop

Be careful, though: the sum function and the np.sum function are not identical, which can sometimes lead to confusion! In particular, their optional arguments have different meanings, and np.sum is aware of multiple array dimensions, as we will see in the following section.

Minimum and Maximum

Similarly, Python has built-in min and max functions, used to find the minimum value and maximum value of any given array:

 In   [   5   ]:    min   (   big_array   ),    max   (   big_array   )   

Out[5]: (1.1717128136634614e-06, 0.9999976784968716)

NumPy’s corresponding functions have similar syntax, and again operate much more quickly:

 In   [   6   ]:    np   .   min   (   big_array   ),    np   .   max   (   big_array   )   

Out[6]: (1.1717128136634614e-06, 0.9999976784968716)

 In   [   7   ]:    %   timeit    min   (   big_array   )           %   timeit    np   .   min   (   big_array   )   

10 loops, best of 3: 82.3 ms per loop 1000 loops, best of 3: 497 µs per loop

For min , max , sum , and several other NumPy aggregates, a shorter syntax is to use methods of the array object itself:

 In   [   8   ]:    print   (   big_array   .   min   (),    big_array   .   max   (),    big_array   .   sum   ())   

1.17171281366e-06 0.999997678497 499911.628197

Whenever possible, make sure that you are using the NumPy version of these aggregates when operating on NumPy arrays!

 In   [   9   ]:    M    =    np   .   random   .   random   ((   3   ,    4   ))           print   (   M   )   

[[ 0.8967576   0.03783739  0.75952519  0.06682827]  [ 0.8354065   0.99196818  0.19544769  0.43447084]  [ 0.66859307  0.15038721  0.37911423  0.6687194 ]]

 In   [   10   ]:    M   .   sum   ()   

Out[10]: 6.0850555667307118

 In   [   11   ]:    M   .   min   (   axis   =   0   )   

Out[11]: array([ 0.66859307,  0.03783739,  0.19544769,  0.06682827])

 In   [   12   ]:    M   .   max   (   axis   =   1   )   

Out[12]: array([ 0.8967576 ,  0.99196818,  0.6687194 ])

Example: What Is the Average Height of US Presidents?

Aggregates available in NumPy can be extremely useful for summarizing a set of values. As a simple example, let’s consider the heights of all US presidents. This data is available in the file president_heights.csv , which is a simple comma-separated list of labels and values:

 In   [   13   ]:    !   head    -   4    data   /   president_heights   .   csv   

order,name,height(cm) 1,George Washington,189 2,John Adams,170 3,Thomas Jefferson,189

We’ll use the Pandas package, which we’ll explore more fully in , to read the file and extract this information (note that the heights are measured in centimeters):

 In   [   14   ]:    import    pandas    as    pd            data    =    pd   .   read_csv   (   'data/president_heights.csv'   )            heights    =    np   .   array   (   data   [   'height(cm)'   ])            print   (   heights   )   

[189 170 189 163 183 171 185 168 173 183 173 173 175 178 183 193 178 173  174 183 183 168 170 178 182 180 183 178 182 188 175 179 183 193 182 183  177 185 188 188 182 185]

Now that we have this data array, we can compute a variety of summary statistics:

 In   [   15   ]:    print   (   "Mean height:       "   ,    heights   .   mean   ())            print   (   "Standard deviation:"   ,    heights   .   std   ())            print   (   "Minimum height:    "   ,    heights   .   min   ())            print   (   "Maximum height:    "   ,    heights   .   max   ())   

Mean height:        179.738095238 Standard deviation: 6.93184344275 Minimum height:     163 Maximum height:     193

Note that in each case, the aggregation operation reduced the entire array to a single summarizing value, which gives us information about the distribution of values. We may also wish to compute quantiles:

 In   [   16   ]:    print   (   "25th percentile:   "   ,    np   .   percentile   (   heights   ,    25   ))            print   (   "Median:            "   ,    np   .   median   (   heights   ))            print   (   "75th percentile:   "   ,    np   .   percentile   (   heights   ,    75   ))   

25th percentile:    174.25 Median:             182.0 75th percentile:    183.0

We see that the median height of US presidents is 182 cm, or just shy of six feet.

Of course, sometimes it’s more useful to see a visual representation of this data, which we can accomplish using tools in Matplotlib (we’ll discuss Matplotlib more fully in ). For example, this code generates the chart shown in :

 In   [   17   ]:    %   matplotlib    inline            import    matplotlib.pyplot    as    plt            import    seaborn   ;    seaborn   .   set   ()     # set plot style   

 In   [   18   ]:    plt   .   hist   (   heights   )            plt   .   title   (   'Height Distribution of US Presidents'   )            plt   .   xlabel   (   'height (cm)'   )            plt   .   ylabel   (   'number'   );   

Figure 2-3. Histogram of presidential heights

These aggregates are some of the fundamental pieces of exploratory data analysis that we’ll explore in more depth in later chapters of the book.

Introducing Broadcasting

Recall that for arrays of the same size, binary operations are performed on an element-by-element basis:

 In   [   1   ]:    import    numpy    as    np   

 In   [   2   ]:    a    =    np   .   array   ([   0   ,    1   ,    2   ])           b    =    np   .   array   ([   5   ,    5   ,    5   ])           a    +    b   

Out[2]: array([5, 6, 7])

Broadcasting allows these types of binary operations to be performed on arrays of different sizes — for example, we can just as easily add a scalar (think of it as a zero-dimensional array) to an array:

 In   [   3   ]:    a    +    5   

Out[3]: array([5, 6, 7])

We can think of this as an operation that stretches or duplicates the value 5 into the array [5, 5, 5] , and adds the results. The advantage of NumPy’s broadcasting is that this duplication of values does not actually take place, but it is a useful mental model as we think about broadcasting.

We can similarly extend this to arrays of higher dimension. Observe the result when we add a one-dimensional array to a two-dimensional array:

 In   [   4   ]:    M    =    np   .   ones   ((   3   ,    3   ))           M   

Out[4]: array([[ 1.,  1.,  1.],                [ 1.,  1.,  1.],                [ 1.,  1.,  1.]])

 In   [   5   ]:    M    +    a   

Out[5]: array([[ 1.,  2.,  3.],                [ 1.,  2.,  3.],                [ 1.,  2.,  3.]])

Here the one-dimensional array a is stretched, or broadcast, across the second dimension in order to match the shape of M .

While these examples are relatively easy to understand, more complicated cases can involve broadcasting of both arrays. Consider the following example:

 In   [   6   ]:    a    =    np   .   arange   (   3   )           b    =    np   .   arange   (   3   )[:,    np   .   newaxis   ]            print   (   a   )           print   (   b   )   

[0 1 2] [[0]  [1]  [2]]

 In   [   7   ]:    a    +    b   

Out[7]: array([[0, 1, 2],                [1, 2, 3],                [2, 3, 4]])

Just as before we stretched or broadcasted one value to match the shape of the other, here we’ve stretched both a and b to match a common shape, and the result is a two-dimensional array! The geometry of these examples is visualized in .

Figure 2-4. Visualization of NumPy broadcasting

The light boxes represent the broadcasted values: again, this extra memory is not actually allocated in the course of the operation, but it can be useful conceptually to imagine that it is.

Rules of Broadcasting

Broadcasting in NumPy follows a strict set of rules to determine the interaction between the two arrays:

• Rule 1: If the two arrays differ in their number of dimensions, the shape of the one with fewer dimensions is padded with ones on its leading (left) side.
• Rule 2: If the shape of the two arrays does not match in any dimension, the array with shape equal to 1 in that dimension is stretched to match the other shape.
• Rule 3: If in any dimension the sizes disagree and neither is equal to 1, an error is raised.

To make these rules clear, let’s consider a few examples in detail.

 In   [   8   ]:    M    =    np   .   ones   ((   2   ,    3   ))           a    =    np   .   arange   (   3   )   

 In   [   9   ]:    M    +    a   

Out[9]: array([[ 1.,  2.,  3.],                [ 1.,  2.,  3.]])

 In   [   10   ]:    a    =    np   .   arange   (   3   )   .   reshape   ((   3   ,    1   ))            b    =    np   .   arange   (   3   )   

 In   [   11   ]:    a    +    b   

Out[11]: array([[0, 1, 2],                 [1, 2, 3],                 [2, 3, 4]])

 In   [   12   ]:    M    =    np   .   ones   ((   3   ,    2   ))            a    =    np   .   arange   (   3   )   

 In   [   13   ]:    M    +    a   

 ---------------------------------------------------------------------------  ValueError                                Traceback (most recent call last)  <ipython-input-13-9e16e9f98da6> in <module>() ----> 1 M + a   ValueError: operands could not be broadcast together with shapes (3,2) (3,)

 In   [   14   ]:    a   [:,    np   .   newaxis   ]   .   shape   

Out[14]: (3, 1)

 In   [   15   ]:    M    +    a   [:,    np   .   newaxis   ]   

Out[15]: array([[ 1.,  1.],                 [ 2.,  2.],                 [ 3.,  3.]])

 In   [   16   ]:    np   .   logaddexp   (   M   ,    a   [:,    np   .   newaxis   ])   

Out[16]: array([[ 1.31326169,  1.31326169],                 [ 1.69314718,  1.69314718],                 [ 2.31326169,  2.31326169]])

Broadcasting in Practice

Broadcasting operations form the core of many examples we’ll see throughout this book. We’ll now take a look at a couple simple examples of where they can be useful.

 In   [   17   ]:    X    =    np   .   random   .   random   ((   10   ,    3   ))   

 In   [   18   ]:    Xmean    =    X   .   mean   (   0   )            Xmean   

Out[18]: array([ 0.53514715,  0.66567217,  0.44385899])

 In   [   19   ]:    X_centered    =    X    -    Xmean   

 In   [   20   ]:    X_centered   .   mean   (   0   )   

Out[20]: array([  2.22044605e-17,  -7.77156117e-17,  -1.66533454e-17])

Plotting a two-dimensional function

One place that broadcasting is very useful is in displaying images based on two-dimensional functions. If we want to define a function z = f (x , y ), broadcasting can be used to compute the function across the grid:

 In   [   21   ]:    # x and y have 50 steps from 0 to 5            x    =    np   .   linspace   (   0   ,    5   ,    50   )            y    =    np   .   linspace   (   0   ,    5   ,    50   )[:,    np   .   newaxis   ]             z    =    np   .   sin   (   x   )    **    10    +    np   .   cos   (   10    +    y    *    x   )    *    np   .   cos   (   x   )   

We’ll use Matplotlib to plot this two-dimensional array (these tools will be discussed in full in ):

 In   [   22   ]:    %   matplotlib    inline           import    matplotlib.pyplot    as    plt   

 In   [   23   ]:    plt   .   imshow   (   z   ,    origin   =   'lower'   ,    extent   =   [   0   ,    5   ,    0   ,    5   ],                       cmap   =   'viridis'   )            plt   .   colorbar   ();   

The result, shown in , is a compelling visualization of the two-dimensional function.

Figure 2-5. Visualization of a 2D array

Example: Counting Rainy Days

Imagine you have a series of data that represents the amount of precipitation each day for a year in a given city. For example, here we’ll load the daily rainfall statistics for the city of Seattle in 2014, using Pandas (which is covered in more detail in ):

 In   [   1   ]:    import    numpy    as    np           import    pandas    as    pd            # use Pandas to extract rainfall inches as a NumPy array           rainfall    =    pd   .   read_csv   (   'data/Seattle2014.csv'   )[   'PRCP'   ]   .   values           inches    =    rainfall    /    254     # 1/10mm -> inches           inches   .   shape   

Out[1]: (365,)

The array contains 365 values, giving daily rainfall in inches from January 1 to December 31, 2014.

As a first quick visualization, let’s look at the histogram of rainy days shown in , which was generated using Matplotlib (we will explore this tool more fully in ):

 In   [   2   ]:    %   matplotlib    inline           import    matplotlib.pyplot    as    plt           import    seaborn   ;    seaborn   .   set   ()     # set plot styles   

 In   [   3   ]:    plt   .   hist   (   inches   ,    40   );   

Figure 2-6. Histogram of 2014 rainfall in Seattle

This histogram gives us a general idea of what the data looks like: despite its reputation, the vast majority of days in Seattle saw near zero measured rainfall in 2014. But this doesn’t do a good job of conveying some information we’d like to see: for example, how many rainy days were there in the year? What is the average precipitation on those rainy days? How many days were there with more than half an inch of rain?

Comparison Operators as ufuncs

In we introduced ufuncs, and focused in particular on arithmetic operators. We saw that using + , - , * , / , and others on arrays leads to element-wise operations. NumPy also implements comparison operators such as < (less than) and > (greater than) as element-wise ufuncs. The result of these comparison operators is always an array with a Boolean data type. All six of the standard comparison operations are available:

 In   [   4   ]:    x    =    np   .   array   ([   1   ,    2   ,    3   ,    4   ,    5   ])   

 In   [   5   ]:    x    <    3     # less than   

Out[5]: array([ True,  True, False, False, False], dtype=bool)

 In   [   6   ]:    x    >    3     # greater than   

Out[6]: array([False, False, False,  True,  True], dtype=bool)

 In   [   7   ]:    x    <=    3     # less than or equal   

Out[7]: array([ True,  True,  True, False, False], dtype=bool)

 In   [   8   ]:    x    >=    3     # greater than or equal   

Out[8]: array([False, False,  True,  True,  True], dtype=bool)

 In   [   9   ]:    x    !=    3     # not equal   

Out[9]: array([ True,  True, False,  True,  True], dtype=bool)

 In   [   10   ]:    x    ==    3     # equal   

Out[10]: array([False, False,  True, False, False], dtype=bool)

It is also possible to do an element-by-element comparison of two arrays, and to include compound expressions:

 In   [   11   ]:    (   2    *    x   )    ==    (   x    **    2   )   

Out[11]: array([False,  True, False, False, False], dtype=bool)

As in the case of arithmetic operators, the comparison operators are implemented as ufuncs in NumPy; for example, when you write x < 3 , internally NumPy uses np.less(x, 3) . A summary of the comparison operators and their equivalent ufunc is shown here:

Just as in the case of arithmetic ufuncs, these will work on arrays of any size and shape. Here is a two-dimensional example:

 In   [   12   ]:    rng    =    np   .   random   .   RandomState   (   0   )            x    =    rng   .   randint   (   10   ,    size   =   (   3   ,    4   ))            x   

Out[12]: array([[5, 0, 3, 3],                 [7, 9, 3, 5],                 [2, 4, 7, 6]])

 In   [   13   ]:    x    <    6   

Out[13]: array([[ True,  True,  True,  True],                 [False, False,  True,  True],                 [ True,  True, False, False]], dtype=bool)

In each case, the result is a Boolean array, and NumPy provides a number of straightforward patterns for working with these Boolean results.

Working with Boolean Arrays

Given a Boolean array, there are a host of useful operations you can do. We’ll work with x , the two-dimensional array we created earlier:

 In   [   14   ]:    print   (   x   )   

[[5 0 3 3]  [7 9 3 5]  [2 4 7 6]]

 In   [   15   ]:    # how many values less than 6?            np   .   count_nonzero   (   x    <    6   )   

Out[15]: 8

 In   [   16   ]:    np   .   sum   (   x    <    6   )   

Out[16]: 8

 In   [   17   ]:    # how many values less than 6 in each row?            np   .   sum   (   x    <    6   ,    axis   =   1   )   

Out[17]: array([4, 2, 2])

 In   [   18   ]:    # are there any values greater than 8?            np   .   any   (   x    >    8   )   

Out[18]: True

 In   [   19   ]:    # are there any values less than zero?            np   .   any   (   x    <    0   )   

Out[19]: False

 In   [   20   ]:    # are all values less than 10?            np   .   all   (   x    <    10   )   

Out[20]: True

 In   [   21   ]:    # are all values equal to 6?            np   .   all   (   x    ==    6   )   

Out[21]: False

 In   [   22   ]:    # are all values in each row less than 8?            np   .   all   (   x    <    8   ,    axis   =   1   )   

Out[22]: array([ True, False,  True], dtype=bool)

 In   [   23   ]:    np   .   sum   ((   inches    >    0.5   )    &    (   inches    <    1   ))   

Out[23]: 29

 inches    >    (   0.5    &    inches   )    <    1   

 In   [   24   ]:    np   .   sum   (   ~   (    (   inches    <=    0.5   )    |    (   inches    >=    1   )    ))   

Out[24]: 29

 In   [   25   ]:    print   (   "Number days without rain:      "   ,    np   .   sum   (   inches    ==    0   ))            print   (   "Number days with rain:         "   ,    np   .   sum   (   inches    !=    0   ))            print   (   "Days with more than 0.5 inches:"   ,    np   .   sum   (   inches    >    0.5   ))            print   (   "Rainy days with < 0.1 inches  :"   ,    np   .   sum   ((   inches    >    0   )    &                                                            (   inches    <    0.2   )))   

Number days without rain:       215 Number days with rain:          150 Days with more than 0.5 inches: 37 Rainy days with < 0.1 inches  : 75

Boolean Arrays as Masks

In the preceding section, we looked at aggregates computed directly on Boolean arrays. A more powerful pattern is to use Boolean arrays as masks, to select particular subsets of the data themselves. Returning to our x array from before, suppose we want an array of all values in the array that are less than, say, 5:

 In   [   26   ]:    x   

Out[26]: array([[5, 0, 3, 3],                 [7, 9, 3, 5],                 [2, 4, 7, 6]])

We can obtain a Boolean array for this condition easily, as we’ve already seen:

 In   [   27   ]:    x    <    5   

Out[27]: array([[False,  True,  True,  True],                 [False, False,  True, False],                 [ True,  True, False, False]], dtype=bool)

Now to select these values from the array, we can simply index on this Boolean array; this is known as a masking operation:

 In   [   28   ]:    x   [   x    <    5   ]   

Out[28]: array([0, 3, 3, 3, 2, 4])

What is returned is a one-dimensional array filled with all the values that meet this condition; in other words, all the values in positions at which the mask array is True .

We are then free to operate on these values as we wish. For example, we can compute some relevant statistics on our Seattle rain data:

 In   [   29   ]:    # construct a mask of all rainy days    rainy    =    (   inches    >    0   )     # construct a mask of all summer days (June 21st is the 172nd day)    summer    =    (   np   .   arange   (   365   )    -    172    <    90   )    &    (   np   .   arange   (   365   )    -    172    >    0   )     print   (   "Median precip on rainy days in 2014 (inches):   "   ,          np   .   median   (   inches   [   rainy   ]))    print   (   "Median precip on summer days in 2014 (inches):  "   ,          np   .   median   (   inches   [   summer   ]))    print   (   "Maximum precip on summer days in 2014 (inches): "   ,          np   .   max   (   inches   [   summer   ]))    print   (   "Median precip on non-summer rainy days (inches):"   ,          np   .   median   (   inches   [   rainy    &    ~   summer   ]))   

Median precip on rainy days in 2014 (inches):    0.194881889764 Median precip on summer days in 2014 (inches):   0.0 Maximum precip on summer days in 2014 (inches):  0.850393700787 Median precip on non-summer rainy days (inches): 0.200787401575

By combining Boolean operations, masking operations, and aggregates, we can very quickly answer these sorts of questions for our dataset.

The difference is this: and and or gauge the truth or falsehood of entire object , while & and | refer to bits within each object .

When you use and or or , it’s equivalent to asking Python to treat the object as a single Boolean entity. In Python, all nonzero integers will evaluate as True . Thus:

 In   [   30   ]:    bool   (   42   ),    bool   (   0   )   

Out[30]: (True, False)

 In   [   31   ]:    bool   (   42    and    0   )   

Out[31]: False

 In   [   32   ]:    bool   (   42    or    0   )   

Out[32]: True

When you use & and | on integers, the expression operates on the bits of the element, applying the and or the or to the individual bits making up the number:

 In   [   33   ]:    bin   (   42   )   

Out[33]: '0b101010'

 In   [   34   ]:    bin   (   59   )   

Out[34]: '0b111011'

 In   [   35   ]:    bin   (   42    &    59   )   

Out[35]: '0b101010'

 In   [   36   ]:    bin   (   42    |    59   )   

Out[36]: '0b111011'

Notice that the corresponding bits of the binary representation are compared in order to yield the result.

When you have an array of Boolean values in NumPy, this can be thought of as a string of bits where 1 = True and 0 = False , and the result of & and | operates in a similar manner as before:

 In   [   37   ]:    A    =    np   .   array   ([   1   ,    0   ,    1   ,    0   ,    1   ,    0   ],    dtype   =   bool   )            B    =    np   .   array   ([   1   ,    1   ,    1   ,    0   ,    1   ,    1   ],    dtype   =   bool   )            A    |    B   

Out[37]: array([ True,  True,  True, False,  True,  True], dtype=bool)

Using or on these arrays will try to evaluate the truth or falsehood of the entire array object, which is not a well-defined value:

 In   [   38   ]:    A    or    B   

 ---------------------------------------------------------------------------  ValueError                                Traceback (most recent call last)  <ipython-input-38-5d8e4f2e21c0> in <module>() ----> 1 A or B   ValueError: The truth value of an array with more than one element is...

Similarly, when doing a Boolean expression on a given array, you should use | or & rather than or or and :

 In   [   39   ]:    x    =    np   .   arange   (   10   )            (   x    >    4   )    &    (   x    <    8   )   

Out[39]: array([False, False, ...,  True,  True, False, False], dtype=bool)

Trying to evaluate the truth or falsehood of the entire array will give the same ValueError we saw previously:

 In   [   40   ]:    (   x    >    4   )    and    (   x    <    8   )   

 ---------------------------------------------------------------------------  ValueError                                Traceback (most recent call last)  <ipython-input-40-3d24f1ffd63d> in <module>() ----> 1 (x > 4) and (x < 8)   ValueError: The truth value of an array with more than one element is...

So remember this: and and or perform a single Boolean evaluation on an entire object, while & and | perform multiple Boolean evaluations on the content (the individual bits or bytes) of an object. For Boolean NumPy arrays, the latter is nearly always the desired operation .