Scientific python: NumPy, SciPy and Matplotlib¶

Learning objectives¶

Learn what NumPy arrays are
Learn basic array manipulations
Learn what vectorial code is
Quick overview of SciPy
Learn how to do a simple 2D plot and decorate it
Learn how to combine plots into a single figure

SciPy ecosystem¶

From their website:

SciPy (pronounced “Sigh Pie”) is a Python-based ecosystem of open-source software for mathematics, science, and engineering. In particular, these are some of the core packages:

NumPy

SciPy library

Matplotlib

IPython

Sympy

Pandas

…and much more.

NumPy¶

NumPy is a key building block used by SciPy.
NumPy supports:
- Multidimensional arrays (ndarray)
- Matrices and linear algebra operations
- Random number generation
- Fourier transforms
- Polynomials
- Tools for integrating with Fortran/C libraries
NumPy provides fast precompiled functions for numerical routines.
Official website: https://www.numpy.org/

“Array programming with NumPy”¶

https://www.nature.com/articles/s41586-020-2649-2

Abstract: Array programming provides a powerful, compact and expressive syntax for accessing, manipulating and operating on data in vectors, matrices and higher-dimensional arrays. NumPy is the primary array programming library for the Python language. It has an essential role in research analysis pipelines in fields as diverse as physics, chemistry, astronomy, geoscience, biology, psychology, materials science, engineering, finance and economics. For example, in astronomy, NumPy was an important part of the software stack used in the discovery of gravitational waves and in the first imaging of a black hole. Here we review how a few fundamental array concepts lead to a simple and powerful programming paradigm for organizing, exploring and analysing scientific data. NumPy is the foundation upon which the scientific Python ecosystem is constructed. It is so pervasive that several projects, targeting audiences with specialized needs, have developed their own NumPy-like interfaces and array objects. Owing to its central position in the ecosystem, NumPy increasingly acts as an interoperability layer between such array computation libraries and, together with its application programming interface (API), provides a flexible framework to support the next decade of scientific and industrial analysis.

NumPy Arrays overview¶

Core (or Standard) Python Library provides lists and 1D arrays (array.array)
- Lists are general containers for objects
- Arrays are 1D containers for objects of the same type
- Limited functionality
- Some memory and performance overhead associated with these structures
NumPy provides multidimensional arrays (numpy.ndarray)
- Can store many elements of the same data type in multiple dimensions
- cf. Fortran/C/C++ arrays
- More functionality than Core Python e.g. many conveninent methods for array manipulation
- Efficient storage and execution
Extensive online documentation !

Creating 1D arrays¶

# Convention: Import numpy and alias/rename it as np
import numpy as np

There are many ways to create 1d array, e.g.:

# Create array from Python list.
a = np.array([-1, 0, 1])
print(a)

# Create array using the array class copy constructor.
b = np.array(a)
print(b)

While these might look the same when printed above, all Numpy arrays are of type ‘ndarray’:

print(type(b))
li = [1, 2, 3]
print(type(li))

Arrays can also be created using numpy functions, e.g.

# https://numpy.org/doc/stable/reference/generated/numpy.arange.html
# arange for arrays (like using range for lists)
a = np.arange( -2, 6, 2 )
print(a)

# https://numpy.org/doc/stable/reference/generated/numpy.linspace.html
# linspace to create sample step points in an interval
a = np.linspace(-10, 10, 4) 
print(a)

# Contrast this with arange.
print(np.arange(-10, 10, 4))

Exercise 2.1¶

Add comments explaining what the following NumPy functions do.

b = np.zeros(3)
print(b)

c = np.ones(3)
print(c)

d = np.eye(3)
print(d)

Array attributes¶

As part of the array structure, NumPy keeps track of metadata for the array as “attributes”

# Taking "a" from the previous example
a = np.linspace(-10, 10, 5) 

