NumPy#

We always import NumPy using the standard np abbreviation.

import numpy as np

Regular arrays#

To learn some different methods for creating NumPy arrays, we will make the following arrays. (The next one is written as a list, but our primary goal is to make the NumPy array version.)

Example 1#

  • [0, 0, 0, …, 0] (length 10)

Here we use np.zeros. Many Matlab functions have an equivalent version in NumPy, and zeros is one such example.

np.zeros(10)
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

Just as an aside, there is also a very elegant way to make a list (not a NumPy array) of 10 zeros using what is called list comprehension. (We will introduce list comprehension more systematically later in the course.) The underscore _ in the following is just a way of naming a variable that we will never refer to again; it could be replaced by any variable name and the result would be the same.

Here is the list comprehension way to make a list of 10 zeros.

[0 for _ in range(10)]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

One other way of making this list is based on the following, which uses + to concatenate two lists.

[0]+[0]
[0, 0]

So if we want to add [0] to itself 10 times, we can use *.

[0]*10
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Example 2#

  • A five row, seven column matrix containing all 4s.

Here we’re trying to make a 5x7 NumPy array of all zeros (then we will add 4 to it). The following is a common mistake I make. The error message is not too helpful the first time you see it.

np.zeros(5,7)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [6], in <cell line: 1>()
----> 1 np.zeros(5,7)

TypeError: Cannot interpret '7' as a data type

Notice in the documentation that the first argument for shape is supposed to be an integer or a tuple of integers. The second positional argument is called dtype, and NumPy is trying to interpret our second input 7 as a dtype. That was our mistake; we should have given a single input argument as a tuple, not two separate input arguments.

help(np.zeros)
Help on built-in function zeros in module numpy:

zeros(...)
    zeros(shape, dtype=float, order='C', *, like=None)
    
    Return a new array of given shape and type, filled with zeros.
    
    Parameters
    ----------
    shape : int or tuple of ints
        Shape of the new array, e.g., ``(2, 3)`` or ``2``.
    dtype : data-type, optional
        The desired data-type for the array, e.g., `numpy.int8`.  Default is
        `numpy.float64`.
    order : {'C', 'F'}, optional, default: 'C'
        Whether to store multi-dimensional data in row-major
        (C-style) or column-major (Fortran-style) order in
        memory.
    like : array_like, optional
        Reference object to allow the creation of arrays which are not
        NumPy arrays. If an array-like passed in as ``like`` supports
        the ``__array_function__`` protocol, the result will be defined
        by it. In this case, it ensures the creation of an array object
        compatible with that passed in via this argument.
    
        .. versionadded:: 1.20.0
    
    Returns
    -------
    out : ndarray
        Array of zeros with the given shape, dtype, and order.
    
    See Also
    --------
    zeros_like : Return an array of zeros with shape and type of input.
    empty : Return a new uninitialized array.
    ones : Return a new array setting values to one.
    full : Return a new array of given shape filled with value.
    
    Examples
    --------
    >>> np.zeros(5)
    array([ 0.,  0.,  0.,  0.,  0.])
    
    >>> np.zeros((5,), dtype=int)
    array([0, 0, 0, 0, 0])
    
    >>> np.zeros((2, 1))
    array([[ 0.],
           [ 0.]])
    
    >>> s = (2,2)
    >>> np.zeros(s)
    array([[ 0.,  0.],
           [ 0.,  0.]])
    
    >>> np.zeros((2,), dtype=[('x', 'i4'), ('y', 'i4')]) # custom dtype
    array([(0, 0), (0, 0)],
          dtype=[('x', '<i4'), ('y', '<i4')])

Here is the correct way to make a 5x7 NumPy array of all zeros. Notice the extra set of parentheses around (5,7). The inner parentheses make a tuple, and the outer parentheses go together with the zeros function.

np.zeros((5,7))
array([[0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0.]])

Once we have the array of all zeros, it’s easy to get the array of all 4s. We just add 4, and NumPy automatically “broadcasts” the 4 to all of the entries in the array.

np.zeros((5,7))+4
array([[4., 4., 4., 4., 4., 4., 4.],
       [4., 4., 4., 4., 4., 4., 4.],
       [4., 4., 4., 4., 4., 4., 4.],
       [4., 4., 4., 4., 4., 4., 4.],
       [4., 4., 4., 4., 4., 4., 4.]])

Here is a similar strategy, using the ones function instead of the zeros function.

np.ones((5,7))
array([[1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1.]])
np.ones((5,7))*4
array([[4., 4., 4., 4., 4., 4., 4.],
       [4., 4., 4., 4., 4., 4., 4.],
       [4., 4., 4., 4., 4., 4., 4.],
       [4., 4., 4., 4., 4., 4., 4.],
       [4., 4., 4., 4., 4., 4., 4.]])

If we would rather have integers on the inside (notice how the decimal point goes away), we can specify int as the dtype.

np.ones((5,7), int)*4
array([[4, 4, 4, 4, 4, 4, 4],
       [4, 4, 4, 4, 4, 4, 4],
       [4, 4, 4, 4, 4, 4, 4],
       [4, 4, 4, 4, 4, 4, 4],
       [4, 4, 4, 4, 4, 4, 4]])

Example 3#

  • [0, …, 100] (length 5, evenly distributed).

In Matlab, this would be made using linspace, and NumPy has its own version of linspace.

np.linspace(0,100,5)
array([  0.,  25.,  50.,  75., 100.])

Example 4#

  • The 3x3 matrix \(\begin{pmatrix} 0 & 0 & 0 \\ 1 & 1 & 1 \\ 2 & 2 & 2 \end{pmatrix}\)

Here is a first way to make this matrix. It’s again using broadcasting. For example, when i is 2, we broadcast 2 to the second row. (The second row is specified by arr[2].) We will later see ways to make this sort of matrix avoiding the for loop.

arr = np.zeros((3,3))
for i in range(3):
    arr[i] = i
arr
array([[0., 0., 0.],
       [1., 1., 1.],
       [2., 2., 2.]])

Example 5#

  • The 2x5 matrix \( \begin{pmatrix} 2 & 5 & 8 & 11 & 14 \\ 17 & 20 & 23 & 26 & 29 \end{pmatrix} \)

We start with a nested for loop approach to make this matrix. We will print out the current version of the array each step through the for loops. Notice especially how arr[i,j] refers to the entry in the i-th row and the j-th column. The following is almost correct, but it doesn’t get the bottom row right.

arr = np.zeros((2,5))
for i in range(2):
    for j in range(5):
        arr[i,j] = 2+3*j
        print(arr)
