3  Homogeneous structures

Exercise 3.1 - Introduction to NumPy

Creation. The size and type of the elements of a NumPy array must be known in advance. The first way to create a NumPy array is to construct an array filled with zeros by specifying the size and type:

array0 = np.zeros(3, dtype=int) # vector of 3 integers
array1 = np.zeros((2,4), dtype=float) # array of floats of size 2x4
array2 = np.zeros((2,2), dtype=complex) # square matrix of complex numbers of size 2x2
array3 = np.zeros((5,6,4)) # three-dimensional array of floats

The second way is to pass the data directly:

array4 = np.array([1,4,5]) # vector of integers
array5 = np.array([[1.1,2.2,3.3,4.4],[1,2,3,4]]) # matrix of floats of size 2x4
array6 = np.array([[1+1j,0.4],[3,1.5]]) # matrix of complex numbers of size 2x2

NumPy will then determine itself the type and the size of the array. Note that it is possible to force the type:

array0 = np.array([1,4,5], dtype=complex) # vector of complex numbers

The type of the elements of the NumPy array array1 can be determined by array1.dtype. The size of this array is given by array1.shape. The following commands are used to access the array elements:

array4[1] # return 4
array5[1,3] # return 4.0

np.float64(4.0)

Note that the indices start at 0 and not at 1. NumPy arrays are mutable in the sense that the data can be modified while keeping the same type and size:

array0[1] = 4
array1[1,3] = 3.3
array3[3,4,2] = 3

Slicing. Slicing allows you to access certain parts of a table:

array4[2:3] # return the elements of indices between 2 and 3
array1[0,:] # return the first row of array1
array1[:,-1] # return the last column of array1
array3[3,3:5,1:4] # return the corresponding sub-matrix

array([[0., 0., 0.],
       [0., 3., 0.]])

Iteration. It is possible to iterate an array on its first dimension, for example, to return the sum of the rows:

for i in array5:
    print(np.sum(i))

11.0
10.0

a. Study the documentation for the arange function and use this function to generate the vectors \((5,6,7,8,9)\) and \((3,5,7,9)\).

b. Study the documentation for the function linspace and use it to generate 10 evenly spaced numbers in the interval \([2,5]\).

c. Read the documentation for the reshape function and perform the following transformations in succession:

\[(1,2,3,4,5,6)\to\begin{pmatrix}1 & 2\\ 3 & 4\\ 5 & 6 \end{pmatrix}\to\begin{pmatrix}1 & 2 & 3\\ 4 & 5 & 6 \end{pmatrix}\to\begin{pmatrix}1 & 4\\ 2 & 5\\ 3 & 6 \end{pmatrix}\]

Exercise 3.2 - Operations on arrays

The basic arithmetic operations on NumPy arrays are performed element by element:

mat1 = np.array([[1,2.5,3],[5,6.1,8],[3,2,5]])
mat2 = np.array([[1,0.5,0],[0,0.9,8],[2,0,0]])
mat1 + mat2 # return the sum element by element
mat1 * mat2 # return the product element by element (not the matrix product)
10*mat1**2 # return 10 times the square of the elements of mat1

array([[ 10. ,  62.5,  90. ],
       [250. , 372.1, 640. ],
       [ 90. ,  40. , 250. ]])