# Examine key array attributes
print(a)
print("Dimensions ", a.ndim)   # Number of dimensions
print("Shape      ", a.shape)  # number of elements in each dimension
print("Size       ", a.size)   # total number of elements
print("Data type  ", a.dtype)  # data type of element, 64 bit float (IEEE 754) by default

Basic data type can be specified at creation:

a = np.array( [1.1,2.2,3.3], np.float32)
print(a)
print("Data type", a.dtype)

a = np.array([1,2,3,4], np.float64)
print(a)
print("Data type", a.dtype)

Multi-dimensional arrays¶

There are many different ways to create N-dimensional arrays. A two-dimensional array or matrix can be created from, e.g., list of lists

mat = np.array( [[1,2,3], [4,5,6]] )
print(mat)
print("Dimensions: ", mat.ndim)
print("Size:       ", mat.size)
print("Shape:      ", mat.shape)

pprint can be used for fancier display of multi-dimensional arrays:

from pprint import pprint

pprint(mat)

You can create 2d arrays with complex elements by specifying the data type

alist = [[1, 2, 3], [4, 5, 6]]
mat = np.array(alist, complex)
pprint(mat)

Exercise 2.2¶

Work out the shape of the resulting arrays before executing the following cells
(HINT: length of the innermost list gives the size of the rightmost shape index)

i = np.array( [[1,1,1], [2,2,2], [3,3,3], [4,4,4]] ) 
pprint(i)

print("i", i.shape)
print(i.ndim)
print(i.size)

j = np.array( [[[1,1],[2,2],[3,3],[4,4]] , [[1,1],[2,2],[3,3],[4,4]], [[1,1],[2,2],[3,3],[4,4]]] )
print("j", j.shape)
print(j.ndim)
print(j.size)

Accessing arrays¶

Basic indexing and slicing can be used to access array elements, as we know from lists

# a[start:stop:stride] (not inclusive of stop)
a = np.arange(8)     # another function for creating arrays
print("a", a)
print(a[2:6])
print("a[0:7:3]", a[0:7:3])
print("a[:5:2]", a[:5:2])

Negative (or modulo) indices are valid:

# To, e.g., access the last element
print(a[-1])

Exercise 2.3¶

Can you guess the output of the following cell?

print(a[2:-3:2])

For multi-dimensional arrays, tuples or index notations can be used.

# Basic indexing of a 3d array
c = np.array([[[1,2],[3,4]],[[5,6],[7,8]]])
pprint(c)

print(c.shape)
print("c[1][0][1]:", c[1][0][1]) # using index notation
print("c[1,0,1]:", c[1,0,1])     # using a tuple (more performant)

If the number of indices given is less than the number of axes, missing axes are taken as complete slices

c = np.array([[1,2,3],[5,6,8]])
pprint(c)

print(c.shape)
print("c[1][2]:", c[1][2]) # using index notation
print("c[1,2]:", c[1,2])     # using a tuple (more performant)

c = np.array([[[1,2],[3,4]],[[5,6],[7,8]]])
pprint(c)

print(c[1]) 
print(c[1,0])   
print(c[1,0,...]) # can use elipsis (3 dots) for missing indices

Slicing tuple objects can also be created to extract/manipulate indices. A slicing object takes the form slice(start, stop, step). Therefore, to take the slice a[2:-3:2] we could also write

my_slice = (slice(2, -3, 2))
print(a[my_slice])

Note that a negative step can also be used to invert an array. Take for example

d = np.arange(9).reshape(3,3)
print(d)

This can be inverted along both axes via

print(d[::-1, ::-1])

or equivalently

print(d[(slice(None, None, -1), slice(None, None, -1))])

Array copies¶

Simple assignment creates references or ‘shallow’ copies of arrays

a = np.array( [-2,6,2] )
print("a", a)
b = a
a[0] = 20

# b points to a, hence printing a and b will produce the same result
print(a)
print(b)

Use copy() to create a true or ‘deep’ copy

c = a.copy()
# check c really is an independent copy of a
c[2] = 42
print(a)
print(c)

Views from slices¶