[[2. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]
[[2. 5. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]
[[2. 5. 8. 0. 0.]
 [0. 0. 0. 0. 0.]]
[[ 2.  5.  8. 11.  0.]
 [ 0.  0.  0.  0.  0.]]
[[ 2.  5.  8. 11. 14.]
 [ 0.  0.  0.  0.  0.]]
[[ 2.  5.  8. 11. 14.]
 [ 2.  0.  0.  0.  0.]]
[[ 2.  5.  8. 11. 14.]
 [ 2.  5.  0.  0.  0.]]
[[ 2.  5.  8. 11. 14.]
 [ 2.  5.  8.  0.  0.]]
[[ 2.  5.  8. 11. 14.]
 [ 2.  5.  8. 11.  0.]]
[[ 2.  5.  8. 11. 14.]
 [ 2.  5.  8. 11. 14.]]

Here is a final version of this approach. We make the code a little more stable by replacing some of the numbers with variables, for example, we replace the 3 above with the variable step.

cols = 5
step = 3
arr = np.zeros((2,cols))
for i in range(2):
    for j in range(cols):
        arr[i,j] = 2+step*j+cols*step*i
print(arr)
[[ 2.  5.  8. 11. 14.]
 [17. 20. 23. 26. 29.]]

Here is an easier way to make the same NumPy array. Here is a first step: np.arange is the NumPy analogue of range from Python (or the colon function from Matlab). (We’ll see later that np.arange has the advantage over range that it can use non-integer step sizes.)

np.arange(2,29,3)
array([ 2,  5,  8, 11, 14, 17, 20, 23, 26])

Remember that the right endpoint is usually not included in Python. Since we want our array to end at 29, we should list some number bigger than 29 (and less than or equal to 32) as our right endpoint.

np.arange(2,30,3)
array([ 2,  5,  8, 11, 14, 17, 20, 23, 26, 29])

Now we will use the reshape method to turn this length-10 one-dimensional NumPy array into a 2x5 two-dimensional NumPy array.

np.arange(2,30,3).reshape((2,5))
array([[ 2,  5,  8, 11, 14],
       [17, 20, 23, 26, 29]])

Random numbers#

  • Make a length 10 NumPy array of random integers between 0 (inclusive) and 40 (exclusive).

You will often see code like the following in code examples, but this is actually the “old” way of getting random numbers in NumPy.

np.random.randint(0, 40, size=10)
array([30, 35, 37,  3,  0, 18, 12, 17, 38, 16])

Notice the portion of the documentation:

New code should use the integers method of a default_rng() instance instead; please see the :ref:random-quick-start

For the rest of the course, we will use this default_rng method.

help(np.random.randint)
Help on built-in function randint:

randint(...) method of numpy.random.mtrand.RandomState instance
    randint(low, high=None, size=None, dtype=int)
    
    Return random integers from `low` (inclusive) to `high` (exclusive).
    
    Return random integers from the "discrete uniform" distribution of
    the specified dtype in the "half-open" interval [`low`, `high`). If
    `high` is None (the default), then results are from [0, `low`).
    
    .. note::
        New code should use the ``integers`` method of a ``default_rng()``
        instance instead; please see the :ref:`random-quick-start`.
    
    Parameters
    ----------
    low : int or array-like of ints
        Lowest (signed) integers to be drawn from the distribution (unless
        ``high=None``, in which case this parameter is one above the
        *highest* such integer).
    high : int or array-like of ints, optional
        If provided, one above the largest (signed) integer to be drawn
        from the distribution (see above for behavior if ``high=None``).
        If array-like, must contain integer values
    size : int or tuple of ints, optional
        Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
        ``m * n * k`` samples are drawn.  Default is None, in which case a
        single value is returned.
    dtype : dtype, optional
        Desired dtype of the result. Byteorder must be native.
        The default value is int.
    
        .. versionadded:: 1.11.0
    
    Returns
    -------
    out : int or ndarray of ints
        `size`-shaped array of random integers from the appropriate
        distribution, or a single such random int if `size` not provided.
    
    See Also
    --------
    random_integers : similar to `randint`, only for the closed
        interval [`low`, `high`], and 1 is the lowest value if `high` is
        omitted.
    random.Generator.integers: which should be used for new code.
    
    Examples
    --------
    >>> np.random.randint(2, size=10)
    array([1, 0, 0, 0, 1, 1, 0, 0, 1, 0]) # random
    >>> np.random.randint(1, size=10)
    array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
    
    Generate a 2 x 4 array of ints between 0 and 4, inclusive:
    
    >>> np.random.randint(5, size=(2, 4))
    array([[4, 0, 2, 1], # random
           [3, 2, 2, 0]])
    
    Generate a 1 x 3 array with 3 different upper bounds
    
    >>> np.random.randint(1, [3, 5, 10])
    array([2, 2, 9]) # random
    
    Generate a 1 by 3 array with 3 different lower bounds
    
    >>> np.random.randint([1, 5, 7], 10)
    array([9, 8, 7]) # random
    
    Generate a 2 by 4 array using broadcasting with dtype of uint8
    
    >>> np.random.randint([1, 3, 5, 7], [[10], [20]], dtype=np.uint8)
    array([[ 8,  6,  9,  7], # random
           [ 1, 16,  9, 12]], dtype=uint8)

The default_rng approach to random numbers is a good example of object-oriented programming, because we are making a new object, which we call rng here, and we are using methods on that object to produce random numbers.

rng = np.random.default_rng()
type(rng)
numpy.random._generator.Generator

The note above said to use the integers method, which is what we do here. This is the modern way (as of Summer 2022) to make random numbers in NumPy.

arr = rng.integers(0, 40, size=10)
arr
array([10, 36, 39, 25, 23, 30,  5, 20, 25, 12])

For the rest of this section, we will see a variety of examples of different sorts of methods we can use from our random number generator rng object.

  • Choose 6 of those numbers (with replacement) and put them into a NumPy array.

rng.choice(arr, size=6)
array([23, 30, 23, 30, 30, 39])
  • Make a three row, five column NumPy array of random real numbers between -1 and 4.

The following is what I would have expected, but the random method does not allow us to specify a low and a high value for getting random floats. Instead random will always get us random numbers between 0 and 1, and if we want our random numbers in a different range, then we need to rescale the outputs.

rng.random(-1, 4, size=(3,5))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [11], in <cell line: 1>()
----> 1 rng.random(-1, 4, size=(3,5))

File _generator.pyx:238, in numpy.random._generator.Generator.random()

TypeError: random() got multiple values for keyword argument 'size'

Notice that the outputs are between 0 and 1.

help(rng.random)
Help on built-in function random:

random(...) method of numpy.random._generator.Generator instance
    random(size=None, dtype=np.float64, out=None)
    
    Return random floats in the half-open interval [0.0, 1.0).
    
    Results are from the "continuous uniform" distribution over the
    stated interval.  To sample :math:`Unif[a, b), b > a` multiply
    the output of `random` by `(b-a)` and add `a`::
    
      (b - a) * random() + a
    
    Parameters
    ----------
    size : int or tuple of ints, optional
        Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
        ``m * n * k`` samples are drawn.  Default is None, in which case a
        single value is returned.
    dtype : dtype, optional
        Desired dtype of the result, only `float64` and `float32` are supported.
        Byteorder must be native. The default value is np.float64.
    out : ndarray, optional
        Alternative output array in which to place the result. If size is not None,
        it must have the same shape as the provided size and must match the type of
        the output values.
    
    Returns
    -------
    out : float or ndarray of floats
        Array of random floats of shape `size` (unless ``size=None``, in which
        case a single float is returned).
    
    Examples
    --------
    >>> rng = np.random.default_rng()
    >>> rng.random()
    0.47108547995356098 # random
    >>> type(rng.random())
    <class 'float'>
    >>> rng.random((5,))
    array([ 0.30220482,  0.86820401,  0.1654503 ,  0.11659149,  0.54323428]) # random
    
    Three-by-two array of random numbers from [-5, 0):
    
    >>> 5 * rng.random((3, 2)) - 5
    array([[-3.99149989, -0.52338984], # random
           [-2.99091858, -0.79479508],
           [-1.23204345, -1.75224494]])
arr = rng.random((3,5))
arr
array([[0.14373169, 0.04151085, 0.85160258, 0.5625753 , 0.77434184],
       [0.09323801, 0.63691547, 0.83783101, 0.60707909, 0.69066845],
       [0.65647018, 0.45466071, 0.77161217, 0.53020513, 0.71808796]])

The random numbers in arr span from 0 to 1 (a distance of 1 apart). We want to go between -1 and 4, so we eventually want a width of 5. That is why we will multiply by 5.

arr*5
array([[0.71865844, 0.20755425, 4.25801288, 2.81287651, 3.87170919],
       [0.46619004, 3.18457733, 4.18915507, 3.03539543, 3.45334223],
       [3.28235092, 2.27330357, 3.85806084, 2.65102563, 3.59043979]])

The above array goes between 0 and 5, but we want it to go between -1 and 4, so we subract 1.

arr*5-1
array([[-0.28134156, -0.79244575,  3.25801288,  1.81287651,  2.87170919],
       [-0.53380996,  2.18457733,  3.18915507,  2.03539543,  2.45334223],
       [ 2.28235092,  1.27330357,  2.85806084,  1.65102563,  2.59043979]])

Here is a way to make this array all at once (of course the specific random numbers will be different).

5*rng.random(size=(3,5)) - 1
array([[-0.01973026,  1.61155096,  2.9866675 ,  1.71899874, -0.01933733],
       [ 1.36928355,  2.32786313, -0.00448403, -0.82777187,  0.89141526],
       [ 3.13028253, -0.38060258,  3.05444416,  0.42173001,  1.3847107 ]])
  • Make a length 10 NumPy array of random numbers following a normal distribution with mean 2 and standard deviation 0.1.

As just a last example in this section, here is how we can make normally distributed random numbers.

help(rng.normal)
Help on built-in function normal:

normal(...) method of numpy.random._generator.Generator instance
    normal(loc=0.0, scale=1.0, size=None)
    
    Draw random samples from a normal (Gaussian) distribution.
    
    The probability density function of the normal distribution, first
    derived by De Moivre and 200 years later by both Gauss and Laplace
    independently [2]_, is often called the bell curve because of
    its characteristic shape (see the example below).
    
    The normal distributions occurs often in nature.  For example, it
    describes the commonly occurring distribution of samples influenced
    by a large number of tiny, random disturbances, each with its own
    unique distribution [2]_.
    
    Parameters
    ----------
    loc : float or array_like of floats
        Mean ("centre") of the distribution.
    scale : float or array_like of floats
        Standard deviation (spread or "width") of the distribution. Must be
        non-negative.
    size : int or tuple of ints, optional
        Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
        ``m * n * k`` samples are drawn.  If size is ``None`` (default),
        a single value is returned if ``loc`` and ``scale`` are both scalars.
        Otherwise, ``np.broadcast(loc, scale).size`` samples are drawn.
    
    Returns
    -------
    out : ndarray or scalar
        Drawn samples from the parameterized normal distribution.
    
    See Also
    --------
    scipy.stats.norm : probability density function, distribution or
        cumulative density function, etc.
    
    Notes
    -----
    The probability density for the Gaussian distribution is
    
    .. math:: p(x) = \frac{1}{\sqrt{ 2 \pi \sigma^2 }}
                     e^{ - \frac{ (x - \mu)^2 } {2 \sigma^2} },
    
    where :math:`\mu` is the mean and :math:`\sigma` the standard
    deviation. The square of the standard deviation, :math:`\sigma^2`,
    is called the variance.
    
    The function has its peak at the mean, and its "spread" increases with
    the standard deviation (the function reaches 0.607 times its maximum at
    :math:`x + \sigma` and :math:`x - \sigma` [2]_).  This implies that
    :meth:`normal` is more likely to return samples lying close to the
    mean, rather than those far away.
    
    References
    ----------
    .. [1] Wikipedia, "Normal distribution",
           https://en.wikipedia.org/wiki/Normal_distribution
    .. [2] P. R. Peebles Jr., "Central Limit Theorem" in "Probability,
           Random Variables and Random Signal Principles", 4th ed., 2001,
           pp. 51, 51, 125.
    
    Examples
    --------
    Draw samples from the distribution:
    
    >>> mu, sigma = 0, 0.1 # mean and standard deviation
    >>> s = np.random.default_rng().normal(mu, sigma, 1000)
    
    Verify the mean and the variance:
    
    >>> abs(mu - np.mean(s))
    0.0  # may vary
    
    >>> abs(sigma - np.std(s, ddof=1))
    0.0  # may vary
    
    Display the histogram of the samples, along with
    the probability density function:
    
    >>> import matplotlib.pyplot as plt
    >>> count, bins, ignored = plt.hist(s, 30, density=True)
    >>> plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) *
    ...                np.exp( - (bins - mu)**2 / (2 * sigma**2) ),
    ...          linewidth=2, color='r')
    >>> plt.show()
    
    Two-by-four array of samples from N(3, 6.25):
    
    >>> np.random.default_rng().normal(3, 2.5, size=(2, 4))
    array([[-4.49401501,  4.00950034, -1.81814867,  7.29718677],   # random
           [ 0.39924804,  4.68456316,  4.99394529,  4.84057254]])  # random