Most of the mathematical functions defined by NumPy (see https://numpy.org/doc/stable/reference/routines.math.html) are also performed element by element:

np.cos(mat1) # return the cosine element by element of mat1
np.exp(mat1) # return the exponential element by element of mat1

array([[2.71828183e+00, 1.21824940e+01, 2.00855369e+01],
       [1.48413159e+02, 4.45857770e+02, 2.98095799e+03],
       [2.00855369e+01, 7.38905610e+00, 1.48413159e+02]])

The matrix product can be performed in one of three ways:

np.dot(mat1,mat2)
mat1.dot(mat2)
mat1 @ mat2

array([[ 7.  ,  2.75, 20.  ],
       [21.  ,  7.99, 48.8 ],
       [13.  ,  3.3 , 16.  ]])

a. Given a vector \((v_0,v_1,\dots,v_{n-1})\), the discrete derivative of this vector is defined by the vector \((d_0,d_1,\dots,d_{n-2})\) given by \(d_i = v_{i+1}-v_{i}\) for \(i=0,1,\dots,n-2\).

Write a function diff_list that computes the discrete derivative of a list and a function diff_np that does the same operation but on NumPy vectors using slicing.

b. Let a_list and a_np be, respectively, a list and an array of 1 000 elements drawn at random in the interval \([0,1]\):

a_list = [np.random.random() for _ in range(1000)]
a_np = np.random.random(1000)

Compare the execution time of diff_list(a_list) and diff_np(a_np).

%%time
result = diff_list(a_list)

CPU times: user 42 μs, sys: 3 μs, total: 45 μs
Wall time: 48.4 μs

Exercise 3.3 - Vandermonde matrix

For \(p, n\in \mathbb{N}^*\) and \(\boldsymbol{x}= (x_1, \ldots, x_p)\) a vector of size \(p\), the corresponding Vandermonde matrix is defined by:

\[V(\boldsymbol{x},n)=\begin{pmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} & x_1^n \\ 1 & x_2 & x_2^2 & \cdots & x_2^{n-1} & x_2^n \\ \vdots & \vdots& \vdots & \ddots & \vdots & \vdots\\ 1 & x_{p-1} & x_{p-1}^2 & \cdots & x_{p-1}^{n-1} & x_{p-1}^n \\ 1 & x_p & x_p^2 & \cdots & x_p^{n-1} & x_p^n \end{pmatrix}.\]

a. Write a function that constructs the matrix \(V(\boldsymbol{x},n)\) element by element using a double loop.

b. After establishing a relationship to write the \(k\)-th column of \(V(\boldsymbol{x},n)\) solely as a function of \(\boldsymbol{x}\) and \(k\), write a second function that constructs the matrix \(V(\boldsymbol{x},n)\) column by column using this relationship.

c. After establishing a relationship between the \(k\)-th column of \(V(\boldsymbol{x},n)\), its \((k-1)\)-th column, and the vector \(\boldsymbol{x}\), write a third function that constructs the matrix \(V(\boldsymbol{x},n)\) column by column using this relationship.

d. Compare the execution times of these three functions for \(n=150\), \(p=100\), and \(\boldsymbol{x}\) generated randomly.

Exercise 3.4 - Array indexing !

Slicing allows you to select blocks in an array, but it is also possible to select disparate elements using an array as indexing:

a = np.arange(12)**2 # array of perfect squares
i = np.array([1,3,8,5]) # array of indices
a[i] # array of the elements with indices i

array([ 1,  9, 64, 25])

Note that it is also possible to index by an array of higher dimension. The result is then an array of the same shape as the index:

j = np.array([[3,4],[9,7]]) # two-dimensional array of indices
a[j] # select the elements with indices j

array([[ 9, 16],
       [81, 49]])

For a multi-dimensional array:

b = np.array([[0,1,2,3],[4,5,6,7],[8,9,10,11]])
i = np.array([0,1,2,2]) # array of first indices
j = np.array([1,0,3,1]) # array of second indices
b[i,j] # select the elements of indices ij

array([ 1,  4, 11,  9])

Finally, it is possible to index an array by an array of Booleans:

c = np.array([[0,1,2,3],[4,5,6,7],[8,9,10,11]])
cond = (c >= 5) # array of Booleans defined by True if >= 5 and False otherwise
c[cond] = 5 # assign the value 5 to all entries greater than 5

For the following, we consider the numbers:

[0.9602, -0.99, 0.2837, 0.9602, 0.7539, -0.1455, -0.99, -0.9111, 0.9602, -0.1455, -0.99, 0.5403, -0.99, 0.9602, 0.2837, -0.99, 0.2837, 0.9602]

[0.9602,
 -0.99,
 0.2837,
 0.9602,
 0.7539,
 -0.1455,
 -0.99,
 -0.9111,
 0.9602,
 -0.1455,
 -0.99,
 0.5403,
 -0.99,
 0.9602,
 0.2837,
 -0.99,
 0.2837,
 0.9602]

as the results of a measurement made every 0.1 second at times between 2 and 3.7 seconds.

a. Since the measurements are supposed to be positive, change the data to 0 when the values are negative.

b. Calculate the times for which the previous measurements are maximum.

c. For each maximum measure, return the previous measure, the maximum measure, and the next measure. If the previous or the next measure does not exist, replace them with np.nan.