A “view” is an array that refers to another array’s data (like references). You can create a view on an array by selecting a slice of an array. No data is copied when a view is created.

a = np.array([[1,2,3],[4,5,6],[7,8,9]]) 
pprint(a)

# Can assign a slice to a variable and change the array referred to 
# by the slice
s = a[2:3, 1:3]
print(s)
s[:,:] = -2
print("s", s)
print("a", a)

Exercise 2.4¶

Change all the elements with values 7,8,10,11 in matrix m below (i.e. bottom right corner elements) to 1000 using a slice

m = np.array([[0,1,2],[3,4,5],[6,7,8],[9,10,11]])
pprint(m)

# solution...

Reshaping arrays¶

The shape of an array can be modified, and/or its size changed:

a = np.arange(6)
print("a = ")
pprint(a)
print("a.shape = ", a.shape, "\n")

# modifying the shape attribute (not a copy) requires
# that the size remains the same
a.shape = (3,2)
pprint(a)

# Or can alter the size and shape of the array with resize().
# https://numpy.org/doc/stable/reference/generated/numpy.resize.html
mat = np.arange(6)
print("mat = ", mat)

mat1 = np.resize(mat, (3, 2))
print("mat1 = ", mat1)

mat2 = np.resize(mat, (3, 9))
print("mat2 = ", mat2)

base can be used to check if arrays share the same data (i.e. they are not copies):

mat1.base is mat

Exercise 2.5¶

What does the reshape() method do?

Vectorization and operations on arrays¶

Vectorization is why numpy arrays are great. It allows element-wise operations (avoids writing loops!). What output do the following cells give?

a = np.arange(10).reshape([2,5])
b = np.arange(10).reshape([2,5])
print(a)
print(b)

# Try these
-0.1*a
print(-0.1*a)

print(a+b)

# This is NOT matrix multiplication!
print(b)
a*b

# Use dot product for vector/matrix multiplication
# Note: .T gives you transpose of matrix i.e. reshapes it
a.dot(b.T)

Be careful: the type of data elements matters here!

Vectorization also works with functions!

def f(x):
    return x**3

x = np.array([1,2,3,4,5,6,7,8,9])
y = f(x)

print(y)

Exercise 2.6¶

Let $A$ be the two-dimensional array $$ \mathbf{A} = \left\lbrack\begin{array}{ccc} 2 & 1 & 0\cr -1 & 2 & 1\cr 0 & -1 & 2 \end{array}\right\rbrack $$

A = np.array([[2.,1.,0.], [-1.,2.,1.], [0.,-1.,2.]])
pprint(A)

Implement and apply the function $$ f(x) = x^3 + xe^x + 1 $$ to $A$. Check that you get what you expect.

# Solution...

Manipulating arrays¶

There are many methods for manipulating arrays (reshaping, joining, splitting, inserting, …). Check the documentation.

E.g.,

concatenate((a1,a2),axis=0)
split(a, indices_or_sections, axis=0)
flatten
ravel(a)
stack(arrays[, axis])
tile(a, reps)
repeat(a, repeats[, axis])
unique(ar[, return_index, return_inverse, ...])
trim_zeros(filt[, trim])
fill(scalar)
xv, yv = meshgrid(x,y)

Exercise 2.7¶

See what arrays you can create from some of the functions listed above.

Fancy indexing¶

Advanced or fancy indexing lets you do more than simple indexing.

p = np.array([[0, 1, 2],[3, 4, 5],[6, 7,8],[9, 10, 11]]) 
pprint(p)

rows = [0,0,3,3]   # indices for rows
cols = [0,2,0,2]   # indices for columns
q=p[rows,cols]
print(q)

Fancy indexing returns a copy (not a view like slicing)

# ... check if a is a view or a copy
q[0]=1000
print(q)
print(p)

Exercise 2.8¶

Use base to check whether q is a copy. Do the same for a simple indexed slice of p (e.g. p[1:2,3:4])

Logical expressions and boolean ‘masks’ can be used to find indices of elements of interest e.g.