Notice that the result of the following does seem to be clustered around 2. That makes sense because we had mean 2 with a relatively low standard deviation (0.1).

rng.normal(2, 0.1, size=10)
array([1.97074003, 2.10991763, 1.98318881, 2.00287628, 2.10707454,
       2.02365598, 1.89343959, 1.97254061, 1.97359416, 1.96238154])

The array is displayed almost like it’s two-dimensional, but notice that there is only one set of brackets. It really is one-dimensional, as we can verify by checking its shape attribute.

rng.normal(2, 0.1, size=10).shape
(10,)

Changing rows and columns#

Here we’re going to start with (a slight variant of) the code from an earlier example.

arr = np.zeros((4,4))
for i in range(4):
    arr[i] = i
arr
array([[0., 0., 0., 0.],
       [1., 1., 1., 1.],
       [2., 2., 2., 2.],
       [3., 3., 3., 3.]])
arr = np.zeros((4,4))
arr
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])

If we use indexing of the form arr[i] where arr is a two-dimensional array, this arr[i] will represent the i-th row of arr.

arr[2]
array([0., 0., 0., 0.])

Here we set the second row (where as always we start counting from 0) to be all 13s.

arr[2] = 13
arr
array([[ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [13., 13., 13., 13.],
       [ 0.,  0.,  0.,  0.]])

Accessing a column is slightly harder. You should read arr[:,2] as indicating “every row, the second column”.

arr[:,2] = -4
arr
array([[ 0.,  0., -4.,  0.],
       [ 0.,  0., -4.,  0.],
       [13., 13., -4., 13.],
       [ 0.,  0., -4.,  0.]])

The same sort of indexing with a colon could be used with rows. Here arr[1,:] is the same as arr[1]. The : in this case can be read as, “every column”. In the following, we go to row 1 and set all of its entries to be 1.

arr[1,:] = 1
arr
array([[ 0.,  0., -4.,  0.],
       [ 1.,  1.,  1.,  1.],
       [13., 13., -4., 13.],
       [ 0.,  0., -4.,  0.]])

In the following, we are storing the second column with the variable name v. This column is a one-dimensional NumPy array.

v = arr[:, 2]

Notice how v is displayed horizontally. Looking at v, there is no evidence that this corresponds to a column in our array.

v
array([-4.,  1., -4., -4.])

Here is how a length-1 tuple is displayed in Python, in this case, the length-1 tuple that has only the number 4 in it is written (4,).

v.shape
(4,)

The following code is trying to assign the 3rd row to have the values [2,10]. In this case, we get an error message. The error message says that we cannot broadcast from shape (2,) to shape (4,). In this case, [2,10] has shape (2,) and arr[3] has shape (4,). Broadcasting is a difficult and important notion in NumPy.

arr[3] = [2,10]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [16], in <cell line: 1>()
----> 1 arr[3] = [2,10]

ValueError: could not broadcast input array from shape (2,) into shape (4,)

Here is an example of two different shapes where we can broadcast from one to the other. In this case, [1,3,4,7] has shape (4,) and arr[:] stands for all the entries in the entire array arr, which has shape (4,4). These two shapes are compatible with respect to broadcasting.

Two shapes in NumPy are compatible with respect to broadcasting if, starting from the right-most dimension (the last dimension), the two dimensions are either equal or one of them is equal to 1.

arr[:] = [1,3,4,7]
arr
array([[1., 3., 4., 7.],
       [1., 3., 4., 7.],
       [1., 3., 4., 7.],
       [1., 3., 4., 7.]])
w = np.array([1,3,4,7])
arr.shape
(4, 4)
w.shape
(4,)

Here is a for-loop version of creating an array arr.

arr = np.zeros((3,3))
for i in range(3):
    arr[i] = i
arr
array([[0., 0., 0.],
       [1., 1., 1.],
       [2., 2., 2.]])

Here is a way that uses broadcasting followed by taking the matrix transpose.

arr = np.zeros((3,3))
arr[:] = np.arange(3)
arr.T
array([[0., 0., 0.],
       [1., 1., 1.],
       [2., 2., 2.]])

Here is what arr looked like before we took the matrix transpose. We broadcast from np.arange(3) to the 3x3 array arr. If we want to actually change arr, we need to reassign arr after taking the transpose, using arr = arr.T.

arr
array([[0., 1., 2.],
       [0., 1., 2.],
       [0., 1., 2.]])
arr = np.zeros((3,3))
arr[:] = np.arange(3)
arr = arr.T
arr
array([[0., 0., 0.],
       [1., 1., 1.],
       [2., 2., 2.]])

Here is the code for another way (that avoids taking the transpose). It is a little complicated, so we will break down below what was happening.

arr = np.zeros((3,3))
arr[:] = np.arange(3).reshape((3,1))
arr
array([[0., 0., 0.],
       [1., 1., 1.],
       [2., 2., 2.]])

Let’s look in particular at the right-hand side of arr[:] = np.arange(3).reshape((3,1)).

x = np.arange(3).reshape((3,1))
x
array([[0],
       [1],
       [2]])

Notice that these dimensions are compatible: for the right-most dimensions, we have 1 and 3, which are compatible (since one of them is equal to 1), and in the next dimension, we have 3 and 3, which are compatible (since they are equal).

x.shape
(3, 1)
arr.shape
(3, 3)

If we tried to reshape using (2,1), we would get an error, because we can’t fit 3 numbers into an array of shape (2,1).

y = np.arange(3).reshape((2,1))
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [36], in <cell line: 1>()
----> 1 y = np.arange(3).reshape((2,1))

ValueError: cannot reshape array of size 3 into shape (2,1)

A very convenient abbreviation to specify “just 1 column” is to use .reshape((-1,1)). Then NumPy is automatically going to replace the -1 with whatever positive integer is necessary to fill in all the numbers. In this case, the -1 gets replaced by 3.

z = np.arange(3).reshape((-1,1))
z
array([[0],
       [1],
       [2]])
z.shape
(3, 1)
arr
array([[0., 0., 0.],
       [1., 1., 1.],
       [2., 2., 2.]])

As another example of this -1 syntax, we get the following length-9 one-dimensional NumPy array. NumPy automatically replaced the -1 with a 9 in this case, because there were nine numbers in arr.

arr.reshape(-1)
array([0., 0., 0., 1., 1., 1., 2., 2., 2.])

Broadcasting in NumPy#

Because broadcasting is an important and potentially confusing topic in NumPy, here is an additional video on broadcasting, by Yasmeen Baki.

Timing a computation with a NumPy array vs with a list#

The main goal in this section is to see examples where code written using NumPy can run much faster than code written using base Python.

We’ll make some examples using random numbers, so we start out by instantiating a random number generator object.

rng = np.random.default_rng()

Now we make a length one million list of random integers from 1 to 10 (not including 10). We will be dividing by these numbers, so that’s why we don’t start with 0. The integers method of rng is a NumPy array.

arr = rng.integers(1,10,size=10**6)

Here is the list version of arr.

mylist = list(arr)

Here are the first 5 elements in arr.

arr[:5]
array([3, 7, 8, 2, 1])

Let’s compute 1.5 divided by each of these five numbers.

1.5/arr[:5]
array([0.5       , 0.21428571, 0.1875    , 0.75      , 1.5       ])

If we try to do the same thing with a list, we get an error, saying we can’t combine the float 1.5 with the list mylist using the operator /.

1.5/mylist[:5]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [7], in <cell line: 1>()
----> 1 1.5/mylist[:5]

TypeError: unsupported operand type(s) for /: 'float' and 'list'

Here is one way to make the list version, using a for loop.

newlist = []
for x in mylist[:5]:
    newlist.append(1.5/x)
newlist
[0.5, 0.21428571428571427, 0.1875, 0.75, 1.5]

We haven’t carefully studied list comprehension yet (we will next week), but just for contrast, here is how we can make this same list using list comprehension. The list comprehension is definitely more elegant (often said to be more “Pythonic”) than the for loop version.

[1.5/x for x in mylist[:5]]
[0.5, 0.21428571428571427, 0.1875, 0.75, 1.5]

Let’s see how the timings compare.

%%timeit
1.5/arr
1.12 ms ± 38.9 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
%%timeit
newlist = []
for x in mylist:
    newlist.append(1.5/x)
1.66 s ± 35.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%%timeit
[1.5/x for x in mylist]
1.62 s ± 27.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Notice how the NumPy version is over one thousand times faster than the list versions. My understanding is that there are two important reasons the NumPy version is much faster. One reason is that the NumPy version of division is vectorized, meaning all the divisions are performed at once (no looping is required). The other reason is that the entries in arr are attached to a specific data type, so no time is spent checking the data types of the elements in arr.

Also notice that, even though the list comprehension version is more elegant than the for loop version, the timings are very similar.

Counting using NumPy#

We’ve seen examples of making arrays using NumPy and of performing basic arithmetic with NumPy arrays. In this section, we’ll see some examples of counting using NumPy. For example, if we want to make a probability estimate using the formula “number of successes divided by number of experiments”, then being able to count efficiently will be very helpful.

Our example will involve random numbers. We’ll start with a small number (n = 20) of random numbers, and then later we will increase n. In general, you should imagine that the bigger the sample, the more accurate our probability estimate should be.

rng = np.random.default_rng()
n = 20
  • Make a NumPy array of random integers between -50 and 50 (inclusive). Also make a list with the same numbers.

arr = rng.integers(-50, 51, size=n)

Here is the list version of the array arr. (Be sure to always use a somewhat random variable name like mylist, as opposed to something like list, because list is already a built-in term in Python.)

mylist = list(arr)
  • How many of the integers are strictly less than 10?

Here is a for-loop approach to the above question.

count = 0
for x in mylist:
    if x < 10:
        count += 1
count
13

Here is the list comprehension version of making this count; again, notice how much more concise this is than the for-loop approach. (We will cover list comprehension more systematically next week… for now just think of these list comprehension examples as previews.) We get the same answer 13 as before.

len([x for x in mylist if x < 10])
13

To explain the NumPy approach to counting, it helps to actually look at the values in arr.

arr
array([-20,   7,  32, -46,  -2, -16, -12, -39, -22, -38,  17, -20,  49,
        27, -36,   0,  49,   3,  41,  39])

The expression arr < 10 creates what is called a Boolean array (meaning an array of Trues and Falses). One way to think about this is, it is broadcasting the < 10 portion to all of the entries in arr. Be sure you understand how the following Boolean array relates to the numbers in arr.

arr < 10
array([ True,  True, False,  True,  True,  True,  True,  True,  True,
        True, False,  True, False, False,  True,  True, False,  True,
       False, False])

There are many different possible ways to count True in this Boolean array. I believe the approach which is usually fastest is to use NumPy’s count_nonzero function. (Remember that True corresponds to 1 and that False corresponds to 0, so in a Boolean array, counting nonzero elements is the same as counting how often True occurs.) Again, we get the same answer 13 as above.

np.count_nonzero(arr < 10)
13
  • At what indices do these integers less than 10 occur?

We use the NumPy function np.nonzero to find indices where the value is True (compare this to np.count_nonzero which we found above). Notice that the resulting output is not a NumPy array… the output is a length-1 tuple which contains a NumPy array.

np.nonzero(arr < 10)
(array([ 0,  1,  3,  4,  5,  6,  7,  8,  9, 11, 14, 15, 17]),)
  • Make a new array containing only these integers.

We first show the for-loop approach.

newlist = []
for x in mylist:
    if x < 10:
        newlist.append(x)
newlist
[-20, 7, -46, -2, -16, -12, -39, -22, -38, -20, -36, 0, 3]

Here is the list comprehension aproach. We basically already did this approach above (the only difference is that we included a len wrapper above).

[x for x in mylist if x < 10]
[-20, 7, -46, -2, -16, -12, -39, -22, -38, -20, -36, 0, 3]
arr
array([-20,   7,  32, -46,  -2, -16, -12, -39, -22, -38,  17, -20,  49,
        27, -36,   0,  49,   3,  41,  39])
arr < 10
array([ True,  True, False,  True,  True,  True,  True,  True,  True,
        True, False,  True, False, False,  True,  True, False,  True,
       False, False])

The following is what is called “Boolean indexing” or “Applying a Boolean mask”. You should read arr[arr < 10] similarly to how you would read arr[:5]. Both are indicating which elements in arr to keep. In the case of arr[arr < 10], the part inside the brackets is a Boolean array, and we keep all the elements where True occurs.

arr[arr < 10]
array([-20,   7, -46,  -2, -16, -12, -39, -22, -38, -20, -36,   0,   3])
  • If you pick a random integer from such an array, what’s the probability it is less than 10?

To get an accurate probability estimate, we want to use as large a value of n as can be handled in a reasonably short amount of time; ten million works well in this case.

n = 10**7
arr = rng.integers(-50, 51, size=n)
mylist = list(arr)

Here is the for-loop approach. We think of count as the number of successes and we think of n as the number of experiments.

count = 0
for x in mylist:
    if x < 10:
        count += 1
count/n
0.5942684
%%timeit
count = 0
for x in mylist:
    if x < 10:
        count += 1
count/n
930 ms ± 134 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

The list comprehension takes a similar amount of time (it’s a little faster than the for-loop version).

%%timeit
len([x for x in mylist if x < 10])
717 ms ± 69.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Again, the NumPy version is much faster, over 100 times faster than the list comprehension version (and nearly 200 times faster than the for-loop version).

%%timeit
np.count_nonzero(arr < 10)
5.13 ms ± 235 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

To get a sense for how accurate our probability estimate is, we can run the estimate again with new random numbers and see how similar the probability is.

arr = rng.integers(-50, 51, size=n)

In both cases, we get a probability estimate of approximately 0.594.

np.count_nonzero(arr < 10)/n
0.5939019

There are 101 numbers from -50 to 50 (inclusive), and 60 of those are strictly less than 10, so the exact probability is 60/101. That exact value, 60/101, is very close to our probability estimates, so n = 10**7 in this case was a big enough sample size to get a very accurate probability estimate.

60/101
0.594059405940594

Logic in Python and NumPy#

General rule of thumb:

  • In base Python, use and, or, not

  • In NumPy (or pandas), use &, |, ~ (tilde)

rng = np.random.default_rng()
n = 20

arr = rng.integers(-50, 51, size=n)
mylist = list(arr)
mylist
[-20,
 17,
 -11,
 29,
 -4,
 -50,
 40,
 -36,
 23,
 24,
 -38,
 37,
 6,
 -32,
 -35,
 2,
 16,
 -32,
 -46,
 -47]
  • Find all the entries in arr and in mylist that are strictly between -10 and 10.

Here is the list comprehension solution. We use and because we are doing this in base Python.

[x for x in mylist if (x > -10) and (x < 10)]
[-4, 6, 2]

Here is a more complicated approach, that’s just here as an excuse to use not and or.

[x for x in mylist if not ((x <= -10) or (x >= 10))]
[-4, 6, 2]

Here is the NumPy version.

arr
array([-20,  17, -11,  29,  -4, -50,  40, -36,  23,  24, -38,  37,   6,
       -32, -35,   2,  16, -32, -46, -47])
arr > -10
array([False,  True, False,  True,  True, False,  True, False,  True,
        True, False,  True,  True, False, False,  True,  True, False,
       False, False])
arr < 10
array([ True, False,  True, False,  True,  True, False,  True, False,
       False,  True, False,  True,  True,  True,  True, False,  True,
        True,  True])

If we try to use and spelled out between these two Boolean arrays, we will get an error.

(arr > -10) and (arr < 10)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [10], in <cell line: 1>()
----> 1 (arr > -10) and (arr < 10)

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

We need to use the & to combine these Boolean arrays.

(arr > -10) & (arr < 10)
array([False, False, False, False,  True, False, False, False, False,
       False, False, False,  True, False, False,  True, False, False,
       False, False])

Here we find the values using Boolean indexing.

arr[(arr > -10) & (arr < 10)]
array([-4,  6,  2])

Just for fun, here are some other ways to get the same arrays. We use ~ to negate, and we use the vertical line | (sometimes called “pipe”) for or.

arr[(~(arr <= -10)) & (~(arr >= 10))]
array([-4,  6,  2])

The following is the analogue of the complicated list example from above. Notice how we replace not and or in the list version with ~ and | in this NumPy version.

arr[~((arr <= -10) | (arr >= 10))]
array([-4,  6,  2])

The axis keyword argument#

  • If you roll four distinct 6-sided dice, what is the probability that the biggest value is 5?

rng = np.random.default_rng()

Here is an example of simulating this experiment using a NumPy random number generator.

rng.integers(1, 7, size=4)
array([4, 3, 5, 1])

Let’s try this using ten experiments and a for loop. (We will eventually want to use many more experiments and we will want to avoid thee for loop.) Think of s as standing for the number of successes. In the following, we use np.max instead of the built-in Python function max. In general, if you have a choice between using a NumPy version or a built-in Python version, you should always prefer the NumPy version.

exps = 10
s = 0

for i in range(exps):
    if np.max(rng.integers(1, 7, size=4)) == 5:
        s += 1

s/exps
0.3

What we just did is very slow. Even for just one million experiments, the following took 13 seconds. (It is slow enough that I wanted to use %%time, which runs once, instead of %%timeit, which runs multiple times and takes an average.)

%%time

exps = 10**6
s = 0

for i in range(exps):
    if np.max(rng.integers(1, 7, size=4)) == 5:
        s += 1

s/exps
CPU times: user 13.7 s, sys: 164 ms, total: 13.9 s
Wall time: 14.3 s
0.284671

The key idea to making this computation more efficient is to work with a two-dimensional array, and to think of each row as corresponding to one experiment.

exps = 10
rng.integers(1, 7, size=(exps,4))
array([[1, 2, 3, 5],
       [5, 6, 1, 4],
       [2, 3, 4, 3],
       [5, 5, 5, 4],
       [4, 5, 4, 3],
       [3, 3, 5, 1],
       [6, 4, 2, 5],
       [3, 2, 2, 4],
       [3, 1, 5, 6],
       [2, 3, 4, 2]])

We can’t use np.max directly this context (at least, not without using the axis keyword argument), because np.max will find the maximum of all 40 numbers in this array.

exps = 10
arr = rng.integers(1, 7, size=(exps,4))
print(arr)
np.max(arr)
[[1 4 5 1]
 [6 5 4 4]
 [2 1 3 3]
 [6 6 1 6]
 [6 6 6 6]
 [1 4 2 2]
 [4 1 3 3]
 [5 1 4 5]
 [3 5 3 5]
 [3 4 6 1]]
6

As a first workaround, we could use NumPy’s apply_along_axis function. The documentation shows that the first three input arguments are: a function, an axis, and an array. In our case, the function will be np.max and the array will be arr. The main topic of this video is, what is meant by the axis keyword argument?

help(np.apply_along_axis)
Help on function apply_along_axis in module numpy:

apply_along_axis(func1d, axis, arr, *args, **kwargs)
    Apply a function to 1-D slices along the given axis.
    
    Execute `func1d(a, *args, **kwargs)` where `func1d` operates on 1-D arrays
    and `a` is a 1-D slice of `arr` along `axis`.
    
    This is equivalent to (but faster than) the following use of `ndindex` and
    `s_`, which sets each of ``ii``, ``jj``, and ``kk`` to a tuple of indices::
    
        Ni, Nk = a.shape[:axis], a.shape[axis+1:]
        for ii in ndindex(Ni):
            for kk in ndindex(Nk):
                f = func1d(arr[ii + s_[:,] + kk])
                Nj = f.shape
                for jj in ndindex(Nj):
                    out[ii + jj + kk] = f[jj]
    
    Equivalently, eliminating the inner loop, this can be expressed as::
    
        Ni, Nk = a.shape[:axis], a.shape[axis+1:]
        for ii in ndindex(Ni):
            for kk in ndindex(Nk):
                out[ii + s_[...,] + kk] = func1d(arr[ii + s_[:,] + kk])
    
    Parameters
    ----------
    func1d : function (M,) -> (Nj...)
        This function should accept 1-D arrays. It is applied to 1-D
        slices of `arr` along the specified axis.
    axis : integer
        Axis along which `arr` is sliced.
    arr : ndarray (Ni..., M, Nk...)
        Input array.
    args : any
        Additional arguments to `func1d`.
    kwargs : any
        Additional named arguments to `func1d`.
    
        .. versionadded:: 1.9.0
    
    
    Returns
    -------
    out : ndarray  (Ni..., Nj..., Nk...)
        The output array. The shape of `out` is identical to the shape of
        `arr`, except along the `axis` dimension. This axis is removed, and
        replaced with new dimensions equal to the shape of the return value
        of `func1d`. So if `func1d` returns a scalar `out` will have one
        fewer dimensions than `arr`.
    
    See Also
    --------
    apply_over_axes : Apply a function repeatedly over multiple axes.
    
    Examples
    --------
    >>> def my_func(a):
    ...     """Average first and last element of a 1-D array"""
    ...     return (a[0] + a[-1]) * 0.5
    >>> b = np.array([[1,2,3], [4,5,6], [7,8,9]])
    >>> np.apply_along_axis(my_func, 0, b)
    array([4., 5., 6.])
    >>> np.apply_along_axis(my_func, 1, b)
    array([2.,  5.,  8.])
    
    For a function that returns a 1D array, the number of dimensions in
    `outarr` is the same as `arr`.
    
    >>> b = np.array([[8,1,7], [4,3,9], [5,2,6]])
    >>> np.apply_along_axis(sorted, 1, b)
    array([[1, 7, 8],
           [3, 4, 9],
           [2, 5, 6]])
    
    For a function that returns a higher dimensional array, those dimensions
    are inserted in place of the `axis` dimension.
    
    >>> b = np.array([[1,2,3], [4,5,6], [7,8,9]])
    >>> np.apply_along_axis(np.diag, -1, b)
    array([[[1, 0, 0],
            [0, 2, 0],
            [0, 0, 3]],
           [[4, 0, 0],
            [0, 5, 0],
            [0, 0, 6]],
           [[7, 0, 0],
            [0, 8, 0],
            [0, 0, 9]]])
arr
array([[1, 4, 5, 1],
       [6, 5, 4, 4],
       [2, 1, 3, 3],
       [6, 6, 1, 6],
       [6, 6, 6, 6],
       [1, 4, 2, 2],
       [4, 1, 3, 3],
       [5, 1, 4, 5],
       [3, 5, 3, 5],
       [3, 4, 6, 1]])

Think of axis=1 as plugging in one row at a time, and think of axis=0 as plugging in one column at a time. In the following example, we get 5 as our first value, because 5 is the maximum of the top row. We get 6 as our next value, because 6 is the maximum of the next row, and so on.

np.apply_along_axis(np.max, axis=1, arr=arr)
array([5, 6, 3, 6, 6, 4, 4, 5, 5, 6])

What if we used axis=0 instead? Then we will get four 6s, because there are four columns, and each column has a 6 as its maximum.

np.apply_along_axis(np.max, axis=0, arr=arr)
array([6, 6, 6, 6])

There is a trick for getting “successes divided by experiments” as our probability estimate in this context. We make a Boolean array np.apply_along_axis(np.max, axis=1, arr=arr) == 5, and then we compute its average. The intuition is that you can think of the Boolean array as containing 1s and 0s (even though it actually contains Trues and Falses), and if we take its average, that is the same as “number of Trues divided by total number”.

exps = 10**6

arr = rng.integers(1, 7, size=(exps,4))
(np.apply_along_axis(np.max, axis=1, arr=arr) == 5).mean()
0.284891

This method is somewhat faster than our for loop method from above, but it still takes about 5 seconds, which is still fairly slow for only having one million experiments.

%%time

exps = 10**6

arr = rng.integers(1, 7, size=(exps,4))
(np.apply_along_axis(np.max, axis=1, arr=arr) == 5).mean()
CPU times: user 5.27 s, sys: 82.5 ms, total: 5.36 s
Wall time: 5.5 s
0.284142

The problem above is that np.apply_along_axis is very general, and is not optimized to the specific maximum computation we’re making. We will replace it by using the max method of arr. If we call max() on its own, we get the maximum of the whole array, which is not what we want.

exps = 10

arr = rng.integers(1, 7, size=(exps,4))
print(arr)
arr.max()
[[6 6 5 6]
 [5 1 5 5]
 [4 5 4 5]
 [1 3 2 3]
 [2 4 5 6]
 [1 1 1 5]
 [1 2 1 3]
 [2 5 6 1]
 [5 6 1 1]
 [2 1 2 6]]
6

If we check the documentation for arr.max (notice that you get documentation for this sort of method by calling help(arr.max), not help(arr.max() and not help(max)), we see that this max method also accepts an axis keyword argument, just like np.apply_along_axis.

help(arr.max)
Help on built-in function max:

max(...) method of numpy.ndarray instance
    a.max(axis=None, out=None, keepdims=False, initial=<no value>, where=True)
    
    Return the maximum along a given axis.
    
    Refer to `numpy.amax` for full documentation.
    
    See Also
    --------
    numpy.amax : equivalent function

Here is an example of using the axis keyword argument on a 10-row, 4-column NumPy array.

exps = 10

arr = rng.integers(1, 7, size=(exps,4))
print(arr)
arr.max(axis=1)
[[5 3 6 4]
 [6 3 1 1]
 [6 6 1 5]
 [1 6 6 2]
 [5 6 4 6]
 [3 1 5 2]
 [1 3 5 1]
 [6 2 4 2]
 [6 1 1 3]
 [2 6 1 5]]
array([6, 6, 6, 6, 6, 5, 5, 6, 6, 6])
exps = 10

arr = rng.integers(1, 7, size=(exps,4))
(arr.max(axis=1) == 5)
array([ True, False, False, False, False, False,  True,  True, False,
        True])

If we try to use axis=1 as an input argument to mean in this case, we will get an error, because the array is one-dimensional, so axis=1 does not exist.

exps = 10

arr = rng.integers(1, 7, size=(exps,4))
(arr.max(axis=1) == 5).mean(axis=1)
---------------------------------------------------------------------------
AxisError                                 Traceback (most recent call last)
Input In [21], in <cell line: 4>()
      1 exps = 10
      3 arr = rng.integers(1, 7, size=(exps,4))
----> 4 (arr.max(axis=1) == 5).mean(axis=1)

File ~/miniconda3/envs/math9f22/lib/python3.10/site-packages/numpy/core/_methods.py:168, in _mean(a, axis, dtype, out, keepdims, where)
    164 arr = asanyarray(a)
    166 is_float16_result = False
--> 168 rcount = _count_reduce_items(arr, axis, keepdims=keepdims, where=where)
    169 if rcount == 0 if where is True else umr_any(rcount == 0, axis=None):
    170     warnings.warn("Mean of empty slice.", RuntimeWarning, stacklevel=2)

File ~/miniconda3/envs/math9f22/lib/python3.10/site-packages/numpy/core/_methods.py:76, in _count_reduce_items(arr, axis, keepdims, where)
     74     items = 1
     75     for ax in axis:
---> 76         items *= arr.shape[mu.normalize_axis_index(ax, arr.ndim)]
     77     items = nt.intp(items)
     78 else:
     79     # TODO: Optimize case when `where` is broadcast along a non-reduction
     80     # axis and full sum is more excessive than needed.
     81 
     82     # guarded to protect circular imports

AxisError: axis 1 is out of bounds for array of dimension 1

We could do axis=0, but for a one-dimensional NumPy array, we might as well just not provide an axis argument.

exps = 10

arr = rng.integers(1, 7, size=(exps,4))
(arr.max(axis=1) == 5).mean(axis=0)
0.4

Here is the best version of our computation.

exps = 10**6

arr = rng.integers(1, 7, size=(exps,4))
(arr.max(axis=1) == 5).mean()
0.285484

The code is not only concise but it is also very efficient. Whereas even our NumPy code above took over 5 seconds, this is about one hundred times faster.

%%time
exps = 10**6

arr = rng.integers(1, 7, size=(exps,4))
(arr.max(axis=1) == 5).mean()
CPU times: user 46.5 ms, sys: 5.92 ms, total: 52.4 ms
Wall time: 51.3 ms
0.284911

More on the axis keyword argument#

Here is a supplemental video, by Yasmeen Baki, on the axis keyword argument. If you go on later to use the pandas library (for example, in Math 10), what you’ve learned about axis will become even more relevant.

An exact probability#

We ask the same question as last time, but now we will compute the exact probability.

  • If you roll four distinct 6-sided dice, what is the probability that the biggest value is 5?

Often computing the exact probability (at least using an exhaustive search where we check every possibility) is not feasible, but in this case, there are only \(6^4\) possible values (six different values for each of the four dice), so it’s no problem to compute all of the possibilities.

6**4
1296

Here is a correct but inelegant way to make a list of all possible quadruples of integers from 1 to 6 (including 6).

all_dice = []
for i in range(1,7):
    for j in range(1,7):
        for k in range(1,7):
            for l in range(1,7):
                all_dice.append((i,j,k,l))

As a reality check, this list does have the correct length.

len(all_dice)
1296

Here are the first ten values in all_dice. Do you see why it begins this way, and not for example [(1,1,1,1), (2,1,1,1), ...]?

all_dice[:10]
[(1, 1, 1, 1),
 (1, 1, 1, 2),
 (1, 1, 1, 3),
 (1, 1, 1, 4),
 (1, 1, 1, 5),
 (1, 1, 1, 6),
 (1, 1, 2, 1),
 (1, 1, 2, 2),
 (1, 1, 2, 3),
 (1, 1, 2, 4)]

What we really want to do is compute the Cartesian product of [1,2,3,4,5,6] with itself four times. The standard way to do that in Python is to use the itertools library, which is part of the standard Python libraries (meaning it does not need to be installed). We won’t import the whole module; we’ll just import its product function.

from itertools import product
help(product)
Help on class product in module itertools:

class product(builtins.object)
 |  product(*iterables, repeat=1) --> product object
 |  
 |  Cartesian product of input iterables.  Equivalent to nested for-loops.
 |  
 |  For example, product(A, B) returns the same as:  ((x,y) for x in A for y in B).
 |  The leftmost iterators are in the outermost for-loop, so the output tuples
 |  cycle in a manner similar to an odometer (with the rightmost element changing
 |  on every iteration).
 |  
 |  To compute the product of an iterable with itself, specify the number
 |  of repetitions with the optional repeat keyword argument. For example,
 |  product(A, repeat=4) means the same as product(A, A, A, A).
 |  
 |  product('ab', range(3)) --> ('a',0) ('a',1) ('a',2) ('b',0) ('b',1) ('b',2)
 |  product((0,1), (0,1), (0,1)) --> (0,0,0) (0,0,1) (0,1,0) (0,1,1) (1,0,0) ...
 |  
 |  Methods defined here:
 |  
 |  __getattribute__(self, name, /)
 |      Return getattr(self, name).
 |  
 |  __iter__(self, /)
 |      Implement iter(self).
 |  
 |  __next__(self, /)
 |      Implement next(self).
 |  
 |  __reduce__(...)
 |      Return state information for pickling.
 |  
 |  __setstate__(...)
 |      Set state information for unpickling.
 |  
 |  __sizeof__(...)
 |      Returns size in memory, in bytes.
 |  
 |  ----------------------------------------------------------------------
 |  Static methods defined here:
 |  
 |  __new__(*args, **kwargs) from builtins.type
 |      Create and return a new object.  See help(type) for accurate signature.

Here is an example of getting all possible (ordered) pairs of integers from 1 to 6. It returns an itertools product object, but we can convert that object into a list if we want to see all the pairs at once.

product(range(1,7), repeat=2)
<itertools.product at 0x112cd4900>
list(product(range(1,7), repeat=2))
[(1, 1),
 (1, 2),
 (1, 3),
 (1, 4),
 (1, 5),
 (1, 6),
 (2, 1),
 (2, 2),
 (2, 3),
 (2, 4),
 (2, 5),
 (2, 6),
 (3, 1),
 (3, 2),
 (3, 3),
 (3, 4),
 (3, 5),
 (3, 6),
 (4, 1),
 (4, 2),
 (4, 3),
 (4, 4),
 (4, 5),
 (4, 6),
 (5, 1),
 (5, 2),
 (5, 3),
 (5, 4),
 (5, 5),
 (5, 6),
 (6, 1),
 (6, 2),
 (6, 3),
 (6, 4),
 (6, 5),
 (6, 6)]

Here is another example of using the product function without using the repeat keyword argument. In this case, we’re getting the Cartesian product of the list ["a", "b", "c"] with the list [5, 6]. In total, this product has \(3 \cdot 2 = 6\) elements.

list(product(["a","b","c"], [5, 6]))
[('a', 5), ('a', 6), ('b', 5), ('b', 6), ('c', 5), ('c', 6)]

One more example.

list(product(["a","b","c"], repeat=3))
[('a', 'a', 'a'),
 ('a', 'a', 'b'),
 ('a', 'a', 'c'),
 ('a', 'b', 'a'),
 ('a', 'b', 'b'),
 ('a', 'b', 'c'),
 ('a', 'c', 'a'),
 ('a', 'c', 'b'),
 ('a', 'c', 'c'),
 ('b', 'a', 'a'),
 ('b', 'a', 'b'),
 ('b', 'a', 'c'),
 ('b', 'b', 'a'),
 ('b', 'b', 'b'),
 ('b', 'b', 'c'),
 ('b', 'c', 'a'),
 ('b', 'c', 'b'),
 ('b', 'c', 'c'),
 ('c', 'a', 'a'),
 ('c', 'a', 'b'),
 ('c', 'a', 'c'),
 ('c', 'b', 'a'),
 ('c', 'b', 'b'),
 ('c', 'b', 'c'),
 ('c', 'c', 'a'),
 ('c', 'c', 'b'),
 ('c', 'c', 'c')]

The above example, using repeat, is the exact same as this next example, which explicitly lists the list three times.

list(product(["a","b","c"], ["a","b","c"], ["a","b","c"]))
[('a', 'a', 'a'),
 ('a', 'a', 'b'),
 ('a', 'a', 'c'),
 ('a', 'b', 'a'),
 ('a', 'b', 'b'),
 ('a', 'b', 'c'),
 ('a', 'c', 'a'),
 ('a', 'c', 'b'),
 ('a', 'c', 'c'),
 ('b', 'a', 'a'),
 ('b', 'a', 'b'),
 ('b', 'a', 'c'),
 ('b', 'b', 'a'),
 ('b', 'b', 'b'),
 ('b', 'b', 'c'),
 ('b', 'c', 'a'),
 ('b', 'c', 'b'),
 ('b', 'c', 'c'),
 ('c', 'a', 'a'),
 ('c', 'a', 'b'),
 ('c', 'a', 'c'),
 ('c', 'b', 'a'),
 ('c', 'b', 'b'),
 ('c', 'b', 'c'),
 ('c', 'c', 'a'),
 ('c', 'c', 'b'),
 ('c', 'c', 'c')]

We’re now ready to use product to simulate every possible roll of four distinct six-sided dice.

temp_list = list(product(range(1,7), repeat=4))
len(temp_list)
1296

Let’s convert this to a NumPy array, so we can use some of the earlier methods from this week.

arr = np.array(temp_list)

The array arr has 1296 rows and 4 columns. We think of each row as corresponding to one possible outcome.

arr.shape
(1296, 4)

We want to know, for what proportion of the rows, is the maximum value in that row equal to 5? Because we want to get the maximum of each row, we use the axis=1 keyword argument. (If we wanted to get the maximum of each column, we would instead use axis=0.) We compute the total number of “successes” by using sum, and then we divide by the number of rows, which is arr.shape[0].

(arr.max(axis=1) == 5).sum()/arr.shape[0]
0.2847222222222222

An abbreviation is to apply the mean method to our Boolean array.

(arr.max(axis=1) == 5).mean()
0.2847222222222222

Just for fun, here is another computation of the exact probability. The first number in the difference, \((5/6)^4\), is equal to the probability that none of the dice lands 6, and then we are subtracting the probability that none of the dice lands 5 or 6.

(5/6)**4-(4/6)**4
0.2847222222222223

any and all#

Make a 100 row, 4 column NumPy array arr of random real numbers between 0 and 1.

rng = np.random.default_rng()
arr = rng.random(size=(100,4))
  • Find the sub-array of arr containing all rows in which at least one number is bigger than 0.9.

Here are the first 5 rows of arr.

arr[:5]
array([[0.02679582, 0.39964722, 0.38608334, 0.9557704 ],
       [0.62903247, 0.17981028, 0.48643802, 0.16059557],
       [0.49612872, 0.34270062, 0.00862007, 0.77873829],
       [0.55550408, 0.37485681, 0.20417908, 0.82782464],
       [0.51663663, 0.85479878, 0.48625127, 0.88475327]])

Here is a Boolean array corresponding to which elements are strictly greater than 5.

arr[:5] > 0.9
array([[False, False, False,  True],
       [False, False, False, False],
       [False, False, False, False],
       [False, False, False, False],
       [False, False, False, False]])

For each row, we want to know if any value in that row is True. To compute this, we will use the any method, together with axis=1. The reason we’re using axis=1 is because we are plugging in a full row at a time.

(arr[:5] > 0.9).any(axis=1)
array([ True, False, False, False, False])

If we wanted to count how many of these first 5 rows contain a number bigger than 0.9, we could use sum.

(arr[:5] > 0.9).any(axis=1).sum()
1

The question asked for us to find the rows which contain a number bigger than 0.9, so we will use Boolean indexing. We are plugging in the array([ True, False, False, False, False]) from above into the square brackets, which is like saying, keep the top row, discard the next four rows.

arr[:5][(arr[:5] > 0.9).any(axis=1)]
array([[0.02679582, 0.39964722, 0.38608334, 0.9557704 ]])

We now do the same thing to the full array. This will be a little easier to read, because we can get rid of the [:5] slices.

arr[(arr > 0.9).any(axis=1)]
array([[0.02679582, 0.39964722, 0.38608334, 0.9557704 ],
       [0.89128322, 0.97335335, 0.1821295 , 0.52423414],
       [0.4629412 , 0.71844233, 0.91362938, 0.04770196],
       [0.61030793, 0.15565872, 0.15616296, 0.9572168 ],
       [0.45523301, 0.85912337, 0.34660356, 0.91882733],
       [0.14055923, 0.86443333, 0.95413845, 0.94380163],
       [0.07850217, 0.96370185, 0.21801727, 0.22118187],
       [0.10730651, 0.56976149, 0.95367251, 0.7782287 ],
       [0.95705318, 0.11012904, 0.1380347 , 0.67057015],
       [0.95930433, 0.06501316, 0.67658178, 0.20159313],
       [0.2457492 , 0.67240127, 0.94358665, 0.90993476],
       [0.37975889, 0.92498652, 0.99327127, 0.7283782 ],
       [0.37093209, 0.45237422, 0.91160131, 0.96612184],
       [0.52254037, 0.19969189, 0.94795525, 0.15420271],
       [0.90339743, 0.59821255, 0.62045249, 0.54939576],
       [0.48895101, 0.52105202, 0.96297772, 0.45230437],
       [0.35246537, 0.32415082, 0.86599985, 0.98355888],
       [0.32443071, 0.94130683, 0.86154932, 0.95131523],
       [0.28175352, 0.71585491, 0.91062862, 0.47574106],
       [0.75785393, 0.97896958, 0.89668326, 0.81597407],
       [0.42821129, 0.43666654, 0.92245681, 0.6772787 ],
       [0.69686621, 0.34335879, 0.91188883, 0.73999112],
       [0.02840598, 0.96355532, 0.72067752, 0.77935551],
       [0.85410091, 0.06666216, 0.66487804, 0.93493077],
       [0.13820849, 0.59467061, 0.90928252, 0.64484179],
       [0.08076356, 0.54348621, 0.97440544, 0.7990639 ],
       [0.90641341, 0.55333912, 0.33929642, 0.77887452],
       [0.77605269, 0.45487112, 0.79606863, 0.90207317],
       [0.83243102, 0.34853769, 0.48317369, 0.99294516],
       [0.72748438, 0.95781736, 0.85863108, 0.78340871],
       [0.71794935, 0.97886636, 0.51022616, 0.41182274],
       [0.00766348, 0.94589262, 0.94978998, 0.86756495],
       [0.97386912, 0.84335624, 0.54012798, 0.47553105],
       [0.91223874, 0.48230536, 0.43958204, 0.51024251],
       [0.29992039, 0.89336646, 0.95530599, 0.90846783],
       [0.96276511, 0.46449629, 0.60265653, 0.55037618],
       [0.04358133, 0.68861017, 0.9216942 , 0.00731241],
       [0.57128587, 0.99787681, 0.77452267, 0.19917239],
       [0.91106315, 0.49769224, 0.30259767, 0.13776197],
       [0.76003629, 0.07852721, 0.97851938, 0.77227829]])
  • Find the sub-array of arr containing all rows in which no numbers are between 0.4 and 0.6.

One way to think about this is, we want the rows in which every element in that row is less than 0.4 or bigger than 0.6. (We allow a single row to have both smaller and bigger elements; all we care about is that it does not contain any values that start with 0.4 or 0.5.) Since we want a condition to be true for all the numbers in the row, we use all(axis=1).

((arr < 0.4) | (arr > 0.6)).all(axis=1)
array([ True, False, False, False, False, False, False, False, False,
       False, False,  True, False,  True,  True, False,  True,  True,
       False, False, False,  True,  True,  True,  True,  True,  True,
        True,  True, False, False, False, False, False, False, False,
       False, False,  True, False, False,  True, False,  True,  True,
       False, False, False, False, False,  True,  True,  True, False,
        True,  True, False,  True, False, False, False, False, False,
        True,  True,  True, False, False, False,  True, False, False,
        True, False,  True,  True, False, False, False, False, False,
        True, False, False, False, False, False, False,  True, False,
        True, False, False, False,  True,  True, False, False,  True,
       False])

As a reality check, our Boolean array here is length 100.

len(((arr < 0.4) | (arr > 0.6)).all(axis=1))
100

Here again we’re using Boolean indexing. We’re keeping all the rows of arr which satisfy this condition. As you scroll through this array, notice how all the values that appear are either less than 0.4 or bigger than 0.6.

arr[((arr < 0.4) | (arr > 0.6)).all(axis=1)]
array([[0.02679582, 0.39964722, 0.38608334, 0.9557704 ],
       [0.61030793, 0.15565872, 0.15616296, 0.9572168 ],
       [0.64463239, 0.62136197, 0.83285561, 0.69102624],
       [0.19806438, 0.28789025, 0.16912325, 0.71827727],
       [0.14055923, 0.86443333, 0.95413845, 0.94380163],
       [0.07850217, 0.96370185, 0.21801727, 0.22118187],
       [0.95705318, 0.11012904, 0.1380347 , 0.67057015],
       [0.6703015 , 0.81869015, 0.31666259, 0.67068789],
       [0.74644917, 0.06370154, 0.34252816, 0.06391344],
       [0.06099061, 0.08149432, 0.29021873, 0.85693158],
       [0.80698662, 0.31211537, 0.68117698, 0.78836747],
       [0.95930433, 0.06501316, 0.67658178, 0.20159313],
       [0.2457492 , 0.67240127, 0.94358665, 0.90993476],
       [0.37975889, 0.92498652, 0.99327127, 0.7283782 ],
       [0.35246537, 0.32415082, 0.86599985, 0.98355888],
       [0.32443071, 0.94130683, 0.86154932, 0.95131523],
       [0.75210216, 0.38678934, 0.19343926, 0.65172332],
       [0.75785393, 0.97896958, 0.89668326, 0.81597407],
       [0.69686621, 0.34335879, 0.91188883, 0.73999112],
       [0.02840598, 0.96355532, 0.72067752, 0.77935551],
       [0.81806092, 0.6881485 , 0.19411104, 0.05981003],
       [0.10944381, 0.11567344, 0.78604289, 0.88719194],
       [0.85410091, 0.06666216, 0.66487804, 0.93493077],
       [0.17557637, 0.17629594, 0.32083318, 0.78378524],
       [0.07667857, 0.75976036, 0.13388494, 0.8849134 ],
       [0.07788747, 0.17609963, 0.82659974, 0.23594508],
       [0.83519757, 0.04403135, 0.76341093, 0.12392983],
       [0.88621369, 0.82607382, 0.24463174, 0.02825251],
       [0.85579081, 0.00571833, 0.24021269, 0.70925784],
       [0.07239563, 0.0329126 , 0.23356791, 0.25781394],
       [0.72748438, 0.95781736, 0.85863108, 0.78340871],
       [0.00766348, 0.94589262, 0.94978998, 0.86756495],
       [0.29992039, 0.89336646, 0.95530599, 0.90846783],
       [0.72087445, 0.00735632, 0.66339639, 0.33424491],
       [0.04358133, 0.68861017, 0.9216942 , 0.00731241],
       [0.36901017, 0.06002681, 0.32674171, 0.26870498],
       [0.76003629, 0.07852721, 0.97851938, 0.77227829]])

From binary to decimal, version 1#

\[\begin{split} \begin{pmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 1 & 1 & 1 \end{pmatrix} \mapsto \begin{pmatrix} 10 \\ 4 \\ 15 \end{pmatrix} \end{split}\]
  • Write a function to_bin which takes in an m-by-n NumPy array of 0s and 1s, and as output returns a length m NumPy array of the corresponding integers, where we think of each row as representing the binary digits of an integer.

For example, we think of the top row as representing \(1 \cdot 2^3 + 0 \cdot 2^2 + 1 \cdot 2^1 + 0 \cdot 2^0 = 10\).

(If I were writing this again, I would call the function from_bin instead of to_bin, but I’ll stick with to_bin to match what’s in the video.)

Let’s start by getting a random array that can be our test input.

rng = np.random.default_rng()
arr = rng.integers(2, size=(10,4))

We assign m to be the number of rows and n to be the number of columns, using what is called tuple unpacking.

m,n = arr.shape
m
10
n
4

We want to multiply the left-most column by \(2^3\), the next column by \(2^2\), and so on. Here are the exponents. The middle -1 is telling NumPy where to stop (stop before -1), and the right-most -1 is telling NumPy what the step size should be (go down by 1 each step).

np.arange(n-1, -1, -1)
array([3, 2, 1, 0])

Here are the coefficients we want to multiply the columns by.

2**np.arange(n-1, -1, -1)
array([8, 4, 2, 1])

Using broadcasting, it’s very easy to make the multiplication. And to add up along the rows, we use sum(axis=1).

(arr*(2**np.arange(n-1, -1, -1))).sum(axis=1)
array([12,  7, 11, 15,  4,  1,  0, 12, 14,  9])
def to_bin(arr):
    _,n = arr.shape
    return (arr*(2**np.arange(n-1, -1, -1))).sum(axis=1)
to_bin(arr)
array([12,  7, 11, 15,  4,  1,  0, 12, 14,  9])

Here we check on a smaller array.

to_bin(np.array([[1,0],[1,1]]))
array([2, 3])

From binary to decimal, version 2#

\[\begin{split} \begin{pmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 1 & 1 & 1 \end{pmatrix} \mapsto \begin{pmatrix} 10 \\ 4 \\ 15 \end{pmatrix} \end{split}\]
  • Write a new function, to_bin2, which takes in an m-by-n NumPy array of 0s and 1s, and as output returns a length m NumPy array of the corresponding integers, where we think of each row as representing the binary digits of an integer.

rng = np.random.default_rng()
arr = rng.integers(2, size=(10,4))
arr
array([[0, 1, 1, 1],
       [0, 0, 1, 0],
       [1, 0, 0, 1],
       [1, 0, 1, 1],
       [0, 1, 1, 1],
       [1, 0, 0, 1],
       [0, 1, 0, 1],
       [0, 0, 1, 0],
       [1, 1, 0, 0],
       [0, 0, 1, 0]])

Here is the “usual” way of converting from a string to an integer.

int("0111")
111

But we can also tell Python that the string should be thought of as being written in binary. For example, 0111 in binary should be thought of as \(4+2+1 = 7\). In the following, we’re specifying that the number is being written in “base-2”.

int("0111", 2)
7

Our “usual” number system is written in base-10.

int("0111", 10)
111

Here we make a length-4 NumPy array.

z = np.array([0,1,1,1])

Or it’s probably more elegant to just take the top row of arr.

z = arr[0]
z
array([0, 1, 1, 1])

How can we convert that into the string “0111”? I feel like there might be an easier way that I’m missing, but this is what I’ve come up with. As a preliminary step, here is a reminder of list comprehension which we won’t officially cover in Math 9 until next week. In the following, we use list comprehension to put the integers 0, 1, 1, 1 into a list.

[x for x in z]
[0, 1, 1, 1]

Here we use list comprehension to put the strings “0”, “1”, “1”, “1” into a list. It is the same as the previous command, except that we convert each integer to a string using str.

[str(x) for x in z]
['0', '1', '1', '1']

I had initially thought the following would concatenate these strings (which would make sense, because "0"+"1" concatenates the strings to produce the string "01"), but it doesn’t work.

sum([str(x) for x in z])
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [12], in <cell line: 1>()
----> 1 sum([str(x) for x in z])

TypeError: unsupported operand type(s) for +: 'int' and 'str'

The following is the suggested way in Python to concatenate all the strings in a list together.

''.join([str(x) for x in z])
'0111'

We want to apply that idea to every row in arr. Let’s define a helper function to use. Later we will see a way to write this sort of short function in a more concise way, using “lambda functions”. (Caution. A common programming mistake is to use something like print instead of return. That’s very incorrect. If you want your function to have an output, you need to use return.)

def helper(z):
    return ''.join([str(x) for x in z])

Let’s test out our helper function on the 3rd row of arr.

helper(arr[3])
'1011'
arr[3]
array([1, 0, 1, 1])
arr
array([[0, 1, 1, 1],
       [0, 0, 1, 0],
       [1, 0, 0, 1],
       [1, 0, 1, 1],
       [0, 1, 1, 1],
       [1, 0, 0, 1],
       [0, 1, 0, 1],
       [0, 0, 1, 0],
       [1, 1, 0, 0],
       [0, 0, 1, 0]])

Here is a version helper2 that treats the row as if it corresponds to binary digits. We could do this all in one line, but it is more readable to break it up into two steps.

def helper2(z):
    temp = ''.join([str(x) for x in z])
    return int(temp, 2)
helper2(arr[3])
11

Here we’re using apply_along_axis again. The function we want to apply is helper2, and we are applying along rows, so we use axis=1.

def to_bin2(A):
    return np.apply_along_axis(helper2, axis=1, arr=A)

This doesn’t work because A is not defined.

to_bin2(A)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [21], in <cell line: 1>()
----> 1 to_bin2(A)

NameError: name 'A' is not defined
arr
array([[0, 1, 1, 1],
       [0, 0, 1, 0],
       [1, 0, 0, 1],
       [1, 0, 1, 1],
       [0, 1, 1, 1],
       [1, 0, 0, 1],
       [0, 1, 0, 1],
       [0, 0, 1, 0],
       [1, 1, 0, 0],
       [0, 0, 1, 0]])

Notice how, for example, the top row corresponds to 7, the next row corresponds to 2, and so on.

to_bin2(arr)
array([ 7,  2,  9, 11,  7,  9,  5,  2, 12,  2])

Raising errors#

  • Edit the to_bin function so that it raises an error if the input is not a NumPy array and it also raises an error if an entry is not 0 or 1.

rng = np.random.default_rng()
arr = rng.integers(2, size=(10,4))
arr
array([[1, 1, 1, 1],
       [0, 0, 1, 1],
       [1, 0, 1, 1],
       [1, 1, 1, 0],
       [0, 0, 1, 1],
       [0, 1, 0, 1],
       [1, 1, 0, 1],
       [1, 1, 1, 1],
       [1, 1, 0, 0],
       [0, 0, 1, 1]])
type(arr)
numpy.ndarray

Here is the usual way to check the type of a variable, using the function isinstance. It doesn’t literally check that the type is exactly a NumPy array, but what this checks is usually good enough. (Technically, it checks if arr is a NumPy array or is a subclass of the class NumPy array.)

isinstance(arr, np.ndarray)
True
isinstance(arr, list)
False

Here is how we raise an error in Python. In this case, we will raise a TypeError. Don’t worry about learning the many different error types in Python; for Math 9, the important thing is to know the general syntax of how to raise an error.

def to_bin(arr):
    if not isinstance(arr, np.ndarray):
        raise TypeError("Input should be a NumPy array")
    _,n = arr.shape
    return (arr*(2**np.arange(n-1, -1, -1))).sum(axis=1)
to_bin(arr)
array([15,  3, 11, 14,  3,  5, 13, 15, 12,  3])

Notice how the error message we wrote, “Input should be a NumPy array”, shows up in the following error message.

to_bin(list(arr))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [10], in <cell line: 1>()
----> 1 to_bin(list(arr))

Input In [8], in to_bin(arr)
      1 def to_bin(arr):
      2     if not isinstance(arr, np.ndarray):
----> 3         raise TypeError("Input should be a NumPy array")
      4     _,n = arr.shape
      5     return (arr*(2**np.arange(n-1, -1, -1))).sum(axis=1)

TypeError: Input should be a NumPy array

Now we will work on the second portion, namely, checking that all of the entries are 0 or 1.

arr == 0
array([[False, False, False, False],
       [ True,  True, False, False],
       [False,  True, False, False],
       [False, False, False,  True],
       [ True,  True, False, False],
       [ True, False,  True, False],
       [False, False,  True, False],
       [False, False, False, False],
       [False, False,  True,  True],
       [ True,  True, False, False]])
arr == 1
array([[ True,  True,  True,  True],
       [False, False,  True,  True],
       [ True, False,  True,  True],
       [ True,  True,  True, False],
       [False, False,  True,  True],
       [False,  True, False,  True],
       [ True,  True, False,  True],
       [ True,  True,  True,  True],
       [ True,  True, False, False],
       [False, False,  True,  True]])

The array is acceptable if each entry is zero or one, so we use the vertical line | for “or”.

(arr == 0) | (arr == 1)
array([[ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True]])

Because we are checking each entry (not each entry in a given row), we use the all method with no axis argument.

((arr == 0) | (arr == 1)).all()
True
def to_bin(arr):
    if not isinstance(arr, np.ndarray):
        raise TypeError("Input should be a NumPy array")
    if not ((arr == 0) | (arr == 1)).all():
        raise ValueError("All entries should be 0 or 1")
    _,n = arr.shape
    return (arr*(2**np.arange(n-1, -1, -1))).sum(axis=1)
to_bin(arr)
array([15,  3, 11, 14,  3,  5, 13, 15, 12,  3])
arr2 = arr

We haven’t emphasized the following type of assignment much, but here is how you set the 5th row, 2nd column to be 6.

arr2[5,2] = 6
arr2
array([[1, 1, 1, 1],
       [0, 0, 1, 1],
       [1, 0, 1, 1],
       [1, 1, 1, 0],
       [0, 0, 1, 1],
       [0, 1, 6, 1],
       [1, 1, 0, 1],
       [1, 1, 1, 1],
       [1, 1, 0, 0],
       [0, 0, 1, 1]])

Here this raises the error we chose, because one of the entries is 6.

to_bin(arr2)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [20], in <cell line: 1>()
----> 1 to_bin(arr2)

Input In [15], in to_bin(arr)
      3     raise TypeError("Input should be a NumPy array")
      4 if not ((arr == 0) | (arr == 1)).all():
----> 5     raise ValueError("All entries should be 0 or 1")
      6 _,n = arr.shape
      7 return (arr*(2**np.arange(n-1, -1, -1))).sum(axis=1)

ValueError: All entries should be 0 or 1

A big “gotcha” in Python is that our assignment above also changed the original arr array. (NumPy does this to save memory space.)

to_bin(arr)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [21], in <cell line: 1>()
----> 1 to_bin(arr)

Input In [15], in to_bin(arr)
      3     raise TypeError("Input should be a NumPy array")
      4 if not ((arr == 0) | (arr == 1)).all():
----> 5     raise ValueError("All entries should be 0 or 1")
      6 _,n = arr.shape
      7 return (arr*(2**np.arange(n-1, -1, -1))).sum(axis=1)

ValueError: All entries should be 0 or 1

Notice how arr got changed, even though we never changed it ourselves, we only changed arr2.

arr
array([[1, 1, 1, 1],
       [0, 0, 1, 1],
       [1, 0, 1, 1],
       [1, 1, 1, 0],
       [0, 0, 1, 1],
       [0, 1, 6, 1],
       [1, 1, 0, 1],
       [1, 1, 1, 1],
       [1, 1, 0, 0],
       [0, 0, 1, 1]])

If you want to make sure the original array doesn’t get changed, use the copy method.

arr3 = arr.copy()
arr3
array([[1, 1, 1, 1],
       [0, 0, 1, 1],
       [1, 0, 1, 1],
       [1, 1, 1, 0],
       [0, 0, 1, 1],
       [0, 1, 6, 1],
       [1, 1, 0, 1],
       [1, 1, 1, 1],
       [1, 1, 0, 0],
       [0, 0, 1, 1]])

Here we set all the even-indexed columns in row 6 to be equal to 99.

arr3[6,::2] = 99
arr3
array([[ 1,  1,  1,  1],
       [ 0,  0,  1,  1],
       [ 1,  0,  1,  1],
       [ 1,  1,  1,  0],
       [ 0,  0,  1,  1],
       [ 0,  1,  6,  1],
       [99,  1, 99,  1],
       [ 1,  1,  1,  1],
       [ 1,  1,  0,  0],
       [ 0,  0,  1,  1]])

Because we used the copy method when we created arr3, the original array arr did not get changed by our assignment.

arr
array([[1, 1, 1, 1],
       [0, 0, 1, 1],
       [1, 0, 1, 1],
       [1, 1, 1, 0],
       [0, 0, 1, 1],
       [0, 1, 6, 1],
       [1, 1, 0, 1],
       [1, 1, 1, 1],
       [1, 1, 0, 0],
       [0, 0, 1, 1]])