## Overview and Learning Objectives

Python has several useful built-in packages as well as additional external packages. One such package is NumPy, which adds support for multi-dimensional arrays and matrices, along with a number of mathematical functions to operate on these structures. This unit covers array objects as the most fundamental NumPy data structure along with important NumPy array operations. Furthermore, we discuss numerical types, methods for type conversion, and constants (including the constants nan and inf) offered by NumPy. The three exercises included at the end of this unit cover key aspects needed throughout the PCP notebooks. Therefore, we encourage students to work through these exercises carefully. In particular, we recapitulate the mathematical concept of matrices and matrix multiplication, before we cover these aspects in Exercise 2 from an implementation perspective. We believe that understanding both—the mathematical formulation of abstract concepts and their realization in a programming context—is key in engineering education. This philosophy also forms the basis of the PCP notebooks to follow.

## NumPy Module¶

As said above, NumPy is a Python library used for working with arrays. In principle, one can also use the Python concept of lists to model an array. However, processing such lists is usually slow. Being mostly written in C or C++ and storing arrays at continuous places in memory (unlike lists), NumPy can process and compute with arrays in a very efficient way. This is the reason why NumPy is so powerful. In order to use NumPy, one needs to install the numpy package. This is what we already did when we created the PCP Environment, which contains a numpy-version. Furthermore, we need to use the Python import statement to get access to the NumPy functions. It is convenient to assign a short name to a frequently-used package (for example np in the case of numpy). This short name is used as a prefix when calling a function from the package. In the following code cell, we import numpy as np. To find out what the numpy-module contains, we use the Python command dir(np) to return a list all properties and methods of the specified module.

In [1]:
import numpy as np
list_numpy_names = dir(np)
print(list_numpy_names)

## NumPy Arrays¶

The fundamental NumPy data structure is an array object, which represents a multidimensional, homogeneous array of fixed-size items. Each array has a shape, a type, and a dimension. One way to create a NumPy array is to use the array-function. In the following code cells, we provide various examples and discuss basic properties of NumPy arrays.

In [2]:
import numpy as np

x = np.array([1, 2, 3, 3])
print('x = ', x)
print('Shape:', x.shape)
print('Type:', x.dtype)
print('Dimension:', x.ndim)
x =  [1 2 3 3]
Shape: (4,)
Type: int64
Dimension: 1

Note that, in this example, x.shape produces a one-element tuple, which is encoded by (4,) for disambiguation. (The object (4) would be an integer of type int rather than a tuple.) Two-dimensional arrays (also called matrices) are created like follows:

In [3]:
x = np.array([[1, 2, 33], [44, 5, 6]])
print('x = ', x, sep='\n')
print('Shape:', x.shape)
print('Type:', x.dtype)
print('Dimension:', x.ndim)
x =
[[ 1  2 33]
[44  5  6]]
Shape: (2, 3)
Type: int64
Dimension: 2

There are a couple of NumPy functions for creating arrays:

In [4]:
print('Array of given shape and type, filled with zeros: ', np.zeros(2))
print('Array of given shape and type, filled with integer zeros: ', np.zeros(2, dtype='int'))
print('Array of given shape and type, filled with ones: ', np.ones((2, 3)), sep='\n')
print('Evenly spaced values within a given interval: ', np.arange(2, 8, 2))
print('Random values in a given shape: ', np.random.rand(2, 3),  sep='\n')
print('Identity matrix: ', np.eye(3),  sep='\n')
Array of given shape and type, filled with zeros:  [0. 0.]
Array of given shape and type, filled with integer zeros:  [0 0]
Array of given shape and type, filled with ones:
[[1. 1. 1.]
[1. 1. 1.]]
Evenly spaced values within a given interval:  [2 4 6]
Random values in a given shape:
[[0.97071674 0.17657838 0.02515419]
[0.04013127 0.6477375  0.50146295]]
Identity matrix:
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]

## Array Reshaping¶

Keeping the total number of entries, there are various ways for reshaping an array. Here are some examples:

In [5]:
x = np.arange(2 * 3 * 4)
print('x =', x)
print('Shape:', x.shape)

y = x.reshape((3, 8))
print('y = ', y, sep='\n')
print('Shape:', y.shape)

z = np.reshape(x, (3, 2, 4))
print('z = ', z, sep='\n')
print('Shape:', z.shape)