# Find indices of elements with value less than zero
m = np.array( [[0,-1,4,20,99],[-3,-5,6,7,-10]] )
print(m)
print(m[ m < 0 ])

Exercise 2.9¶

Challenge yourself to anticipate the output of the following code before running it.

a = np.arange(10)
print(a)
mask = np.ones(len(a), dtype=bool)
mask[[0,2,4]] = False  # set certain mask values to False
result = a[mask]
print(result)

Random number generation¶

Numpy provides utilities for random number generation

# Create an array of 10 random real numbers
a = np.random.ranf(10)
print(a)

# Create a 2d array (5x5) reshaped matrix from a 1d array of (25) 
# random ints between 0 and 5 (not inclusize)
a = np.random.randint(0,high=5,size=25).reshape(5,5)   
pprint(a)

# Generate sample from normal distribution
# (mean=0, standard deviation=1)
s = np.random.standard_normal((5,5))
pprint(s)

Exercise 2.10¶

Explore other ways of generating random numbers. What other distributions can you sample?

File operations¶

Numpy provides an easy way to save data to text file and to read structured data.

# Generate an array of 5 random real numbers
pts = 5
x = np.arange(pts)
y = np.random.random(pts)
print(x)
print(y)

# data format specifiers: d = int, f = float, e = exponential
np.savetxt('savedata.txt', np.stack((x,y),axis=1), header='DATA', \
           footer='END', fmt='%d %1.4f')
!cat savedata.txt

# Reload data to an array
p = np.loadtxt('savedata.txt')
print(p)

More flexibility is offered with genfromtext() (query ?np.genfromtext)

p = np.genfromtxt('savedata.txt', skip_header=2, skip_footer=1)
print(p)

However, be aware that using plain text files for data is considered bad practice due to the lack of any form of standardization. For example, next week you will be introduced to Pandas which is designed to handle structured data.

Exercise 2.11¶

What do numpy.save() and numpy.load() do ?

Linear algebra with numpy.linalg¶

Numpy provides some linear algebra capabilities, from matrix-vector product to matrix inversion and system solution

Simple matrix vector product within numpy core:

A = np.array([[1,2,3],[4,5,6],[7,8,8]])
b = np.array([1,2,1])

print(np.dot(A,b))

The numpy.linalg module is necessary for more complex operations.

import numpy.linalg as la

Usual quantities can be computed from arrays:

n = la.norm(b)
print(n)

n = la.norm(A)
print(n)

d = la.det(A)
print(d)

And it is possible to solve linear systems, using low level C/Fortran code:

la.solve(A,b)

or to invert matrices (which is generally not a good thing to do due to computational cost!)

A_inv = la.inv(A)
print(A_inv)

The eigen decomposition (of a square matrix) can also be computed:

eival, eivec = la.eig(A)
print(eival)
pprint(eivec)

Performance¶

Python has a convenient timing function called timeit.

Can use this to measure the execution time of small code snippets.

From python: import timeit and supply code snippet as a string
From ipython: can use magic command %timeit

By default, %timeit loops (repeats) over your code 3 times and outputs the best time. It also tells you how many iterations it ran the code per loop. You can specify the number of loops and the number of iterations per loop.

%timeit -n <iterations> -r <repeats>  <code_snippet>

See

%timeit? for more information
https://docs.python.org/2/library/timeit.html

Exercise 2.12¶

Here are some timeit experiments for you to try. Which methods are faster?

# Accessing a 2d array
nd = np.arange(100).reshape((10,10))

# accessing element of 2d array
%timeit -n 10000000 -r 3 nd[5][5]
%timeit -n 10000000 -r 3 nd[5,5]

# Multiplying two vectors
x = np.arange(10E7)
%timeit -n 1 -r 10 x*x
%timeit -n 1 -r 10 x**2

# Comparing range functions and iterating in loops
# Note: some distributions may see overflow in range() example

size = int(1E6)

%timeit for x in range(size): x**2

