Numpy (numerical Python) is a library designed for performing scientific computing in Python.
In this notebook, we will introduce numpy arrays, a data structure introduced in numpy for working with vectors and matrices. We will explore how to create them, how to manipulate them, and how to use them for efficient numerical calculations using numpy functions.
Learning objectives for this notebook:
np.zeros(), np.linspace(), np.random.random()time()# Initialisation code for the notebook
import numpy as np
import matplotlib.pyplot as plt
import sys
sys.path.append('resource/asnlib/public/')
sys.path.append('../answer_checking/')
from validate_answers import *
Until now, the variable types we have been working with in Python represent relatively simple data types:
int: integer numbersfloat: floating point numberscomplex: complex-valued floating point numbersbool: boolean "truth" values (which can have True and False)str: strings list: list of variablesThe first four are really actually very simple data types, but actually the last one is more complicated than it looks. The list is a vector like-variable type. However, unlike physical vectors, it cannot be multiplied, subtracted, etc.
Here, we will introduce a new datatype that is handy for Physics calculations and that comes from the Python software package numpy called numpy arrays.
What are numpy arrays?
Numpy arrays are a way to work in Python with not just single numbers, but a whole bunch of numbers. With numpy arrays these numbers can be manipulated just like you do in, for example, linear algebra when one works with vectors:
https://en.wikipedia.org/wiki/Row_and_column_vectors
For example, a (column) vector $\vec{a}$ with 5 entries might look like this:
$$ \vec{a} = \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \end{bmatrix} $$In linear algebra you are used to manipulate these vectors, this can be done in a similar way with numpy arrays. We will use numpy arrays extensively in Python as vectors, like above, but also for storing, manipulating, and analyzing datasets (like a column of an excel spreadsheet).
To use numpy arrays, we first need to import the numpy library, which we will do using the shortened name "np" (to save typing):
import numpy as np
Now that we have imported numpy, we can use functions in numpy to create a numpy array. A simple way to do this is to use the function np.array() to make a numpy array from a comma-separated list of numbers in square brackets:
a = np.array([1,2,3,4,5])
print(a)
Note that numpy does not make a distinction between row vectors and column vectors: there are just vectors.
One useful thing about arrays is that you can access the elements of the array using square brackets:
a[n] will give you the n-th element of the array a.
This process of extracting a single element from the array is called indexing.
Note that here we encounter for the fist time what is know as the python counting from zero convention. What is the counting from zero convention? In the example above, we created a array:
a = np.array([1,2,3,4,5])The first element of a is 1. You might think that if you want to access the first element of a, you would use the notation a[1]. Right?
Let's try it:
print(a[1])
WRONG! Why? Because the makers of Python decided to start counting from zero: the first element of tuple a is actually a[0].
(This is a long-standing discussion among computer scientists, and the convention is different in many different languages. There are advantages and disadvantages of both, and even essays written about it...but in any case, Python chose to start arrays at zero.)
This also helps better understand the range() function: for example, to loop over all the elements in a, I can use this code:
for i in range(len(a)):
n = a[i]
print('a[%d] is %d' % (i,n))
Here the len function returns the length of the array a. As we saw before, Python has very smart for loops that can automatically iterate over many types of objects, which means we can also print out all the elements of our array like this:
for n in a:
print(n)
In Python, if you try to index beyond the end of the array, you will get an error:
a[5]
(Remember: indexing starts at zero!)
Python also has a handy feature: negative indices count backwards from the end, and index -1 corresponds to the last element in the array!
a[-1]
a[-2]
We can also use indexing to change the values of elements in our array:
print(a)
a[2] = -1
print(a)
Exercise 1 Set the first three, and the last two, entries of the following array to zero:
a = np.array([1,2,3,4,5,6,7,8,9,10,11,32,55,78,22,99,55,33.2,55.77,99,101.3])
#some code to set the first three and last two entries to zero
print(a)
Numpy arrays also support a special type of indexing called "slicing" that does not just return a single element of an array, but instead returns a whole part of array.
To do this, we put not just a single number inside my square brackets, but instead two numbers, separated by a colon :
a[n:m] will return a new tuple that consist of all the elements in a, starting at element n and ending at element m-1.
Let's look at a concrete example:
a = np.array(range(10))
print(a)
print(a[0:5])
The notation a[0:5] has "sliced" out the first five elements of the array.
With slicing, you can also leave off either n or m from the slice: if leave off n it will default to n=0, and if you leave off m, it will default to the end of the array (also the same as m=-1 in Python indexing):
print(a[:5])
print(a[5:])
Also handy: you can can have Python slice an array with a "step" size that is more than one by adding another : and a number after that. Find out its operation using:
print(a[0:10:2])
Fun, you can also use negative steps:
print(a[-1:-11:-1])
And finally, unlike indexing, Python is a bit lenient if you slice off the end of an array:
print(a[0:20])
Exercise 2: Slicing can also be used to set multiple values in an array at the same time. Use slicing to set first 10 entries of the array below to zero in one line of code.
a = np.array(range(20))+1
print(a)
# some code that sets the first 20 entries to zero
print(a)
An advantage of using numpy arrays for scientific computing is that way they behave under mathematical operations. In particular, they very often do exactly what we would want them to do if they were a vector:
a = np.array([1,2,3,4,5])
print(2*a)
print(a+a)
print(a+1)
print(a-a)
print(a/2)
print(a*a)
print(a**2)
Exercise 3: Generate an array y containing a sequence of the first 20 powers of 2 in a numpy array (starting at $2^0$).
Your output should be an array $[2^0, 2^1, 2^2, 2^3, ...]$.
(Hint: Start with a numpy array created using an appropriate range function that makes an array [0,1,2,3,...])
###
### Your code here
###
answer_2a_3_1 = np.copy(y)
question = "answer_2a_3"
num = 1
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
What about if I multiply two vectors together?
In mathematics, if I multiply two vectors, what I get depends on if I use the "dot product" or the "outer product" for my multiplication:
https://en.wikipedia.org/wiki/Row_and_column_vectors#Operations
The "dot product" corresponds to multiplying a column vector by a row vector to produce a single number. The "outer product" (also called the "tensor product") corresponds to multiplying the column vector by the row vector to make a matrix.
Question: If I type a*a, or more generally a*b, does Python use the inner or outer product?
It turns out: it uses neither! In Python, the notation a*a produces what is commonly called the "element-wise" product: specifically,
a*b = [a[0]*b[0] a[1]*b[1] a[2]*b[2] ...]
(Mathematically, this has a fancy name called the Hadamard product, but as you can see, despite the fancy name, it's actually very simple...)
We can see this in action here:
What if I actually want the dot product or the outer product? For that, Python has functions np.dot() and np.outer():
print(np.dot(a,a))
b=np.outer(a,a)
print(b)
Note that for matrix-matrix multiplication you can use dot for calculating the correct matrix, however matmul is preferred.
print(np.matmul(b,b))
Finally, adding somewhat to the confusion, correctly performing matrix-vector multiplication requires again the use of dot. See below:
print(np.dot(a,b))
This confusion of what method to use is mainly caused by the fact that arrays do not have the correct shape when specified as vectors. There are a few ways to solve this, one is to convert arrays to matrix objects with the numpy matrix class.
Pretty much all operators work with numpy arrays, even comparison operators, which can sometimes be very handy:
print(a)
print(a>3)
In numpy, there are also several handy functions for automatically creating arrays:
a = np.zeros(10)
print(a)
a = np.ones(10)
print(a)
np.linspace¶To automatically generate an array with linerly increasing values you can take np.linspace():
https://docs.scipy.org/doc/numpy/reference/generated/numpy.linspace.html
It takes three arguments: the starting number, the ending number, and the number of points.
This is a bit like the range function we saw before, but allows you to pick the total number of points, automatically calculating the (non-integer) step size you need:
a = np.linspace(0,20,40)
print(a)
print()
print("Length is: ", len(a))
print("Step size is: ", a[1]-a[0])
Note that if we wanted to have a step size of exactly 0.5, we need a total of 41 points:
a = np.linspace(0,20,41)
print(a)
print()
print("Length is: ", len(a))
print("Step size is: ", a[1]-a[0])
Exercise 4: Generate an array x that runs from -2 to 1 with 20 points using linspace.
###
### Your code here
###
answer_2a_4_1 = np.copy(x)
question = "answer_2a_4"
num = 1
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
np.arange()¶If we want to have more control on the exact spacing, we can use the np.arange() function. It is like range(), asking you for the start, stop, and step size:
a = np.arange(0,20,0.5)
print(a)
print()
print("Length is: ", len(a))
print("Step size is: ", a[1]-a[0])
Here, we already see a small quirk of arange: it will stop once the next point it calculates is < (not <=) to the stop point. If we want to get a range that stops at 20.0, we need to make the stop point any number a bit bigger than 20 (but smaller than our step size):
a = np.arange(0,20.00000001,0.5)
print(a)
print()
print("Length is: ", len(a))
print("Step size is: ", a[1]-a[0])
For this reason, np.arange() is not used very often, and mostly np.linspace() is used. There are also several other useful functions, such as np.geomspace(), which produces geometrically spaced points (such that they are evenly spaced on a log scale).
Exercise 5: Generate a numpy array z that has a first element with value 60 and last element 50 and takes steps of -0.5 between the values.
###
### Your code here
###
answer_2a_5_1 = np.copy(z)
question = "answer_2a_5"
num = 1
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
Numpy can also generate arrays of random numbers (we will see more of this later in the course when we deal with Monte Carlo methods):
a = np.random.random(40)
print(a)
This will generate uniform random numbers on the range of 0 to 1, but there are also several other random number generator functions that can make normally distributed random numbers, or random integers, and more.
Exercise 6: Generate a numpy array that contains 300 random grades that have a distribution of a bell-shaped curve that might represent the final grades of the students in this course, with an average grade of 7.5 and a standard deviation of 1. Make sure your grades are rounded to a half point.
(You may find it useful have to look at the help of the np.random.normal() function, and at your answer from Assignment 1 for the rounding.)
# Your solution here
Although we will not use them too much in this course, we can also use numpy also supports two-dimensional (or N-dimensional) numpy arrays, that can represent matrices. To make a 2D numpy array, you can use the zeros() function, for example, but with a two-entry list of numbers specifying the size N and M of the matrix:
m = np.zeros([10,10])
print(m)
For two dimensional matrices, the usual function len() is not enough to tell us about the shape of our matrix. Instead, we can use a property of the numpy matrix itself called its shape:
print(len(m))
print(m.shape)
Indexing two dimensional arrays works by using commas inside square brackets to specify the index of the first and second dimensions:
a = np.array(range(1,6))
m = np.outer(a,a)
print("Initial matrix:")
print(m)
# First index in the row number (counting from zero), second index in the column number
m[2,3] = -1
print("\nAfter changing entry [2,3] to -1:")
print(m)
You can also use slicing to to assign values to an array from a vector, which can be a handy way to enter a matrix by hand:
m = np.zeros([3,3])
m[0,:] = [1,2,3]
m[1,:] = [4,5,6]
m[2,:] = [7,8,9]
print(m)
Similarly, slicing also can be used to extract rows, columns, or blocks of the matrix:
# A row
print(m[1,:])
# A column
print(m[:,1])
# Extract a block as a sub-matrix
print(m[1:,1:])
There are several functions for making matrices (which we will see more of in this course when we are dealing with linear algebra):
https://docs.scipy.org/doc/numpy/reference/routines.array-creation.html
including this one which is used often:
# The identity matrix
print(np.eye(10))
# A band diagonal matrix
print(np.eye(10,k=-1))
Exercise 7: Use Python to calculate the multiplication of matrices m1 and m2:
and store the result in matrix m3. A matrix-matrix multiplication is not like the pointwise multiplication that you are familiar with for numpy arrays . Numpy has the function matmul() to multiply them.
m1 = np.zeros([3,3])
m2 = np.zeros([3,3])
###
### Your code here
###
answer_2a_7_1 = np.copy(m1)
answer_2a_7_2 = np.copy(m2)
answer_2a_7_3 = np.copy(m3)
question = "answer_2a_7"
num = 3
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
We can use for loops to iterate through numpy arrays and perform calculations.
For example, this code will calculate the average value of all the numbers in an array:
a = np.linspace(1,20,20)
avg = 0
for x in a:
avg += x
avg /= len(a)
print("Average is", avg)
Because this is a common operation, there is a function built into numpy np.average() that can do this for you:
print("Average is", np.average(a))
This is very handy: it saves us loads of typing! From the function name, it is also easy to understand what you are doing, making the code clearer and easier to read. However, the purpose of numpy functions is not only to save lots of typing: they also can often perform calculations MUCH faster than if you do program the calculation yourself with a for loop, as we will see in the next section.
Python also has many other useful functions for performing calculations using arrays:
# Calculate the standard deviation
print(np.std(a))
# The square root
print(np.sqrt(a))
# Numpy also has max and min, here is an example of min
a = np.linspace(-10,10,100)
print(np.min(a**2))
Good question for you to think about: why is this not zero? And what would I have to change above to get the code to return zero?
In addition to finding the minimum value in a vector, the function argmin can tell you where (what index number) the minimum is:
# Find the index number of the minimum of the array
i = np.argmin(a**2)
print(i)
print((a**2)[i])
Note also here that we used round brackets () around the a**2 in the print statement to be able to then index the resulting array a**2 array using square brackets [].
You can find the full list of mathematical numpy functions on the documentation website:
https://docs.scipy.org/doc/numpy/reference/routines.math.html
and the full list of all functions in the reference guide:
Exercise 8: Make a array x that runs from 0 to 4 with 20 points, and calculate an array y whose entries are equal to the square root of the entries in x.
###
### Your code here
###
answer_2a_8_1 = np.copy(y)
question = "answer_2a_8"
num = 1
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
In the first example above, we showed two ways of calculating an average: one using a for loop, and one using the numpy function.
Functionally, they are equivalent: they do exactly the same thing.
A curious feature of Python is that if you use functions instead of coding loops yourself, often things are MUCH MUCH faster.
To show this quantitatively, we will use the time library to calculate the time it takes to find the average of a pretty big array using both techniques:
# The time() function from the time library will return a floating point number representing
# the number of seconds since January 1, 1970, 00:00:00 (UTC), with millisecond or even microsecond
# precision
#
# We will use this to make a note of the starting time and the ending time,
# and then print out the time difference
from time import time
# A pretty big array, 50 million random numbers
a = np.random.random(int(50e6))
# Set timer
t1 = time()
# Caluclate the average
avg = 0
for x in a:
avg += x
avg /= len(a)
t2 = time()
t_forloop = t2-t1
print("The 'for' loop took %.3f seconds" % (t2-t1))
t1 = time()
avg = np.average(a)
t2 = time()
t_np = t2-t1
print("Numpy took %.3f seconds" % (t2-t1))
# Now let's compare them
print("\nNumpy was %.1f times faster!" % (t_forloop/t_np))
Why is numpy so much faster? The reason is that Python is an interpreted language. In each of the steps of the for loop, the Python kernel reads in the next step it has to do, translates that into an instruction for your computer processor, asks the computer to perform the step, gets the result back, reads in the next step, translates that into a processor instruction, sends that as an instruction to the computer processor, etc, etc.
If we did the same test in a compiled programing language like C, there would be no difference if we used a library function or if we wrote our own for loop.
When you use smart functions in Python libraries, like (many of) those in numpy, numpy will actually use an external library compiled in a language like C or Fortran that is able to send all of the calculation in one step to your computer processor, and in one step, get all the data back. This makes Python nearly as fast as a compiled language like C or Fortran, as long as you are smart in how you use it and avoid having "manual" for loops for large or long calculations.
(For small calculations, python coded for loops are perfectly fine and very handy!)
In the language of interpreted programmers, finding smart ways of getting what you need done using "compiled library functions" is often referred to as vectorisation.
Note that even normal mathematical operators are actually "vectorized functions" when they operate:
# This is actually a vectorized 'for' loop, it involves multiplying 50 million numbers by 5
b = 5*a
Here is a nice example of a vectorized way of counting the number of times the number '5' occurs in a random sample of 100 integers between 0 and 20:
nums = np.random.randint(0,21,100)
print("There are %d fives" % np.sum(nums == 5))
To see how this works, we can look at the intermediate steps:
nums = np.random.randint(0,21,100)
print(nums)
print(nums == 5)
print("There are %d fives" % np.sum(nums == 5))
Note that in this case, np.sum() will convert the bool value True into 1 and False into 0 for calculating the sum, according the the standard convertion of bool types to int types. You can see this in action if you want using the function astype() that is built into numpy arrays:
print(nums == 5)
print((nums == 5).astype('int'))