print('Element z[0, 1, 2] = ', z[0, 1, 2])
x = [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23]
Shape: (24,)
y =
[[ 0  1  2  3  4  5  6  7]
[ 8  9 10 11 12 13 14 15]
[16 17 18 19 20 21 22 23]]
Shape: (3, 8)
z =
[[[ 0  1  2  3]
[ 4  5  6  7]]

[[ 8  9 10 11]
[12 13 14 15]]

[[16 17 18 19]
[20 21 22 23]]]
Shape: (3, 2, 4)
Element z[0, 1, 2] =  6

NumPy allows for giving one of the new shape parameters as -1. In this case, NumPy automatically figures out the unknown dimension. Note the difference between the shape (6,) and the shape (6,1) in the following example.

In [6]:
x = np.arange(6)
print(f'Shape: {x.shape}; dim: {x.ndim}')
x = x.reshape(-1, 2)
print(f'Shape: {x.shape}; dim: {x.ndim}')
x = x.reshape(-1, 1)
print(f'Shape: {x.shape}; dim: {x.ndim}')
Shape: (6,); dim: 1
Shape: (3, 2); dim: 2
Shape: (6, 1); dim: 2

## Array Operations¶

There are many ways to compute with NumPy arrays. Many operations look similar to the ones when computing with numbers. Applied to arrays, these operations are conducted in an element-wise fashion:

In [7]:
x = np.arange(5)
print('x = ', x)
print('x + 1 =', x + 1)
print('x * 2 =', x * 2)
print('x * x =', x * x)
print('x ** 3 =', x ** 3)
print('x / 4 =', x / 4)
print('x // 4 =', x // 4)
print('x > 2 =', x > 2)
x =  [0 1 2 3 4]
x + 1 = [1 2 3 4 5]
x * 2 = [0 2 4 6 8]
x * x = [ 0  1  4  9 16]
x ** 3 = [ 0  1  8 27 64]
x / 4 = [0.   0.25 0.5  0.75 1.  ]
x // 4 = [0 0 0 0 1]
x > 2 = [False False False  True  True]

Here are some examples for computing with matrices of the same shape. It is important to understand the difference between element-wise multiplication (using the operator *) and usual matrix multiplication (using the function np.dot):

In [8]:
a = np.arange(0, 4).reshape((2, 2))
b = 2 * np.ones((2, 2))
print('a = ', a, sep='\n')
print('b = ', b, sep='\n')
print('a + b = ', a + b, sep='\n')
print('a * b (element-wise multiplication) = ', a * b, sep='\n')
print('np.dot(a, b) (matrix multiplication) = ', np.dot(a, b), sep='\n')
a =
[[0 1]
[2 3]]
b =
[[2. 2.]
[2. 2.]]
a + b =
[[2. 3.]
[4. 5.]]
a * b (element-wise multiplication) =
[[0. 2.]
[4. 6.]]
np.dot(a, b) (matrix multiplication) =
[[ 2.  2.]
[10. 10.]]

Note that arrays and lists may behave in a completely different way. For example, using the operator + leads to the following results:

In [9]:
a = np.arange(4)
b = np.arange(4)
print(a + b, type(a + b))

a = list(a)
b = list(b)
print(a + b, type(a + b))
[0 2 4 6] <class 'numpy.ndarray'>
[0, 1, 2, 3, 0, 1, 2, 3] <class 'list'>

The sum of an array's elements can be computed along the different dimensions specified by the parameter axis. This is illustrated by the following example:

In [10]:
x = np.arange(6).reshape((2, 3))
print('x = ', x, sep='\n')
print('Total sum:', x.sum())
print('Column sum: ', x.sum(axis=0))
print('Row sum:', x.sum(axis=1))
x =
[[0 1 2]
[3 4 5]]
Total sum: 15
Column sum:  [3 5 7]
Row sum: [ 3 12]

There are many ways for accessing and manipulating arrays. Here are some examples:

In [11]:
x = np.arange(6).reshape((2, 3))
print('x = ', x, sep='\n')
print('Element in second column of second row:', x[1, 1])
print('Boolean encoding of element positions with values larger than 1: ', x > 1, sep='\n')
print('All elements larger than 1:', x[x > 1])
print('Second row:', x[1, :])
print('Second column:', x[:, 1])
x =
[[0 1 2]
[3 4 5]]
Element in second column of second row: 4
Boolean encoding of element positions with values larger than 1:
[[False False  True]
[ True  True  True]]
All elements larger than 1: [2 3 4 5]
Second row: [3 4 5]
Second column: [1 4]

## NumPy Type Conversion¶

In the PCP Notebook on Python Basics, we have already discussed the standard numeric types int, float, and complex offered by Python (and the function type() to identify the type of a variable). The NumPy package offers many more numerical types and methods for type conversion. We have already seen how to obtain the type of a numpy array using dtype. One can create or convert a variable to a specific type using astype. Some examples can be found in the next code cell.

In [12]:
a = np.array([-1,2,3])
print(f'a = {a}; dtype: {a.dtype}')
b = a.astype(np.float64)
print(f'b = {b}; dtype: {b.dtype}')
c = a.astype(np.int64)
print(f'c = {c}; dtype: {c.dtype}')
d = a.astype(np.uint8)
print(f'd = {d}; dtype: {d.dtype}')
e = np.array([1,2,3]).astype(np.complex64)
print(f'e = {e}; dtype: {e.dtype}')
a = [-1  2  3]; dtype: int64
b = [-1.  2.  3.]; dtype: float64
c = [-1  2  3]; dtype: int64
d = [255   2   3]; dtype: uint8
e = [1.+0.j 2.+0.j 3.+0.j]; dtype: complex64

Specifying the exact type is often important when using packages such as numba for optimizing machine code at runtime. In the following example, we give an example where a wrong initialization leads to an error (or warning with some nan results) when computing the square root of a negative number.

In [13]:
print('=== Initialization with \'int32\' leading to an error ===', flush=True)
x = np.arange(-2, 2)
print(x, x.dtype)
x = np.sqrt(x)
print(x)

print('=== Initialization with \'complex\' ===', flush=True)
x = np.arange(-3, 3, dtype='complex')
print(x, x.dtype)
x = np.sqrt(x)
print(x)
=== Initialization with 'int32' leading to an error ===
[-2 -1  0  1] int64
[nan nan  0.  1.]
=== Initialization with 'complex' ===
[-3.+0.j -2.+0.j -1.+0.j  0.+0.j  1.+0.j  2.+0.j] complex128
[0.        +1.73205081j 0.        +1.41421356j 0.        +1.j
0.        +0.j         1.        +0.j         1.41421356+0.j        ]
/tmp/ipykernel_114076/2135919221.py:4: RuntimeWarning: invalid value encountered in sqrt
x = np.sqrt(x)

## NumPy Constants¶

NumPy offers several constants that are convenient to compute with. In the following, we give some examples:

In [14]:
print(f'Archimedes constant Pi: {np.pi}; type: {type(np.pi)}')
print(f'Euler’s constant, base of natural logarithms: {np.e}; type: {type(np.e)}')
print(f'Floating point representation of (positive) infinity: {np.inf}; type: {type(np.inf)}')
print(f'Floating point representation of (negative) infinity: {np.NINF}; type: {type(np.NINF)}')
print(f'floating point representation of Not a Number (NaN): {np.nan}; type: {type(np.nan)}')
Archimedes constant Pi: 3.141592653589793; type: <class 'float'>
Euler’s constant, base of natural logarithms: 2.718281828459045; type: <class 'float'>
Floating point representation of (positive) infinity: inf; type: <class 'float'>
Floating point representation of (negative) infinity: -inf; type: <class 'float'>
floating point representation of Not a Number (NaN): nan; type: <class 'float'>

In particular, the constants nan and inf can be convenient to avoid case distinctions. However, computing with such constants can also be a bit tricky as the following examples show.

In [15]:
a = 10
b = np.inf
c = -np.inf
d = np.nan

print(f'a = {a}, b = {b}, c = {c}, d = {d}')
print('a + b =', a + b)
print('a * b =', a * b)
print('a + c =', a + c)
print('a - c =', a - c)
print('a + d =', a + d)
print('b + c =', b + c)
print('b * c =', b * c)
print('b / c =', b / c)
print('b + d =', b + d)
a = 10, b = inf, c = -inf, d = nan
a + b = inf
a * b = inf
a + c = -inf
a - c = inf
a + d = nan
b + c = nan
b * c = -inf
b / c = nan
b + d = nan

NumPy offers functions such as np.where and np.isinf to check for special constants. In the following, the class np.errstate is used to suppress warning that are usually output when dividing by zero.

In [16]:
print('Test element-wise for positive or negative infinity:', np.isinf([np.inf, np.NINF, np.nan]))

a = np.arange(-2, 3)
print('a = ', a)
with np.errstate(divide='ignore', invalid='ignore'):
b = a / 0
print('b = a / 0 =', b)

ind_inf = np.isinf(b)
ind_inf_pos = np.where(ind_inf)
print('Indices with inf values:', ind_inf, ind_inf_pos)

ind_nan = np.isnan(b)
ind_nan_pos = np.where(ind_nan)
print('Indices with nan values:', ind_nan, ind_nan_pos)
Test element-wise for positive or negative infinity: [ True  True False]
a =  [-2 -1  0  1  2]
b = a / 0 = [-inf -inf  nan  inf  inf]
Indices with inf values: [ True  True False  True  True] (array([0, 1, 3, 4]),)
Indices with nan values: [False False  True False False] (array([2]),)

## Exercises and Results¶

In [17]:
import libpcp.numpy
show_result = True

Exercise 1: NumPy Array Manipulations
• Create a NumPy array a with ascending natural numbers in the interval $[10, 20]=\{10,11,\ldots,20\}$ (using np.arange).
• Set all entries of a to zero where a$\leq13$ and a$>16$.
• Extend the resulting array a with a NumPy array containing the numbers of the interval $[4,6]$ and store the result in a variable b (using np.append).
• Remove duplicate values of b (using np.unique) and store the result in a variable c. Note that the result is automatically sorted in ascending order.
• Sort c in descending order and store the result in a variable d. Explore and discuss various options to do this including the slicing method c[::-1], the function reversed(), as well as the NumPy functions np.sort and np.flip.
In [18]:
#<solution>
#</solution>
In [19]:
libpcp.numpy.exercise_numpy_array(show_result=show_result)
[10 11 12 13 14 15 16 17 18 19 20]
[ 0  0  0  0 14 15 16  0  0  0  0]
[ 0  0  0  0 14 15 16  0  0  0  0  4  5  6]
[ 0  4  5  6 14 15 16]
[16 15 14  6  5  4  0] Type: <class 'numpy.ndarray'>
[16 15 14  6  5  4  0] Type: <class 'numpy.ndarray'>
[16 15 14  6  5  4  0] Type: <class 'numpy.ndarray'>
<reversed object at 0x7f6b89f84a60> Type: <class 'reversed'>
<reversed object at 0x7f6b89f848b0> Type: <class 'reversed'>

Recap: Matrices and Basic Operations
Let $N,M\in\mathbb{N}$ be two positive integers. An $(N\times M)$-matrix $\mathbf{A}$ is a rectangular array of entries arranged in rows and columns. For example, if entries are real numbers, one also writes $\mathbf{A}\in\mathbb{R}^{N\times M}$. Let $a_{nm}\in\mathbb{R}$ denote the entry of $\mathbf{A}$ being in the $n^{\mathrm{th}}$ row and $m^{\mathrm{th}}$ column for $n\in[1:N]=\{1,2,...,N\}$ and $m\in[1:M]$. Then, one also often writes $\mathbf{A}=[a_{nm}]$. A vector is a matrix where either $N=1$ (then called row vector) or $M=1$ (then called column vector).

When computing with vectors and matrices, one has to pay attention that the dimensions $N$ and $M$ of the operands match properly. For example, for a matrix summation or matrix subtraction, the operands must be of same dimensions or one of them needs to be a scalar (in this case the operations are applied per entry). The multiplication of matrices refers to a matrix multiplication. For an $(N\times M)$-matrix $\mathbf{A} = [a_{nm}]$ and a $(M\times P)$-matrix $\mathbf{B} = [b_{mp}]$, the product matrix $\mathbf{C} = \mathbf{A}\mathbf{B}$ with entries $\mathbf{C}=[c_{np}]$ is defined by $c_{np} = \sum_{m=1}^M a_{nm}b_{mp}$ for $n\in[1:N]$ and $p\in[1:P]$. In other words, the entry $c_{np}$ is the inner product (sometimes also called scalar product) of $n^{\mathrm{th}}$ row of $\mathbf{A}$ and $p^{\mathrm{th}}$ column of $\mathbf{B}$. This calculation rule is illustrated by the following figure:

Exercise 2: Matrix Operations
• Construct a matrix $\mathbf{B} = \begin{bmatrix}2 & 2 & 2 & 2\\2 & 2 & 2 & 2\\0 & 0 & 0 & 0\\\end{bmatrix}$ just using the NumPy functions np.zeros, np.ones, and np.vstack. Check the matrix dimensions using np.shape.
• Find the row and column index of the maximum entry of the matrix $\mathbf{D} = \begin{bmatrix}2 & 0 & 2\\-1 & 5 & 10\\-3 & 0 & 9\\\end{bmatrix}$ using the functions np.max, np.argmax and np.unravel_index. Why is it not sufficient to use np.argmax?
• Given a row vector $\mathbf{v} = [3\;2\;1]$ and a column vector $\mathbf{w} = [6\;5\;4]^T$, compute $\mathbf{v}\mathbf{w}$ and $\mathbf{w}\mathbf{v}$ using np.dot and np.outer. What is the difference between np.multiply, np.dot, and np.outer?
• Given $\mathbf{A} = \begin{bmatrix}1 & 2\\3 & 5\end{bmatrix}$, $\mathbf{v} = \begin{bmatrix}1\\4\end{bmatrix}$, compute $\mathbf{A}^{-1}$ and $\mathbf{A}^{-1}\mathbf{v}$ (using np.linalg.inv).
In [20]:
#<solution>
#</solution>
In [21]:
libpcp.numpy.exercise_matrix_operation(show_result=show_result)
Matrix B:
[[2. 2. 2. 2.]
[2. 2. 2. 2.]
[0. 0. 0. 0.]]
Shape of matrix B: (3, 4)
Maximum of matrix D: 10
Row and column index of maximum entry: (1, 2)
vw =  32
wv =
[[18 12  6]
[15 10  5]
[12  8  4]]
Inverse(A) =
[[-5.  2.]
[ 3. -1.]]
Inverse(A)v =
[ 3. -1.]

Exercise 3: Mathematical NumPy Functions
The NumPy package offers many different mathematical functions. Explore these functions by trying them out on small examples. In particular, you should gain a good understanding of the following concepts.

• Generate a NumPy array v_deg containing the numbers $[0, 30, 45, 60, 90, 180]$. Apply the function np.deg2rad to convert this array (given in degree) to an array v_rad (given in radiants). Then apply the functions np.cos and np.sin and discuss the results.
• Using the the same array as before, apply the exponential function np.exp(1j * v_rad). What is meant by 1j? What is the relation to np.cos and np.sin? Use the functions np.real, np.imag, and np.isclose to make this relation explicit.
• Try out different functions for rounding using the numbers $[-3.1416, -1.5708, 0, 1.5708, 3.1416]$. What is the difference between the functions np.round, np.floor, np.ceil, and np.trunc?
In [22]:
#<solution>
#</solution>
In [23]:
libpcp.numpy.exercise_numpy_math_function(show_result=show_result)
==== Cos and sin function ====
deg:  [  0  30  45  60  90 180]
rad:  [ 0.0000  0.5236  0.7854  1.0472  1.5708  3.1416]
cos:  [ 1.0000  0.8660  0.7071  0.5000  0.0000 -1.0000]
sin:  [ 0.0000  0.5000  0.7071  0.8660  1.0000  0.0000]

==== Exponential function ====
exp:  [ 1.00000000e+00+0.00000000e+00j  8.66025404e-01+5.00000000e-01j
7.07106781e-01+7.07106781e-01j  5.00000000e-01+8.66025404e-01j
6.12323400e-17+1.00000000e+00j -1.00000000e+00+1.22464680e-16j]
exp_real:  [ 1.0000  0.8660  0.7071  0.5000  0.0000 -1.0000]
exp_imag:  [ 0.0000  0.5000  0.7071  0.8660  1.0000  0.0000]
cos == exp_real: [ True  True  True  True  True  True]
sin == exp_imag: [ True  True  True  True  True  True]

==== Rounding options ====
Original numbers:           [-3.1416 -1.5708  0.0000  1.5708  3.1416]
Round to integer:           [-3.0000 -2.0000  0.0000  2.0000  3.0000]
Round to three decimals:    [-3.1420 -1.5710  0.0000  1.5710  3.1420]
Return floor of numbers:    [-4.0000 -2.0000  0.0000  1.0000  3.0000]
Return ceil of numbers:     [-3.0000 -1.0000  0.0000  2.0000  4.0000]
Return truncated numbers:   [-3.0000 -1.0000  0.0000  1.0000  3.0000]