%timeit for x in np.arange(size): x**2

%timeit np.arange(size)**2

# Extra : from the linear algebra package
%timeit -n 1 -r 10 np.dot(x,x)

Exercise 2.13¶

The Taylor series expansion for the trigonometric function $\arctan(y)$ is : $$\arctan ( y) \, = \,y - \frac{y^3}{3} + \frac{y^5}{5} - \frac{y^7}{7} + \dots $$ Now, $\arctan(1) = \pi/4 $, so $$ \pi = 4 \, \Big( 1- \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + ... \Big) $$ We can represent the series expansion using a numpy `Polynomial`, with coefficients: $$0, +1, 0, -1/3, 0, +1/5, 0, -1/7,\ldots\,$$ and use it to approximate $\pi$.

Write some code that calculates an approximation of $\pi$ using polynomials.

Bonus: Find the size of the expansion that results in an approximation of order $10^{-4}$. (You can check the difference with numpy.pi)

# Solution...

Exercise 2.14 : Darts (calculating $\pi$ again)¶

A Monte Carlo method (aka “throwing darts”)¶

Geometry gives us an expression for $\pi$:

if $N_{in}$ is the number of darts falling on the board, and $N_{tot}$ is the total number of trials

\[ \pi \approx 4 N_{in} / N_{tot} \]

Try using numpy arrays to compute the following:

Choose a sample size ntot
Generate an array of random $x$ coordinates $0 \leq x < 1$.
Generate an array of random $y$ coordinates $0 \leq y < 1$.
Count the number for which $x^2 + y^2 < 1$
Compute appromination to $\pi$
Repeat for several values of ntot and print the error.

# Solution...

SciPy library:¶

It provides a wide range of user-friendly routines (often built on NumPy) needed in scientific work. It is organised in several sub-modules:

scipy.constants
scipy.special
scipy.io
scipy.linalg
scipy.sparse
scipy.integrate
scipy.interpolate
scipy.stats
scipy.cluster
scipy.odr
scipy.optimize
scipy.signal
scipy.fftpack
scipy.ndimage
scipy.spatial

Linear algebra¶

scipy.linalg contains additional routines compared to numpy.linalg

import scipy.linalg as sla

Construct an orthonormal basis for the range of A using its SVD:

A = np.array([[2,1,1],[1,2,1],[1,1,2]])
b = np.array([1,2,3])

orth_basis = sla.orth(A)
pprint(orth_basis)

Compute the Schur decomposition of a matrix:

A = np.array([[0, 2, 2], [0, 1, 2], [1, 0, 1]])
T, Z = sla.schur(A)
pprint(T)
pprint(Z)
print("A == Z.T.Z^t: ", np.allclose(A, Z.dot(T).dot(Z.transpose())))

Integration¶

scipy.integrate contains routines to intgrate expressions as well as discrete data:

import scipy.integrate as sint

x2 = lambda x: x**2
sint.quad(x2, 0, 4)

x = np.arange(0, 10)
y = np.arange(0, 10)
y = np.power(x, 1)
sint.simps(y, x)

Interpolation¶

scipy.interpolate provides functions to interpolate data (regression, curve fitting, etc)

import scipy.interpolate as sitp

Find the Lagrange polynomial (polynomial that passes through all data points) :

%matplotlib notebook
import matplotlib.pyplot as plt

# Invent some raw data 
x=np.array([0.5,2.0,4.0,5.0,7.0,9.0])
y=np.array([0.5,0.4,0.3,0.1,0.9,0.8])

lp=sitp.lagrange(x, y)

xx= np.linspace(0.4, 9.1, 100)
fig = plt.figure()
plt.plot(x,y, "o", xx, lp(xx))
plt.show()

Fit a polynomial of arbitrary degree to a data set. Notice that the object returned is not a function this time and needs to be transformed into a function.

poly_coeffs=np.polyfit(x, y, 3)
p3 = np.poly1d(poly_coeffs)
fig = plt.figure()
plt.plot(x, y, "o", xx, p3(xx), 'g', label='Cubic')
plt.show()

