%config InlineBackend.figure_formats = ['svg']
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [5, 3]
In this lecture we will talk more about NumPy and introduce Matplotlib.
ndarray
array typeTo access the numpy
package it is standard to use the import:
import numpy as np
The first thing to note is that when working with ndarrays
: Operations are element-wise
a = np.arange(4.)
b = a + 1
print("type(a) = {}\na = {}\nb = {}".format(type(a), a, b))
c = a + b
print('c =', c)
a = np.arange(4)
b = np.linspace(1., 2.5, num=4)
c = a * b
print(f" a = {a} \n b = {b} \n c = {c}")
Note: The elements in an ndarray
has a common dtype
the results of array operations is promoted to the proper type.
print(a.dtype, b.dtype, c.dtype, sep=', ')
If shapes of operands are incompatible, NumPy will give you an error:
a = np.arange(4)
b = np.linspace(1, 2, num=3)
print(a.shape, b.shape)
c = a * b
If the shapes of arrays does not match, NumPy tries to apply broadcasting rules
.shape[d] == 1
a = np.zeros(shape=(3, 4))
b = np.arange(4.).reshape((1, 4))
print('a =\n', a)
print('b =', b)
print(a.shape, b.shape)
c = a + b # What happens here!?!
print(c)
np.newaxis
¶np.newaxis
(or None
)a = np.zeros(shape=(3, 4))
b = np.arange(4.)
c = a + b[np.newaxis, :]
print(f'a =\n{a}\nb = {b}\nc =\n{c}')
If you like to think in index notation the above operation is equivalent to: $ c_{ij} = a_{ij} + b_{i} $
What is np.newaxis
actually doing?
x = b[np.newaxis, :]
print(b.shape)
print(x.shape)
Take the difference $d_{ij}$ of all pairs of entries in two vectors $a_i$ and $b_j$ $$ d_{ij} \equiv a_i - b_j $$
Solve by applying broadcasting to both arrays
a_i = np.arange(4.)
b_j = np.arange(4., 0., -1.)
print(f'a_i = {a_i}\nb_j = {b_j}')
d_ij = ? # ToDo
In many algorithms you need to work on parts of an array, i.e. only selected entries will be used.
a = np.linspace(1., 10., 7)
index = np.array([0, 2, 3])
print( f"a = {a}" )
print( f"index = {index}" )
print()
print( a[0] ) # Single element
print( a[[0, 2, 3]] )# Indexing with index array - Python list
print( a[index] ) # Indexing with index array - NumPy array
mask = a > 4
print( f"mask = {mask}" )
print( f"a[mask] = {a[mask]}" )
print( a[0:2] ) # Slicing: Elements from (and including) no. 0 to no. 2 (not including)
print( a[::2] ) # Slicing: Every second element from start to end
print( a[:0:-1] ) # Slicing: Every element _except_ the first one, in reverse order
a = np.arange(24).reshape(3, 2, 4)
print( f"a = \n{a}" )
print()
print( "a[0,:,:] = \n{}".format(a[0,:,:]) ) # "block 0"
print()
print( "a[:,1,:] = \n{}".format(a[:,1,:]) ) # "row 1"
print()
print( "a[:,:,2] = \n{}".format(a[:,:,2]) ) # "column 2"
There are many numerical routines available in NumPy
np.sin
, np.cosh
, np.exp
, np.log
, np.log10
, np.sign
, etc.np.max
, np.min
, np.sum
, np.prod
, np.mean
, np.median
, etc.np.cumsum
, np.cumprod
np.trapz
, np.diff
, np.gradient
np.dot
, np.tensordot
, np.einsum
, etc.np.linalg
np.random
np.fft
np.polynomial
see the NumPy reference manual for a complete list
np.sum
, np.prod
)¶a = 0.1 * np.arange(1., 25.).reshape(2, 3, 4)
print(f"a =\n{a}")
Sum all elements: $\sum_{ijk} a_{ijk}$
np.sum(a)
Product of all elements: $\prod_{ijk} a_{ijk}$
np.prod(a)
To work over a subset of axes, use the axis
key-word argument
Sum all elements over 2nd axis: $\sum_{j} a_{ijk}$
np.sum(a, axis=1)
Product of all elements over the 2nd and 3rd axis: $\prod_{jk} a_{ijk}$
np.prod(a, axis=(1, 2))
np.diff
¶The difference between neighbouring elements $$y_i = x_{i+1} - x_i$$
can be computed by: y = np.diff(x, n=1, axis=-1)
n :
number of times to take the differenceaxis :
axis to work on, default -1
x = np.arange(5.)
print(x)
y = np.diff(x)
print(y)
z = np.diff(x, 2)
print(z)
np.gradient
¶Approximate the gradient calculation
$$\frac{df(x_n)}{dx} \approx \frac{f(x_{n+1})- 2f(x_n) + f(x_{n-1})}{2\Delta x}$$
use np.gradient(f , [dx, dy, ....])
Example: $$f(x) = \cos(x) \quad \Rightarrow \quad f'(x) \equiv \frac{d}{dx}f(x) = \frac{d}{dx} \cos(x) = -\sin(x)$$
dx = 0.1
x = np.arange(0, 3*np.pi, dx)
f = np.cos(x)
f_prime = np.gradient(f, dx)
import matplotlib.pyplot as plt
plt.plot(x, f, label='$f(x)$')
plt.plot(x, f_prime, label='$f\'(x)$')
plt.xlabel('$x$')
plt.legend();
np.cumsum
¶For the Riemann sum approximation of an integral $$F(x) = \int_{x_0}^x f(y) \, dy \quad \Rightarrow \quad F(x_n) \approx \Delta x \cdot \sum_{i = 0}^n f(x_i)$$
use np.cumsum(f, axis=None, out=None)
axis
is the axis to work along (default is to flatten array and use all elements)
Example: $f(x) = \cos(x) \Rightarrow F(x) = \sin(x)$
x = np.linspace(0, 4*np.pi, num=400)
dx = x[1] - x[0] # All distances equal
f = np.cos(x)
F = np.cumsum(f) * dx
plt.plot(x, f, label='$f(x)$')
plt.plot(x, F, label='$F(x)$')
plt.xlabel('$x$')
plt.legend();
For integration over an entire interval $$ \int_{x_0}^{x_{-1}} f(x) \, dx $$
numpy.trapz(f, x=None, dx=1.0, axis=-1)
(trapetzoidal rule)x
or dx
!F_trapz = np.trapz(f, x=x)
print("Trapezoidal integration:", F_trapz)
print("Rieman integration: ", F[-1])
np.dot
¶+
, -
, *
, /
etc., work element-wise on NumPy ndarray
snp.dot(a, b)
ndarrays
have .dot
method, see belowThe dot product np.dot
is defined as
a = np.array([1., 2., 3.])
b = np.array([2., 1., 0.])
np.dot(a, b)
a.dot(b)
a * b
np.dot
does matrix-vector multiplication when given one 2D and one 1D array
A = np.arange(6.).reshape( (3,2) )
b = np.arange(2.)
print( "A =\n{}".format(A) )
print( "b = {}".format(b) )
A.dot(b)
np.dot
does matrix multiplication when given two 2D arrays
A = np.arange(9).reshape((3,3))
B = np.random.randint(0,10,(3,3))
print("A =", A, 'B =', B, sep='\n')
A.dot(B)
T
view, e.g. B.T
.conjugate
or .conj
methodsnp.swapaxes
etc.B = np.arange(6).reshape(2, 3)
B
B.T
np.swapaxes(B, 0, 1)
Complex valued arrays and complex conjugate
C = (1 + 1j) * B
C
C.conj()
np.linalg
¶Provides additional operations
norm(A)
: matrix or vector norm (default: 2-norm)det(A)
: matrix determinantx = solve(A, b)
: linear equation system solver: $A \cdot \mathbf{x} = \mathbf{b}$inv(A)
, pinv(A)
: matrix inverse and pseudo inverseQ, R = qr(A)
QR-factorizationeig(A)
: matrix eigen value and eigen vector decompositionsvd(A)
: matrix single-value decompositionExample: inverse
A = np.array([
[1, 0, 1],
[3, 1, 0],
[1, 2, 1]])
A_inv = np.linalg.inv(A)
print(A_inv)
A.dot(A_inv)
For more information use the built-in help!
np.matrix
¶A @ B
Here are the previous examples with np.dot
replaced by @
A = np.arange(6.).reshape((3, 2))
b = np.arange(2.)
A @ b
A = np.arange(9).reshape((3 ,3))
B = np.random.randint(0, 10, size=(3, 3))
A @ B
A = np.array([
[1, 0, 1],
[3, 1, 0],
[1, 2, 1]])
A_inv = np.linalg.inv(A)
A @ A_inv
There are many additional sub-packages and functions in NumPy not covered here, for example the modules
np.fft
,np.random
, andnp.polynomial
,for more information use the built-in help()
and the online documentation.
Visualization is a cornerstone of scientific computing
There are many tools for data visualisation, e.g.
For Python the dominating package is:
Here we will only look into Matplotlib.
Matplotlib can be used in
Matplotlib tries to make easy things easy and hard things possible.
You can generate
with just a few lines of code.
For a sampling, see the screenshots, thumbnail gallery, and examples directory.
matplotlib.pyplot
¶matplotlib.pyplot
is the "main" module for plottingimport matplotlib.pyplot as plt
To understand how to work with Matplotlib we need to know how Matplotlib constructs its figures.
The "basic" object for all figures is the figure
object.
The figure object contains all the axes, plots and axis lables, etc. that you create
The figure object can be explicitly created using plt.figure()
,
However, it is created automatically when using most plotting-commands.
import matplotlib.pyplot as plt
plt.figure()
plt.plot([0, 2, 3], [5, 3, 1]);
plt.show()
The call to plt.plot(...)
creates a line plot (it also automatically creates an axes object and makes it the current axis).
A figure comprise a number of components
(Source: the Matplotlib homepage)
Methods for the figure object, or for individual axes objects include:
plt.figure()
plt.title('Bell progress')
plt.plot([0,2,3], [5,3,1], label='Tim')
plt.plot([0,1,3], [2,3,5], label='Alice')
plt.legend(title='User', loc='best')
plt.xlabel('time')
plt.ylabel('bells')
plt.xlim(right=3.5)
plt.ylim([0., 5.5])
plt.grid(True)
plt.text( 2.5, 3., r'$\Omega = \sum_i \frac{K_i}{2 \cdot B_i}$', size=15 )
plt.yticks(np.arange(6), ['', 'One', '**', 3, r'$\Delta$', '>>>>>'])
plt.tick_params(axis='x', direction='in', color='r', labelcolor='g',
width=2, length=8, top=False)
plt.savefig('my_plot.svg')
plt.show()
r"..."
) with LaTeX-syntax to display symbols etc. (See the simple example above.)You can combine plots with completely different scales in one direction into one plot (or subplot)
by using plt.twinx()
or plt.twiny()
.
plt.figure()
plt.plot([1, 2, 3, 4], [1, 2.5, 3.1, 4], label='Data 1')
plt.ylabel('Data 1 [m]')
plt.legend()
plt.twinx()
plt.scatter( [1, 2, 3, 4], [60, 48, 42, 30], label='Data 2' )
plt.ylabel('Data 2 [kg]')
plt.legend()
plt.xlabel('X')
plt.show()
Colors, markers and linestyles can be specified using a Matlab-like notation for simple cases
"r*-"
to give a red line with *
-markers at the points.color=
keyword-argumentThere are four color specifications:
#0000FF
for blue etc.(r, g, b, alpha)
with values in $a,b,c,\alpha \in [0, 1]$. The last parameter $\alpha$ is optional and controls the transparency."0.0"
and for white "1.0"
(Note: Transparency can also be set by the alpha=
key-word argument)
plt.plot([0, 1, 2], [1.2, 1.5, 1.0], color='OliveDrab')
plt.plot([0, 1, 2], [0.8, 0.4, 0.2], color=(0.2, 0.8, 0.8, 0.5))
plt.plot([0, 1, 2], [0.5, 0.8, 1.3], color="0.6", alpha=0.5, linewidth=5)
plt.show()
-
: solid lines--
: dashed lines-.
: dot-dashed linesdashes=[dash_length, space_length, ...]
: arbitrary dashesx = np.linspace(0, 1, num=10)
a, b, c, d = np.cumsum(np.random.random((4, len(x))), axis=-1)
plt.plot(x, a, '-', label='a')
plt.plot(x, b, '--', label='b')
plt.plot(x, c, '-.', label='c')
plt.plot(x, d, dashes=[6, 1, 1, 1, 1, 1], label='d')
plt.legend(loc='best'); plt.show()
Markers are set in the fmt
argument or the marker=
keyword argument.
marker | description | marker | description | marker | description | marker | description | |||
---|---|---|---|---|---|---|---|---|---|---|
"." | point | "+" | plus | "," | pixel | "x" | cross | |||
"o" | circle | "D" | diamond | "d" | thin_diamond | |||||
"8" | octagon | "s" | square | "p" | pentagon | "*" | star | |||
"|" | vertical line | "_" | horizontal line | "h" | hexagon1 | "H" | hexagon2 | |||
0 | tickleft | 4 | caretleft | "<" | triangle_left | "3" | tri_left | |||
1 | tickright | 5 | caretright | ">" | triangle_right | "4" | tri_right | |||
2 | tickup | 6 | caretup | "^" | triangle_up | "2" | tri_up | |||
3 | tickdown | 7 | caretdown | "v" | triangle_down | "1" | tri_down | |||
"None" | nothing | None |
nothing | " " | nothing | "" | nothing |
It is also possible to construct new markers by giving a list of xy-coordinate tuples (relative to the center of the marker at (0,0)
).
plt.plot([0, 1, 2], [1.2, 1.5, 1.0], 'kD-', markersize=10 )
plt.plot([0, 1, 2], [0.5, 0.8, 1.3], color="green", marker="8", markersize=15 )
plt.plot([0, 1, 2], [0.8, 0.4, 0.2], color='SlateBlue', alpha=0.75,
marker=[(-3, -5), (0, 3), (3, -5), (0, -3)], markersize=25)
plt.show()
There is a wide range of plot functions available in Matplotlib.
plt.plot
and plt.scatter
plt.plot
¶The plt.plot()
is very flexible
plt.plot(x, y)
will plot the values in x
and y
. x
, y
can be single values or iterable items (list, tuples, NumPy arrays etc.) of coordinatesplot(x1, y1, x2, y2, ...)
will plot the sets of x
and y
:s at once in the same figureplot(y)
will plot y with x as integers starting from 0Common keyword arguments to plt.plot
are
color=
, alpha=
,linestyle=
, linewidth=
, fillstyle=
,marker=
, markersize=
,label=
(will be shown by plt.legend
)Note: It is still possible to use Matlab-like colour/marker/linestyle specifications like "rx-."
Another example:
coord = np.array( [[0.05, 0.05], [0.4, 0.3], [0.6, 0.5], [0.7, 0.55], [0.95, 0.61]] )
print(coord)
plt.plot(coord[:, 0], coord[:, 1], color='g', label='data', fillstyle='bottom', marker='D', markersize=10 )
plt.plot(coord[:, 1], coord[:, 0], color='b', label='flipped data', linestyle='--', linewidth=4 )
plt.legend(loc='lower right')
plt.show()
Multi-dimensional y
array example
x = np.linspace(0, 1, num=10)
y = np.cumsum(np.random.random((len(x), 4)), axis=0)
print(x.shape, y.shape)
plt.plot(x, y)
plt.show()
plt.semilogx
, plt.semilogy
, plt.loglog
¶Any axis can be scaled logarithmically by using
plt.semilogx
: scale the x-axis logaritmicallyplt.semilogy
: the y-axisplt.loglog
: scale both axesx = np.linspace(0, 100, num=1000)
plt.loglog(x, x**2, label=r'$x^2$')
plt.loglog(x, x**3, label=r'$x^3$')
plt.legend(loc='upper left'); plt.show()
plt.semilogy(x, x**2, label=r'$x^2$')
plt.semilogy(x, x**3, label=r'$x^3$')
plt.legend(loc='lower right'); plt.show()
plt.errorbar
¶For errorbars (or any other interval you want to show like an errorbar), use
plt.errorbar(x, y, xerr=..., yerr=...)
x = np.linspace(0, 1, num=10)
y = x*(x-1)*(x+1)
xerror = np.random.normal(size=len(x), scale=0.1)
yerr = np.random.normal(size=len(x), scale=0.05)
plt.errorbar(x, -y, xerr=xerror, yerr=yerr, ecolor='r');
plt.errorbar(x, y, xerr=yerr, ecolor='g');
plt.fill_between
¶plt.fill_between(x, y1, y2=)
x = np.linspace(1e-10, 4*np.pi, num=1000)
y1 = np.sin(x + np.pi*0.5)
y2 = np.sin(x) / x
plt.fill_between(x, y1, y2=y2, edgecolor='Purple', alpha=0.5, linewidth=4);
plt.fill_between(x, y2, edgecolor='r', facecolor='SeaGreen', alpha=0.75);
plt.hist
¶plt.hist(y, bins=10, normed=False, histtype= ...)
x = 100 + 15 * np.random.randn(100000)
plt.hist(x, bins=40, density=True, facecolor='g', alpha=0.75);
plt.hist(x, bins=20, histtype='stepfilled', facecolor='PowderBlue', alpha=0.5);
plt.contour
, plt.contourf
¶You can create filled and non-filled contour plots using plt.contourf
and plt.contour
respectively.
cmap=
keyword argumentx = np.linspace(-4*np.pi, 4*np.pi)
X, Y = np.meshgrid(x, x)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R) / R
plt.contour(X, Y, Z, 15, cmap=plt.cm.RdBu)
plt.axis('image');
plt.contourf(X, Y, Z, 15, cmap=plt.cm.summer); plt.axis('image');
plt.axis
, plt.subplot
¶axis
object in Matplotlibplt.plot
x = np.linspace(-1, 1, num=10)
fig = plt.figure()
ax = fig.add_subplot()
ax.plot(x, x**2)
ax.set_xlabel('x');
figure
and axis
objectsplt.subplot
¶ax = plt.subplot(nrows, ncols, plot_number)
: call for all values of plot_number
fig1, (ax1, ax2, ...) = plt.subplots(nrows, ncols)
: call onceBoth these gives a nrows x ncols
grid of subplots.
plt.subplot(1, 2, 1)
plt.plot([0,2,3], [5,3,1])
plt.subplot(1, 2, 2)
plt.plot([0,1,3], [2,4,3])
plt.show()
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot([0,2,3], [5,3,1])
ax2.plot([0,1,3], [2,4,3])
plt.show()
How do we make a figure with
Solution:
x, y_L, y_R = [0,2,3], [5,3,1], [2,4,3]
plt.subplot(1, 2, 1)
plt.plot(x, y_L)
plt.subplot(2, 2, 2)
plt.plot(x, y_R)
plt.subplot(2, 2, 4)
plt.plot(x, y_R)
plt.tight_layout()
Tip: For more advanced subplot layouts use the GridSpec
class
from mpl_toolkits.mplot3d import Axes3D
projection='3d'
keyword argument when creating axes gives a 3D plot%config InlineBackend.figure_formats = ['retina']
theta = np.linspace(-4, 4, num=100)
z = 0.5*theta; r = z**2 + 1
x = r * np.sin(np.pi*theta)
y = r * np.cos(np.pi*theta)
from mpl_toolkits.mplot3d import Axes3D
plt.figure(figsize=(12, 12))
plt.subplot(1,1,1, projection='3d')
plt.plot(x, y, z, label='parametric curve')
plt.legend()
plt.show()
x = np.linspace(-4*np.pi, 4*np.pi, num=40)
X, Y = np.meshgrid(x, x)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R) / R
plt.figure(figsize=(12, 12))
ax = plt.subplot(1,1,1, projection='3d')
ax.plot_wireframe(X, Y, Z, alpha=0.5)
plt.show()
x = np.linspace(-4*np.pi, 4*np.pi, num=100)
X, Y = np.meshgrid(x, x)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R) / R
plt.figure(figsize=(12, 12))
ax = plt.subplot(1,1,1, projection='3d')
ax.plot_surface(X, Y, Z, cmap=plt.cm.rainbow, rstride=1, cstride=1, alpha=.5, linewidth=0)
ax.contour(X, Y, Z, zdir='y', offset=15)
ax.contour(X, Y, Z, zdir='x', offset=-15)
plt.show()
If you want to plot a large number of elements it is somewhat inefficient to plot them one by one using the methods shown above.
To remedy this Matplotlib have a collections submodule that helps in collecting many objects of the same type that can be used efficiently when drawing images.
There are collection-types for polygons, lines, triangular meshes, paths and more.
For example the matplotlib.collections LineCollection(segments, ...) gathers segments of lines for convinient plotting of many line-segments.
Ny = 10
x = np.arange(8)
# Ny lines where y is randoms between 0 and 2 * Ny + 3
from random import random
lines = []
for i in np.arange(Ny):
line = []
for xj in x:
y = random() * 3 + 2 * i
line.append( (xj, y) )
#print(line)
lines.append(line)
from matplotlib.collections import LineCollection
line_segments = LineCollection(lines)
fig = plt.figure(figsize=(15, 8))
ax = fig.gca()
ax.add_collection(line_segments)
ax.set_ylim((0, 2*Ny+3))
ax.set_xlim((0, np.amax(x)))
plt.show()
Have a look at the thumbnail gallery and see if you see what you want to do.
import this