Tutorial 3 -- The NumPy module

The NumPy Module

NumPy is a well-documented, stable module where you will find important tools for scientific computing:

The documentation can be found here: NumPy User Guide

To import the module, we will use the following statement

import numpy as np

Creation/Manipulation of Arrays

NumPy allows you to easily create and manipulate arrays. To create sequences of numbers, NumPy provides the function arange analogous to range that returns arrays instead of lists.

import numpy as np

x = np.arange(0,2.0,0.1)   # creates array x[i]= i*0.1, i=0..19
print(np.size(x)) # print the size (dimension) of x
print(x[0])  # print first element
print(x[1]) # print second element
print(x[19]) # print last element
#print(x[20]) # This will print: IndexError: index 20 is out of bounds for axis 0 with size 20

a = np.array ([[1,2,3], [4,5,6], [7,8,9]]) # create 3x3 array (matrix)
b = 2 * a  # Multiplication of array by scalar: b[i,j] = 2*a[i,j]
c = a + b  # Summation of arrays: c[i,j]= a[i,j]+b[i,j]
d = np.dot(a, b) # Product of a by b: d[i,j]= a[i,k]*b[k,j] (summed over k)
e= a * b      # term by term product:  e[i,j]=a[i,j]*b[i,j] 
              # (this does not give the standard product of matrices)

You can convert a list to a numpy array. Numpy then allows you to apply mathematical functions on the array, such as scalar multiplication, as was already illustrated in the previous example:

Tcelcius = [20.1, 20.8, 21.9, 22.5, 20.9, 20.1] # this is a list of temperature values in Celcius
C = np.array(Tcelcius) # this is now a 1D array
print(C * 9 / 5 + 32)  # convert to Farhenheit

An interesting functionality in NumPy is its ability to slice through arrays to extract subarrays:

import numpy as np

t = np.array([1,2,3,4,5,6])
print(t[1:4])   # this gives t[i], i=1,2,3 (excluding both i=0 and i=4!)
print(t[:4])    # this prints t[i], i=0,1,2,3 (excluding i=4)
print(t[4:])    # this prints t[i], i=4 to the last index of array
print(t[:-1])   # this excludes the last element
print(t[1:-1])  # this excludes the first and last element

Another method to extract subsets of numpy arrays is to use indexation with a mask.

import numpy as np

a = np.arange(6)**2 # this gives a[i]= i**2, i=0,1,2,3,4,5
print(a > 10) # this gives a boolean array, of values True/False
print(a[a > 10]) # this gives a new array by extracting all values a[i]> 10

The copy of a NumPy array does copy the content of the array. Use .copy() to copy the array and its values.

import numpy as np
a = np.zeros((2, 2)) # The function zeros creates an array full of zeros
b = a
a[1, 1] = 10

c = b.copy()
b[1, 1] = 0

# The same rule applies to sliced arrays

d = np.arange(10)
e = d[:5]
d[0] = 10

The NumPy module offers lots of ways to create special arrays, to manipulate arrays, and to create operations on arrays.

import numpy as np
a = np.arange(10)
x= np.sum(a)
m= np.mean(a)

When arange is used with floating point arguments, it is not possible to predict the number of elements obtained, due to the finite floating point precision. Use the function linspace that receives as an argument the number of elements that we want, instead of the increment.

import numpy as np
from numpy import pi
a = np.linspace( 0, 2, 9 ) # creates an array of 9 elements from 0 to 2 (2 is included!)
x = np.linspace( 0, 2*pi, 100 ) # useful to evaluate function at lots of points
f = np.sin(x)  # create the array of elements f[i]= sin(x[i]), i=0..99

Random Number Generators

There is a fast and efficient way to generate pseudo-random numbers by using the random package of the Numpy module:

import numpy as np
a = np.random.random(10)

This will give you an array of 10 random float numbers in the interval $[0,1]$:

[ 0.25565389  0.0676878   0.74177835  0.2258494   0.03242723  0.34176716   0.5329745   0.42407123  0.40579087  0.97520252]
"<"class 'numpy.ndarray'">"

You can also obtain random integers with the NumPy package randint:

import numpy as np

outcome = np.random.randint(1, 7, size=10)

This will generate 10 integers randomly sampled between 1 and 6:

[6 3 1 6 4 5 1 2 5 5]

Application: A Monte Carlo evaluation of $\pi$

A simple to find a numerical approximation of $\pi$ is to use a Monte Carlo simulation: this consists in randomly selecting points $(x_i,y_i), i= 1, \ldots, n$ in the unit square $[0,1]\times [0,1] $ and in determining the ratio $m/n$, where $m$ is number of points that satisfy $x_i^2 + y_i^2 \le 1$. Here is a way to do this with numpy:

import numpy as np
N = 10000
x, y = np.random.random((2, N)) # create a 2xN array of random numbers in the interval [0,1]
mask = x**2 + y**2 < 1 # mask which defines the interior of 1/4 disk
approx_pi= 4*np.mean(mask) # mean(mask) gives the fraction of points inside the 1/4 disk
error = abs( (approx_pi-np.pi)/np.pi )
print('Relative error is',100*error,'%' )

[Previous Tutorial] [Next Tutorial]