Optimisation¶

scipy.optimise provides a very wide range of optimisation methods: Newton-like minimization algorithms, least square methods, root finding…

NB. Examples below are inspired from the official documentation. Interesting examples require setting up complex problems ands exceed the scope of this presentation.

import scipy.optimize as sopt

From a simple 1D problem:

def f(x):
    return -np.exp(-(x - 0.2)**4)
 
result = sopt.minimize_scalar(f)

x_min = result.x
print(x_min)

To a more complex one:

def f(x):   # rosenbrock function
    return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2

def jacobian(x):
    return np.array((-2*.5*(1 - x[0]) - 4*x[0]*(x[1] - x[0]**2), 2*(x[1] - x[0]**2)))

sopt.minimize(f, [2,-1], method="Newton-CG", jac=jacobian) 

Special functions¶

scipy.specials provides a certain number of useful mathematical functions: Airy functions, Bessel functions, elliptic integrals, gamma function, erf, binomial law, etc.

Plotting with Matplotlib¶

From the official documentation:

Matplotlib is a Python 2D plotting library which produces publication quality figures in a variety of hardcopy formats and interactive environments across platforms. […]

Matplotlib tries to make easy things easy and hard things possible. You can generate plots, histograms, power spectra, bar charts, errorcharts, scatterplots, etc., with just a few lines of code. For examples, see the sample plots and thumbnail gallery.

For simple plotting the pyplot module provides a MATLAB-like interface, particularly when combined with IPython. For the power user, you have full control of line styles, font properties, axes properties, etc, via an object oriented interface or via a set of functions familiar to MATLAB users.

In this notebook, we are mostly focused on the matplotlib.pyplot submodule

Importing Matplotlib¶

The most standard:

import matplotlib.pyplot as plt

