Scientific and Engineering Libraries, Part 1¶
Weeks 12 and 13 are about using external scientific and engineering libraries that are not inside core Python, but are popular enough that we should learn them too.
Numpy¶
NumPy is a Python library for easy and efficient numerical computation of multi-dimensional arrays (vectors, matrices, etc.) You can install it via either pip or conda with, for example the following command:
pip install numpy
We first have to import numpy:
import numpy as np
Numpy works with lists (and lists of lists, etc.) and the core data type is something called numpy.ndarray, which we can create using the np.array() construction function:
l = [1, 2, 3]
arr = np.array(l)
type(arr)
numpy.ndarray
But why numpy arrays?¶
We already have lists, but lists are general-purpose containers of data. Numpy arrays are similar, but they are specifically useful for doing intuitive operations to lists of numbers.
As an example, consider we want to multiply each number in an array of numbers by some value:
l = [1, 2, 3]
l * 2
[1, 2, 3, 1, 2, 3]
As we see, this just concatenates the elements of the list. But with numpy arrays:
arr = np.array([1, 2, 3])
arr * 2
array([2, 4, 6])
Numpy arrays are designed in a way so that multiplying numbers within an array, adding arrays of numbers together, or similar common mathematical operations are both easy to write and also very efficient.
Multi-dimensional arrays¶
This is how we can create a 2D array, to represent something like a matrix:
arr = np.array([
[1, 2, 3, 4],
[5, 6, 7, 8]
])
A numpy array generally holds integers or floating point numbers, which can be checked by looking at its dtype:
arr.dtype
dtype('int64')
We can get the dimensions of a numpy array using .shape:
arr.shape # 2 rows, 4 columns
(2, 4)
.ndim gives the number of dimensions:
arr.ndim
2
A numpy array can be reshaped into a different shape using .reshape(new_shape). This will create a new array with the requested shape.
arr.reshape(8,1)
array([[1],
[2],
[3],
[4],
[5],
[6],
[7],
[8]])
arr.reshape(2, 2, 2) # 2x2x2 cube
array([[[1, 2],
[3, 4]],
[[5, 6],
[7, 8]]])
To get the number of elements inside an array, use .size:
arr.size
8
Similarly to lists in Python, we can access elements inside an array using the same square bracket syntax:
print(arr)
[[1 2 3 4] [5 6 7 8]]
arr[0][2]
np.int64(3)
As a special trick, we can access the same index as above using the following syntax:
arr[0, 2] # separated by commas
np.int64(3)
To get all of the elements of the second columns (counting from zero):
arr[:, 2]
array([3, 7])
To get all of the elements of the first row (counting from zero):
arr[1, :]
array([5, 6, 7, 8])
We can also create arrays that are automatically filled with a specific value:
np.zeros((10, 10))
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., 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.],
[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.]])
np.ones((5, 5))
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.]])
np.identity(8) # Create an 8x8 identity matrix
array([[1., 0., 0., 0., 0., 0., 0., 0.],
[0., 1., 0., 0., 0., 0., 0., 0.],
[0., 0., 1., 0., 0., 0., 0., 0.],
[0., 0., 0., 1., 0., 0., 0., 0.],
[0., 0., 0., 0., 1., 0., 0., 0.],
[0., 0., 0., 0., 0., 1., 0., 0.],
[0., 0., 0., 0., 0., 0., 1., 0.],
[0., 0., 0., 0., 0., 0., 0., 1.]])
# Create a matrix filled with random values sampled from normal distribution with mean 0 and standard deviation 1
np.random.normal(size=(4, 4))
array([[ 0.69542603, 0.79422213, 1.28925216, 1.49453288],
[-0.57540532, -0.01859384, 1.11260015, -3.17524228],
[ 0.9653842 , 1.46098111, 0.39489412, -1.55046235],
[-0.23291131, -1.95249727, 0.9793606 , 0.87044473]])
We can create arrays that with .arange, that mimics the behavior of range, as follows:
np.arange(1, 10, 2)
array([1, 3, 5, 7, 9])
list(range(1, 10, 2))
[1, 3, 5, 7, 9]
One difference is that .arange allows us to give a real-valued step size:
np.arange(1, 10, 0.2)
array([1. , 1.2, 1.4, 1.6, 1.8, 2. , 2.2, 2.4, 2.6, 2.8, 3. , 3.2, 3.4,
3.6, 3.8, 4. , 4.2, 4.4, 4.6, 4.8, 5. , 5.2, 5.4, 5.6, 5.8, 6. ,
6.2, 6.4, 6.6, 6.8, 7. , 7.2, 7.4, 7.6, 7.8, 8. , 8.2, 8.4, 8.6,
8.8, 9. , 9.2, 9.4, 9.6, 9.8])
For equally-spaced sampled points, we can also use .linspace:
np.linspace(2, 5, 10) # Sample 10 equally-spaced points from 2 to 5
array([2. , 2.33333333, 2.66666667, 3. , 3.33333333,
3.66666667, 4. , 4.33333333, 4.66666667, 5. ])
Operations between arrays¶
NumPy comes with lots of operations that can be done between arrays. The simplest is addition / subtraction:
a = np.array([
[1, 2],
[3, 4]
])
b = np.array([
[5, 1],
[3, 4]
])
a + b
array([[6, 3],
[6, 8]])
a - b
array([[-4, 1],
[ 0, 0]])
b - b
array([[0, 0],
[0, 0]])
We can also do logical comparison between arrays: These will be done element-wise, i.e, each element in the first array will be compared with the element in the second array that has the same index:
a < b
array([[ True, False],
[False, False]])
a > 2.5
array([[False, False],
[ True, True]])
There's also lots of ready-to-use functions that are in numpy, including:
np.sin(),np.cos(), ...np.min(),np.max(), ...np.mean(),np.std(),np.sum(), ...np.exp(),np.log(), ...
Most of these can also be called as an instance method (a.mean() vs np.mean(a)).
np.sin(0)
np.float64(0.0)
np.sin(np.pi / 2)
np.float64(1.0)
a
array([[1, 2],
[3, 4]])
np.min(a)
np.int64(1)
a.min()
np.int64(1)
a.mean()
np.float64(2.5)
def standardize(a):
return (a - a.mean()) / a.std()
standardize(a)
array([[-1.34164079, -0.4472136 ],
[ 0.4472136 , 1.34164079]])
Some operations can be specified to be done along an axis (i.e, row-wise, column-wise, etc.)
To do that, we specify the axis parameter:
a
array([[1, 2],
[3, 4]])
a.sum() # Default: sum up everything
np.int64(10)
a.sum(axis=0) # Sum up columns
array([4, 6])
a.sum(axis=1) # Sum up rows
array([3, 7])
a.min(axis=1)
array([1, 3])
a = np.arange(1, 37).reshape(6, 6)
a
array([[ 1, 2, 3, 4, 5, 6],
[ 7, 8, 9, 10, 11, 12],
[13, 14, 15, 16, 17, 18],
[19, 20, 21, 22, 23, 24],
[25, 26, 27, 28, 29, 30],
[31, 32, 33, 34, 35, 36]])
.hsplit and .vsplit are for partitioning an array into sub-arrays:
left, right = np.hsplit(a, 2) # Divide into left and right halves
left
array([[ 1, 2, 3],
[ 7, 8, 9],
[13, 14, 15],
[19, 20, 21],
[25, 26, 27],
[31, 32, 33]])
right
array([[ 4, 5, 6],
[10, 11, 12],
[16, 17, 18],
[22, 23, 24],
[28, 29, 30],
[34, 35, 36]])
top, bottom = np.vsplit(a, 2) # Divide into top and bottom halves
top
array([[ 1, 2, 3, 4, 5, 6],
[ 7, 8, 9, 10, 11, 12],
[13, 14, 15, 16, 17, 18]])
bottom
array([[19, 20, 21, 22, 23, 24],
[25, 26, 27, 28, 29, 30],
[31, 32, 33, 34, 35, 36]])
Counterpart to .hsplit and .vsplit are .hstack and .vstack for combining arrays into a single array.
np.hstack([left, right]) # Horizontally stack halves, reconstructing a
array([[ 1, 2, 3, 4, 5, 6],
[ 7, 8, 9, 10, 11, 12],
[13, 14, 15, 16, 17, 18],
[19, 20, 21, 22, 23, 24],
[25, 26, 27, 28, 29, 30],
[31, 32, 33, 34, 35, 36]])
np.vstack([top, bottom]) # Vertically stack parts, reconstructing a
array([[ 1, 2, 3, 4, 5, 6],
[ 7, 8, 9, 10, 11, 12],
[13, 14, 15, 16, 17, 18],
[19, 20, 21, 22, 23, 24],
[25, 26, 27, 28, 29, 30],
[31, 32, 33, 34, 35, 36]])
a
array([[ 1, 2, 3, 4, 5, 6],
[ 7, 8, 9, 10, 11, 12],
[13, 14, 15, 16, 17, 18],
[19, 20, 21, 22, 23, 24],
[25, 26, 27, 28, 29, 30],
[31, 32, 33, 34, 35, 36]])
A for loop with numpy arrays will loop over each row:
for row in a:
print('row:', row)
row: [1 2 3 4 5 6] row: [ 7 8 9 10 11 12] row: [13 14 15 16 17 18] row: [19 20 21 22 23 24] row: [25 26 27 28 29 30] row: [31 32 33 34 35 36]
We can transpose an array by using .T or .transpose(). Transposition makes it that each element at index (row, col) is now at index (col, row):
a
array([[ 1, 2, 3, 4, 5, 6],
[ 7, 8, 9, 10, 11, 12],
[13, 14, 15, 16, 17, 18],
[19, 20, 21, 22, 23, 24],
[25, 26, 27, 28, 29, 30],
[31, 32, 33, 34, 35, 36]])
a.T
array([[ 1, 7, 13, 19, 25, 31],
[ 2, 8, 14, 20, 26, 32],
[ 3, 9, 15, 21, 27, 33],
[ 4, 10, 16, 22, 28, 34],
[ 5, 11, 17, 23, 29, 35],
[ 6, 12, 18, 24, 30, 36]])
This lets us iterate over the rows of the original matrix, for example:
for col in a.T:
print('col:', col)
col: [ 1 7 13 19 25 31] col: [ 2 8 14 20 26 32] col: [ 3 9 15 21 27 33] col: [ 4 10 16 22 28 34] col: [ 5 11 17 23 29 35] col: [ 6 12 18 24 30 36]
Nested for loops allows us the "unravel" an array as follows:
for row in a:
for num in row:
print(num)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
NumPy provides .flat to do the same thing:
for i in a.flat:
print(i)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
Linear algebra with NumPy¶
NumPy also comes with a lot of tools for linear algebra.
As an example, suppose we have the following equations:
2x + 5y + z = 5
x + 2y + 3z = 3
3x + 3y + 5z = 10
We can find x, y, and z by representing the equations with the following matrices:
a = np.array([ # From the coefficients
[2, 5, 1],
[1, 2, 3],
[3, 3, 5]
])
b = np.array([5, 3, 10]) # From the results
np.linalg.solve will provide a solution if it exists:
x = np.linalg.solve(a, b)
x
array([ 3.63157895, -0.47368421, 0.10526316])
We can check the solution by computing the dot product with each row of the first matrix:
for row in a:
print(row.dot(x))
5.0 3.0 9.999999999999998
b
array([ 5, 3, 10])
In general, we can compute the inverse of a matrix (if it exists) with np.linalg.inv:
a
array([[2, 5, 1],
[1, 2, 3],
[3, 3, 5]])
a_inv = np.linalg.inv(a)
a_inv
array([[ 0.05263158, -1.15789474, 0.68421053],
[ 0.21052632, 0.36842105, -0.26315789],
[-0.15789474, 0.47368421, -0.05263158]])
The inverse of a matrix, when multiplied with original matrix, should produce an identity matrix:
np.matmul(a_inv, a)
array([[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 2.22044605e-16, 1.00000000e+00, 2.22044605e-16],
[ 0.00000000e+00, -5.55111512e-17, 1.00000000e+00]])
np.matmul(a, a_inv)
array([[ 1.00000000e+00, -5.55111512e-17, 1.66533454e-16],
[ 0.00000000e+00, 1.00000000e+00, 5.55111512e-17],
[-5.55111512e-17, 4.99600361e-16, 1.00000000e+00]])
There are many more functions that are available, such as:
np.linalg.det()for computing determinantsnp.trace()to get the trace of a matrixnp.linalg.eigh()to compute eigenvalues and eigenvectors for a given matrix.np.linalg.svd(),np.linalg.cholesky(), etc. to get matrix decompositions.
np.linalg.det(a)
np.float64(19.000000000000004)
np.trace(a)
np.int64(9)
np.linalg.eigh(a)
EighResult(eigenvalues=array([-0.35889894, 1. , 8.35889894]), eigenvectors=array([[ 5.54395106e-01, 7.07106781e-01, -4.38914646e-01],
[ 5.54395106e-01, -7.07106781e-01, -4.38914646e-01],
[-6.20719045e-01, -1.88336504e-17, -7.84033078e-01]]))
SciPy¶
SciPy is another external Python library for scientific computing. It has multiple libraries for various scientific computation tools, including but not limited to:
- Signal (and image) processing
- Numerical optimization
- Interpolation, statistics
- Clustering
See the textbook for further details about SciPy.