Reference: Chapter 1 of Computational Nuclear Engineering and Radiological Science Using Python, R. McClarren (2018)
After studying this notebook, completing the activities, and asking questions in class, you should be able to:
So far, we have discussed lists, which are containers that hold any type of Python object such as strings, floats, integers, other lists, etc.
As engineers, we often want to store floating point numbers in vectors and matrices to perform linear algebra calculations. In this class, we will use NumPy. For a less engineering specific tutorial on NumPy, see https://docs.scipy.org/doc/numpy-1.15.0/user/quickstart.html
The basic unit in numpy is a multi-dimensional array:
We'll start by loading the numpy module (a.k.a. library).
# Import the numpy module (a.k.a library) and give it the "nickname" np
# Warning: Do not modify this line. The autograde assumes numpy is loades as 'np'
import numpy as np
We can now write np.
to call functions in the numpy library (a.k.a. module). You'll see we place a few standard import statements at the top of most notebooks in this class.
# Create a vector
a_vector = np.array([1,2,3,4])
# Create a matrix
a_matrix = np.array([(1,2,3),(4,5,6),(7,8,9)])
# Print to the screen
print("The vector",a_vector)
print("The matrix\n",a_matrix)
Arrays have several attributes that you can use to find out information regarding the array you're working with.
#shape tells you the shape
print("The shape of a_vector is ", a_vector.shape)
print("The shape of a_matrix is ", a_matrix.shape)
#ndim tells you the dimensionality of an array
print("The dimension of a_vector is ", a_vector.ndim)
print("The dimension of a_matrix is ", a_matrix.ndim)
#size is the total number of elements =
#the product of the number of elements in each dimension
print("The size of a_vector is ", a_vector.size,"= ",
a_vector.shape[0])
print("The size of a_matrix is ", a_matrix.size,"=",
a_matrix.shape[0],"*",a_matrix.shape[1])
Notice that a_matrix.shape
is a list. We can access the number of rows with a_matrix.shape[0]
and the number of columns with a_matrix.shape[1]
.
# YOUR SOLUTION HERE
You can also re-size an array after you create it:
A = np.array([2,4,6,8])
print("A is now a vector",A,"\n")
A = A.reshape(2,2)
print("A is now a matrix\n",A)
print("\nNow let's shape A back to a vector")
print(A.reshape(4))
Notice how I needed to assign A
with the reshaped array. We'll talk more about this later.
We can resize my_array
by stacking the 1st row followed by the 2nd row followed by the 3rd row:
print("Before:")
print(my_array)
print("\nAfter stacking:")
print(my_array.reshape(6))
We can resize my_array
by stacking the 2nd column under the first column.
my_array.reshape(6,order='F')
I want to highlight two things in the above example:
order='F'
stacks by columns, which seems strange at first glance.Regarding 2., a quick check of the documentation, https://docs.scipy.org/doc/numpy/reference/generated/numpy.reshape.html, reveals that order='C'
is the default. Moreover, 'C'
if shorthand for C-style stacking and 'F'
if short for Fortan-style stacking. C and Fortran are two mature (read "old") and popular programming languages for scientific computing. They use different conventions for storing arrays in memory, which means different stacking styles. The main take away is, why in doubt, peek at the documentation.
The arange
is a Numpy variant of the range
we saw earlier.
#let's make a vector from 0 to 2*pi in intervals of 0.1
dx = 0.1
X = np.arange(0,2*np.pi,dx)
print(X)
Typically, however, the linspace
function is actually what you want.
X = np.linspace(start = 0, stop = 2*np.pi, num = 62)
print(X)
Notice how it starts and ends exactly where I told it to.
If you want an array of zeros or ones, NumPy has a function for you.
zero_vector = np.zeros(10) #vector of length 10
zero_matrix = np.zeros((4,4)) #4 by 4 matrix
print("The zero vector:",zero_vector)
print("The zero matrix\n",zero_matrix)
ones_vector = np.ones(10) #vector of length 10
ones_matrix = np.ones((4,4)) #4 by 4 matrix
print("The ones vector:",ones_vector)
print("The ones matrix\n",ones_matrix)
# YOUR SOLUTION HERE
You can also create an array with random entries between 0 and 1:
random_matrix = np.random.rand(2,3) #random 2 x 3 matrix
print("Here's a random 2 x 3 matrix\n",random_matrix)
print("\nAnother example")
#make a random array between two numbers
print(np.random.uniform(low=-5,high=6,size=(3,3)))
#make random integers
print("\nA third example")
print(np.random.randint(low=1,high=6,size=10))
Finally, you may want to create an identity matrix:
#3 x 3 identity matrix
identity33 = np.identity(3)
print(identity33)
https://docs.scipy.org/doc/numpy/reference/generated/numpy.vstack.html
https://docs.scipy.org/doc/numpy/reference/generated/numpy.hstack.html
$$ M = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix}$$# YOUR SOLUTION HERE
NumPy also defines arithmetic operations mostly in the way that you would expect.
#vector addition
x = np.ones(3) #3-vector of ones
y = 3*np.ones(3)-1 #3-vector of 2's
print(x,"+",y,"=",x+y)
print(x,"-",y,"=",x-y)
Multiplication and division are element-wise, which may not be what you expect, but it is helpful.
x
y = np.array([1.0,2.0,3.0])
print(x,"*",y,"=",x*y)
print(x,"/",y,"=",x/y)
If you want the dot product, you have to use the dot
function:
print(x,".",y,"=",np.dot(x,y))
x = np.array([1,2.5,5])
y = np.array([4,6,7,9.0])
print(x+y)
my_answer = 'ehferbgrweavurebfvbidbv oisrevfbsrwiuevbgirwosgbvisbv . wefigwiuegfiuw'
Matrices work about the same way
silly_matrix = np.array([(1,2,3),(1,2,3),(1,2,3)])
print("The sum of\n",identity33,"\nand\n",
silly_matrix,"\nis\n",identity33+silly_matrix)
Multiplication and division are element-wise:
identity33 * silly_matrix
identity33 / silly_matrix
The dot
function will give you the matrix product:
print("The matrix product of\n",identity33,
"\nand\n",silly_matrix,"\nis\n",np.dot(identity33,silly_matrix))
#matrix times a vector
y = np.array([1.0,2.0,3.0])
print(silly_matrix,"times", y, "is")
print(np.dot(silly_matrix,y))
Writing np.dot
can get tedious. You can also use the @
symbol for matrix multiplication:
print(silly_matrix @ y)
NumPy also provides universal functions that operate on each element of an array. Common mathematical functions are defined in this way.
#recall we defined X as a linspace from 0 to 2pi
print(X)
# Print an empty line
print("\n")
#taking the sin(X) should be one whole sine wave
print(np.sin(X))
import matplotlib.pyplot as plt
#the next line is only needed in iPython notebooks
%matplotlib inline
# Create a plot. We'll learn more about this below
plt.plot(X,np.sin(X));
plt.axis([0,2*np.pi,-1.05,1.05])
plt.show()
import matplotlib.pyplot as plt
#the next line is only needed in iPython notebooks
%matplotlib inline
# Create a plot. We'll learn more about this below
plt.plot(X,np.cos(X));
plt.axis([0,2*np.pi,-1.05,1.05])
plt.show()
When you assign a new variable name to an existing array, it is the same as giving two names for the object.
a = np.array([1.0,2,3,4,5,6])
print(a)
b = a # this will make a and b different names for the same array
b[2] = 2.56 #change b at position 2, also changes a
print("The value of array a is",a)
print("The value of array b is",b)
Warning to MATLAB users: MATLAB would not change a
with the assignment b[2] = 2.56
. Python and MATLAB behave differently in this respect.
When you pass an array to a function, the function does not copy the array, it just assigns that array another name as in the previous example.
a1 = ["apples","oranges"]
b1 = a1
b1[1] = "grapes"
print(a1)
print(b1)
def devious_function(func_array):
func_array[0] = -1.0e6 #changes the value of array passed in
a = np.array([1.0,2,3,4,5,6])
print("Before the function a =",a)
devious_function(a)
print("After the function a =",a)
def less_devious_function(func_float):
func_float = 9.9 #changes local copy of func_float
a = 1
print("Before the function a =",a)
less_devious_function(a)
print("After the function a =",a)
This is different than what we saw previously for passing floats, integers, booleans, and strings to functions.
Why would Python do this, besides to just confuse me!?! you ask.
Imagine a NumPy array has millions of elements. It is best for efficient memory usage if functions do not make multiple copies of millions of elements.
Summary:
.copy()
. These are called primative types in Python..copy()
.To make a real copy of the entire elements of an array you need to use the copy
function on the array.
a = np.array([1.0,2,3,4,5,6])
print(a)
#this will make a and b different copies for the same array
b = a.copy()
b[2] = 2.56 #change b at position 2, will not change a
print("\nThe value of array a is",a)
print("\nThe value of array b is",b)
Oftentimes we'll want to access more than one element of an array at a time. We can do this by slicing.
#bring these guys back
import numpy as np
a_vector = np.array([1,2,3,4])
a_matrix = np.array([(1,2,3),(4,5,6),(7,8,9)])
print("The vector",a_vector)
print("\nThe matrix\n",a_matrix)
print("\na_vector[:] =",a_vector[:]) #single colon gives everything
print("\na_vector[1:3] =",a_vector[1:3]) #print out position 1 to position 2 (same as for lists)
print("\nFor a matrix, we can slice in each dimension")
print("a_matrix[0,:] =",a_matrix[0,:]) #every column in row 0
print("a_matrix[0,1:3] =",a_matrix[0,1:3]) #columns 1 and 2 in row 0
print("a_matrix[:,2] =",a_matrix[:,2]) #every row in column 2
Iteration over a matrix will give you everything in a row:
a_matrix = np.array([(1,2,3),(4,5,6),(7,8,9)])
count = 0
for i in a_matrix:
print("Row",count,"of a_matrix is",i)
count += 1
for column in a_matrix.transpose():
print("Column",count,"of a_matrix is",column)
count += 1
To iterate over every element in the matrix you'll need two for loops:
a_matrix = np.array([(1,2,3),(4,5,6),(7,8,9)])
row_count = 0
col_count = 0
for row in a_matrix:
col_count = 0
for col in row:
print("Row",row_count,"Column",col_count,"of a_matrix is",col)
col_count += 1
row_count += 1
Tip: Let my_list
be a Python list. Remember that my_list.append(1.0)
will add 1.0
to the list my_list
.
N = 3
A = np.random.rand(N,N)
print("A =\n",A)
# create any empty list
diagonal = []
#print the diagonal of the matrix using a single for loop
print("\nThe diagonal elements of A")
# YOUR SOLUTION HERE
Warning: You need to programmically extract the diagonals. You cannot manually copy them. When the autograder runs, it will generate a new A
matrix.
NumPy also has a built-in function to extract the diagonal of a matrix.
print("\nThe diagonal of A as a vector")
print(np.diag(A))
Numpy can handle complex numbers, but you need to tell it that you have complex numbers by adding a parameter to the function that creates your array (or matrix, or vector, etc.).
cArray = np.zeros(5, dtype="complex")
print(cArray)
dtype
tells NumPy that the array will hold complex numbers. It is not needed if it is obvious you are using complex numbers.
cArray2 = np.array([1+1j,-1])
print(cArray2)
You must be careful when you manipulate an array of floats into complex numbers, say by taking the square root
fArray = np.array([-1,1])
print("Wrong way gives", np.sqrt(fArray))
print("Right way gives", np.sqrt(fArray, dtype="complex"))
# YOUR SOLUTION HERE
def np_activity(N):
'''Generates a random N x N matrix, prints to screen,
and returns a vector that is the sum of the columns.
Arguments:
YOU NEED TO FILL THIS IN
Returns:
YOU NEED TO FILL THIS IN
'''
Discussion Question: what would we need to change to access the matrix we created outside our function?