Be careful to anonymous imports (from matplotlib.pyplot import * that will pollute your namespace.

iPython magic can be used to display plots in the browser:

%matplotlib inline

or for fancier dynamic plots

%matplotlib notebook

You can see what magic commands do by typing ?%matplotlib inline. Note that, in that case, you need to import pyplot after the magic command.

Basic plots (re: lecture 3)¶

Let’s prepare some data:

import numpy as np
%matplotlib notebook
import matplotlib.pyplot as plt

x = np.linspace(0, 2*np.pi, 100)
y = np.sin(x)
z = np.sin(x+np.pi/2)

Plotting this data is as simple as:

fig = plt.figure()

plt.plot(x, y)
plt.plot(x, z)
plt.show()

or

fig = plt.figure()

plt.plot(x,y, x, z)
plt.show()

We can specify axes limits and titles, and a global plot title:

fig = plt.figure()

plt.xlim((x.max()*0.25, x.max()*0.75))
plt.xlabel("x")
plt.title("Trigonometry")
plt.plot(x,y, x, z)
plt.show()

as well as legends:

fig = plt.figure()

plt.xlabel("x")
plt.plot(x,y, x, z)
plt.legend(("sin(x)", "sin(x+pi/2)"))
plt.show()

Line style¶

Each curve of the plot can be decorated individually:

x1 = np.linspace(0, 2*np.pi, 10)
y1 = np.sin(x1)
z1 = np.sin(x1+np.pi/2)

fig = plt.figure()

plt.plot(x1, y1, '+-r', label="sin", linewidth=5.0)
plt.plot(x1, z1, 'o--g', label="$\sin(x+\pi/2)$", markersize=7.0)
plt.legend()
plt.show()

Plots can even be annotated !

fig = plt.figure()

plt.plot(x, z)
atext = 'This is a beautiful curve'
arrowtiploc = (1.5, 0.5)
textloc=(3, 0.75)
plt.annotate(atext, xy=arrowtiploc, xytext=textloc,
            arrowprops=dict(width=5,headwidth=15,
            facecolor='red'))
plt.show()

Histograms¶

Other kinds of plots can be made, like histograms. A more complete list of plots can be found here

fig = plt.figure()

#10,000 Uniform random numbers
x2 = np.random.random(10000)
#10,000 Normally distributed random numbers
y2 = np.random.randn(10000)
#Plot both on a histogram with 50 bins
plt.hist(y2, 50)
plt.hist(x2, 50)
plt.show() #Or plt.save g("out.png")

Saving images¶

Sometimes, you might want to save images for a report or a presentation. Instead of taking a screenshot of the plot, you can save it to a file. The file format is determined by the filename extension you supply. The following common formats are supported: .png, .jpg, .pdf and .ps. There are other options to control the size and resolution of the output.

fig = plt.figure()

plt.plot(x, y, 'bx')
plt.savefig("sin.png", dpi=200)
plt.show()

Subplots¶

Sometimes, you will want to group several plots in one figure. Matplotlib provides an flexible way of managing subplots.

    (fig, axes) = plt.subplots(nrows, ncols)

subplots() returns a tuple with a reference to the figure (fig) and an array of references to each subplot axes (axes)

fig, axes = plt.subplots(2,1)
axes[0].plot(x, y, lw=3., c='r')
axes[1].plot(x, z, '--b')
axes[1].set_xlabel("x", fontsize=16)
axes[0].set_ylabel("sin(x)", fontsize=16)
axes[1].set_ylabel("sin(x+pi/2)", fontsize=16)
plt.show()

(fig, axes) = plt.subplots(2, 2, figsize=(10,7))
axes[0,0].plot(x, y, 'r')
axes[0,1].plot(x, y, 'b')
axes[1,0].plot(x, y, 'g')
axes[1,1].plot(x, y, 'k')
plt.show()

Alternatively, you can select a subplot and make it the “current subplot”

(fig, axes) = plt.subplots(2, 2)
plt.subplots_adjust(wspace = 1., hspace=1.5)

plt.sca(axes[0,0])
plt.plot(x, y, 'r')
plt.sca(axes[0,1])
plt.plot(x, y, 'b')
plt.sca(axes[1,0])
plt.plot(x, y, 'g')
plt.sca(axes[1,1])
plt.plot(x, y, 'k')
plt.show()

For more complex subplot layouts plt.subplot2grid can be used. It maps the figure to a grid, and you can define subplots spanning across several cells of the grid. It takes 4 arguments:

subplot2grid(shape,
             location,
             rowspan = 1,
             colspan = 1)

shape is the shape of the grid (number of cells), it needs to be the same for each call of the function within one figure. location is the location of the top-left corner of the subplot. rowspan and colspan define respectively the length and heigth of the plot.

fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(5)

ax1 = plt.subplot2grid((3, 3), (0, 0))
ax2 = plt.subplot2grid((3, 3), (0, 1), colspan=1)
ax3 = plt.subplot2grid((3, 3), (1, 0), colspan=2, rowspan=2)
ax4 = plt.subplot2grid((3, 3), (1, 2), rowspan=2)
# what happens if you uncomment the following line? 
#ax5 = plt.subplot2grid((3, 3), (1, 1), colspan=1, rowspan=1)
ax1.plot(x, y, 'r')
ax2.plot(x, y, 'b')
ax3.plot(x, y, 'g') 
ax4.plot(x, y, 'k')
plt.show()

2D¶

There are several ways to represent multidimensional data. One can choose to represent the data contained in a 2D domain (a matrix) with colors (in a 2D plot), or with surfaces (in a 3D plot).

fig = plt.figure()

N = 100
x = np.linspace(0,2*np.pi,N)
y = np.sin(x); z = np.cos(x)
#Create 2D field from outer product of previous 1D functions
noise = np.random.random(N**2)
u = np.outer(y,z) + noise.reshape(N,N)
plt.contourf(u, 40, cmap=plt.cm.RdYlBu_r)
plt.colorbar()
plt.show()

Exercise 2.15¶

Consider the following functions defined on a 2D domain $\Omega = [0,1]$: $$ a(x) = \tanh(50\sin(\frac{25}{2}\pi x))$$ and $$ f(x) = a(a(|x-x_0|)) + a(a(|x-x_1|)) $$ for any $x \in \Omega$, with $x_0=(0,0)$ and $x_0=(1,1)$

Plot $f$ on $\Omega$. You may want to use numpy.meshgrid.

# Solution...

3D¶

A wide variety of options are available for 3D plots, that are usually contained in specific submodules.

from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.axes3d import get_test_data


fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(15)

ax = fig.add_subplot(1,3,1, projection='3d')

# plot a curve
x = np.linspace(0.,2*np.pi,100)
ax.plot(x, np.cos(x), np.sin(x))

ax = fig.add_subplot(1, 3, 2, projection='3d')

# plot a 3D surface
X = np.arange(-5, 5, 0.25)
Y = np.arange(-5, 5, 0.25)
X, Y = np.meshgrid(X, Y)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R)
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=plt.cm.coolwarm,
                       linewidth=0, antialiased=False)
ax.set_zlim(-1.01, 1.01)
fig.colorbar(surf, shrink=0.5, aspect=10)

ax = fig.add_subplot(1, 3, 3, projection='3d')

# plot a 3D wireframe like in the example mplot3d/wire3d_demo
X, Y, Z = get_test_data(0.05)
ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)

plt.show()

Animations¶

Animations (i.e. series of images) can be created with matplotlib. However, exercise some restraint as they can require significant resources.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation

fig = plt.figure()

def get_field(a, N = 100):
    x = a*np.linspace(0,2*np.pi,N)
    y = np.sin(x); z = np.cos(x)
    out =  np.outer(y,z)
    plt.contourf(out, 40, cmap=plt.cm.RdYlBu_r)

        
ani= matplotlib.animation.FuncAnimation(plt.gcf(), get_field, frames=range(0,10),
                                       interval=500, repeat=False)
plt.show()

import random

ysample = random.sample(range(-50, 50), 100)

xdata = []
ydata = []

fig = plt.figure()

axes = plt.gca()
axes.set_xlim(0, 100)
axes.set_ylim(-50, +50)
line, = axes.plot(xdata, ydata, 'r-')

def update(i):
    xdata.append(i)
    ydata.append(ysample[i])
    line.set_xdata(xdata)
    line.set_ydata(ydata)

ani= matplotlib.animation.FuncAnimation(plt.gcf(), update, frames=100,
                                       interval=100, repeat=False)

plt.show()

Movies can be generated by saving images then calling an external program like ffmpeg, or using the writers function from the animation submodule, which also uses ffmpef under the hood.

Widgets¶

Matplotlib provides widgets, that allow you to create dynamic plots with sliders and things like that.

import matplotlib.widgets as mw

fig = plt.figure()

#Setup initial plot of sine function
x = np.linspace(0, 2*np.pi, 200)
l, = plt.plot(x, np.sin(x))
#Adjust figure to make room for slider
plt.subplots_adjust(bottom=0.15)
axslide = plt.axes([0.15, 0.05, 0.75, 0.03]) 
s = mw.Slider(axslide, 'A value', 0., 5.)
#s = widgets.FloatSlider(axslide, 'A value', 0., 5.)

#Define function
def update(A):
    l.set_ydata(np.sin(A*x))
    plt.draw()
#Bind update function to change in slider
s.on_changed(update)
plt.show()

Summary¶

NumPy: defines nd.array which is an efficient structure for large arrays, matrices and tensors, and functions to manipulate them
SciPy: defines a lot of user-friendly routines useful for scientific codes
Both rely on efficient low level code, and will generally be much faster than if you try to re-implement them yourself
Vectorization !
Matplotlib: allows you to draw and save plots, from the simplest 2D line plot to complicated 3D plots.
Refer to online documentation for a complete list of features: https://docs.scipy.org/doc/

ACSE 1 Modern Programming Methods