Tools - NumPy
NumPy is the fundamental library for scientific computing with Python. NumPy is centered around a powerful N-dimensional array object, and it also contains useful linear algebra, Fourier transform, and random number functions.
First let's make sure that this notebook works both in python 2 and 3:
In [1]:
from __future__ import division, print_function, unicode_literals
Now let's import numpy. Most people import it as np:
In [2]:
import numpy as np
The zeros function creates an array containing any number of zeros:
In [3]:
np.zeros(5)
Out[3]:
It's just as easy to create a 2D array (ie. a matrix) by providing a tuple with the desired number of rows and columns. For example, here's a 3x4 matrix:
In [4]:
np.zeros((3,4))
Out[4]:
(3, 4).
In [5]:
a = np.zeros((3,4))
a
Out[5]:
In [6]:
a.shape
Out[6]:
In [7]:
a.ndim # equal to len(a.shape)
Out[7]:
In [8]:
a.size
Out[8]:
In [9]:
np.zeros((2,3,4))
Out[9]:
In [10]:
type(np.zeros((3,4)))
Out[10]:
In [11]:
np.ones((3,4))
Out[11]:
In [12]:
np.full((3,4), np.pi)
Out[12]:
In [13]:
np.empty((2,3))
Out[13]:
In [14]:
np.array([[1,2,3,4], [10, 20, 30, 40]])
Out[14]:
In [15]:
np.arange(1, 5)
Out[15]:
It also works with floats:
In [16]:
np.arange(1.0, 5.0)
Out[16]:
Of course you can provide a step parameter:
In [17]:
np.arange(1, 5, 0.5)
Out[17]:
However, when dealing with floats, the exact number of elements in the array is not always predictible. For example, consider this:
In [18]:
print(np.arange(0, 5/3, 1/3)) # depending on floating point errors, the max value is 4/3 or 5/3.
print(np.arange(0, 5/3, 0.333333333))
print(np.arange(0, 5/3, 0.333333334))
np.linspaceFor this reason, it is generally preferable to use the linspace function instead of arange when working with floats. The linspace function returns an array containing a specific number of points evenly distributed between two values (note that the maximum value is included, contrary to arange):
In [19]:
print(np.linspace(0, 5/3, 6))
In [20]:
np.random.rand(3,4)
Out[20]:
Here's a 3x4 matrix containing random floats sampled from a univariate normal distribution (Gaussian distribution) of mean 0 and variance 1:
In [21]:
np.random.randn(3,4)
Out[21]:
To give you a feel of what these distributions look like, let's use matplotlib (see the matplotlib tutorial for more details):
In [22]:
%matplotlib inline
import matplotlib.pyplot as plt
In [23]:
plt.hist(np.random.rand(100000), normed=True, bins=100, histtype="step", color="blue", label="rand")
plt.hist(np.random.randn(100000), normed=True, bins=100, histtype="step", color="red", label="randn")
plt.axis([-2.5, 2.5, 0, 1.1])
plt.legend(loc = "upper left")
plt.title("Random distributions")
plt.xlabel("Value")
plt.ylabel("Density")
plt.show()
In [24]:
def my_function(z, y, x):
return x * y + z
np.fromfunction(my_function, (3, 2, 10))
Out[24]:
NumPy first creates three ndarrays (one per dimension), each of shape (2, 10). Each array has values equal to the coordinate along a specific axis. For example, all elements in the z array are equal to their z-coordinate:
[[[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
[[ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]
[[ 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]
[ 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]]]
So the terms x, y and z in the expression x * y + z above are in fact ndarrays (we will discuss arithmetic operations on arrays below). The point is that the function my_function is only called once, instead of once per element. This makes initialization very efficient.
In [25]:
c = np.arange(1, 5)
print(c.dtype, c)
In [26]:
c = np.arange(1.0, 5.0)
print(c.dtype, c)
Instead of letting NumPy guess what data type to use, you can set it explicitly when creating an array by setting the dtype parameter:
In [27]:
d = np.arange(1, 5, dtype=np.complex64)
print(d.dtype, d)
Available data types include int8, int16, int32, int64, uint8|16|32|64, float16|32|64 and complex64|128. Check out the documentation for the full list.
itemsizeThe itemsize attribute returns the size (in bytes) of each item:
In [28]:
e = np.arange(1, 5, dtype=np.complex64)
e.itemsize
Out[28]:
In [29]:
f = np.array([[1,2],[1000, 2000]], dtype=np.int32)
f.data
Out[29]:
In python 2, f.data is a buffer. In python 3, it is a memoryview.
In [30]:
if (hasattr(f.data, "tobytes")):
data_bytes = f.data.tobytes() # python 3
else:
data_bytes = memoryview(f.data).tobytes() # python 2
data_bytes
Out[30]:
Several ndarrays can share the same data buffer, meaning that modifying one will also modify the others. We will see an example in a minute.
In [31]:
g = np.arange(24)
print(g)
print("Rank:", g.ndim)
In [32]:
g.shape = (6, 4)
print(g)
print("Rank:", g.ndim)
In [33]:
g.shape = (2, 3, 4)
print(g)
print("Rank:", g.ndim)
In [34]:
g2 = g.reshape(4,6)
print(g2)
print("Rank:", g2.ndim)
Set item at row 1, col 2 to 999 (more about indexing below).
In [35]:
g2[1, 2] = 999
g2
Out[35]:
The corresponding element in g has been modified.
In [36]:
g
Out[36]:
In [37]:
g.ravel()
Out[37]:
In [38]:
a = np.array([14, 23, 32, 41])
b = np.array([5, 4, 3, 2])
print("a + b =", a + b)
print("a - b =", a - b)
print("a * b =", a * b)
print("a / b =", a / b)
print("a // b =", a // b)
print("a % b =", a % b)
print("a ** b =", a ** b)
Note that the multiplication is not a matrix multiplication. We will discuss matrix operations below.
The arrays must have the same shape. If they do not, NumPy will apply the broadcasting rules.
In [39]:
h = np.arange(5).reshape(1, 1, 5)
h
Out[39]:
Now let's try to add a 1D array of shape (5,) to this 3D array of shape (1,1,5). Applying the first rule of broadcasting!
In [40]:
h + [10, 20, 30, 40, 50] # same as: h + [[[10, 20, 30, 40, 50]]]
Out[40]:
In [41]:
k = np.arange(6).reshape(2, 3)
k
Out[41]:
Let's try to add a 2D array of shape (2,1) to this 2D ndarray of shape (2, 3). NumPy will apply the second rule of broadcasting:
In [42]:
k + [[100], [200]] # same as: k + [[100, 100, 100], [200, 200, 200]]
Out[42]:
Combining rules 1 & 2, we can do this:
In [43]:
k + [100, 200, 300] # after rule 1: [[100, 200, 300]], and after rule 2: [[100, 200, 300], [100, 200, 300]]
Out[43]:
And also, very simply:
In [44]:
k + 1000 # same as: k + [[1000, 1000, 1000], [1000, 1000, 1000]]
Out[44]:
In [45]:
try:
k + [33, 44]
except ValueError as e:
print(e)
Broadcasting rules are used in many NumPy operations, not just arithmetic operations, as we will see below. For more details about broadcasting, check out the documentation.
In [46]:
k1 = np.arange(0, 5, dtype=np.uint8)
print(k1.dtype, k1)
In [47]:
k2 = k1 + np.array([5, 6, 7, 8, 9], dtype=np.int8)
print(k2.dtype, k2)
Note that int16 is required to represent all possible int8 and uint8 values (from -128 to 255), even though in this case a uint8 would have sufficed.
In [48]:
k3 = k1 + 1.5
print(k3.dtype, k3)
The conditional operators also apply elementwise:
In [49]:
m = np.array([20, -5, 30, 40])
m < [15, 16, 35, 36]
Out[49]:
And using broadcasting:
In [50]:
m < 25 # equivalent to m < [25, 25, 25, 25]
Out[50]:
This is most useful in conjunction with boolean indexing (discussed below).
In [51]:
m[m < 25]
Out[51]:
In [52]:
a = np.array([[-2.5, 3.1, 7], [10, 11, 12]])
print(a)
print("mean =", a.mean())
Note that this computes the mean of all elements in the ndarray, regardless of its shape.
Here are a few more useful ndarray methods:
In [53]:
for func in (a.min, a.max, a.sum, a.prod, a.std, a.var):
print(func.__name__, "=", func())
These functions accept an optional argument axis which lets you ask for the operation to be performed on elements along the given axis. For example:
In [54]:
c=np.arange(24).reshape(2,3,4)
c
Out[54]:
In [55]:
c.sum(axis=0) # sum across matrices
Out[55]:
In [56]:
c.sum(axis=1) # sum across rows
Out[56]:
You can also sum over multiple axes:
In [57]:
c.sum(axis=(0,2)) # sum across matrices and columns
Out[57]:
In [58]:
0+1+2+3 + 12+13+14+15, 4+5+6+7 + 16+17+18+19, 8+9+10+11 + 20+21+22+23
Out[58]:
In [59]:
a = np.array([[-2.5, 3.1, 7], [10, 11, 12]])
np.square(a)
Out[59]:
Here are a few more useful unary ufuncs:
In [60]:
print("Original ndarray")
print(a)
for func in (np.abs, np.sqrt, np.exp, np.log, np.sign, np.ceil, np.modf, np.isnan, np.cos):
print("\n", func.__name__)
print(func(a))
In [61]:
a = np.array([1, -2, 3, 4])
b = np.array([2, 8, -1, 7])
np.add(a, b) # equivalent to a + b
Out[61]:
In [62]:
np.greater(a, b) # equivalent to a > b
Out[62]:
In [63]:
np.maximum(a, b)
Out[63]:
In [64]:
np.copysign(a, b)
Out[64]:
In [65]:
a = np.array([1, 5, 3, 19, 13, 7, 3])
a[3]
Out[65]:
In [66]:
a[2:5]
Out[66]:
In [67]:
a[2:-1]
Out[67]:
In [68]:
a[:2]
Out[68]:
In [69]:
a[2::2]
Out[69]:
In [70]:
a[::-1]
Out[70]:
Of course, you can modify elements:
In [71]:
a[3]=999
a
Out[71]:
You can also modify an ndarray slice:
In [72]:
a[2:5] = [997, 998, 999]
a
Out[72]:
In [73]:
a[2:5] = -1
a
Out[73]:
Also, you cannot grow or shrink ndarrays this way:
In [74]:
try:
a[2:5] = [1,2,3,4,5,6] # too long
except ValueError as e:
print(e)
You cannot delete elements either:
In [75]:
try:
del a[2:5]
except ValueError as e:
print(e)
Last but not least, ndarray slices are actually views on the same data buffer. This means that if you create a slice and modify it, you are actually going to modify the original ndarray as well!
In [76]:
a_slice = a[2:6]
a_slice[1] = 1000
a # the original array was modified!
Out[76]:
In [77]:
a[3] = 2000
a_slice # similarly, modifying the original array modifies the slice!
Out[77]:
If you want a copy of the data, you need to use the copy method:
In [78]:
another_slice = a[2:6].copy()
another_slice[1] = 3000
a # the original array is untouched
Out[78]:
In [79]:
a[3] = 4000
another_slice # similary, modifying the original array does not affect the slice copy
Out[79]:
In [80]:
b = np.arange(48).reshape(4, 12)
b
Out[80]:
In [81]:
b[1, 2] # row 1, col 2
Out[81]:
In [82]:
b[1, :] # row 1, all columns
Out[82]:
In [83]:
b[:, 1] # all rows, column 1
Out[83]:
Caution: note the subtle difference between these two expressions:
In [84]:
b[1, :]
Out[84]:
In [85]:
b[1:2, :]
Out[85]:
The first expression returns row 1 as a 1D array of shape (12,), while the second returns that same row as a 2D array of shape (1, 12).
In [86]:
b[(0,2), 2:5] # rows 0 and 2, columns 2 to 4 (5-1)
Out[86]:
In [87]:
b[:, (-1, 2, -1)] # all rows, columns -1 (last), 2 and -1 (again, and in this order)
Out[87]:
If you provide multiple index arrays, you get a 1D ndarray containing the values of the elements at the specified coordinates.
In [88]:
b[(-1, 2, -1, 2), (5, 9, 1, 9)] # returns a 1D array with b[-1, 5], b[2, 9], b[-1, 1] and b[2, 9] (again)
Out[88]:
In [89]:
c = b.reshape(4,2,6)
c
Out[89]:
In [90]:
c[2, 1, 4] # matrix 2, row 1, col 4
Out[90]:
In [91]:
c[2, :, 3] # matrix 2, all rows, col 3
Out[91]:
If you omit coordinates for some axes, then all elements in these axes are returned:
In [92]:
c[2, 1] # Return matrix 2, row 1, all columns. This is equivalent to c[2, 1, :]
Out[92]:
In [93]:
c[2, ...] # matrix 2, all rows, all columns. This is equivalent to c[2, :, :]
Out[93]:
In [94]:
c[2, 1, ...] # matrix 2, row 1, all columns. This is equivalent to c[2, 1, :]
Out[94]:
In [95]:
c[2, ..., 3] # matrix 2, all rows, column 3. This is equivalent to c[2, :, 3]
Out[95]:
In [96]:
c[..., 3] # all matrices, all rows, column 3. This is equivalent to c[:, :, 3]
Out[96]:
In [97]:
b = np.arange(48).reshape(4, 12)
b
Out[97]:
In [98]:
rows_on = np.array([True, False, True, False])
b[rows_on, :] # Rows 0 and 2, all columns. Equivalent to b[(0, 2), :]
Out[98]:
In [99]:
cols_on = np.array([False, True, False] * 4)
b[:, cols_on] # All rows, columns 1, 4, 7 and 10
Out[99]:
In [100]:
b[np.ix_(rows_on, cols_on)]
Out[100]:
In [101]:
np.ix_(rows_on, cols_on)
Out[101]:
If you use a boolean array that has the same shape as the ndarray, then you get in return a 1D array containing all the values that have True at their coordinate. This is generally used along with conditional operators:
In [102]:
b[b % 3 == 1]
Out[102]:
In [103]:
c = np.arange(24).reshape(2, 3, 4) # A 3D array (composed of two 3x4 matrices)
c
Out[103]:
In [104]:
for m in c:
print("Item:")
print(m)
In [105]:
for i in range(len(c)): # Note that len(c) == c.shape[0]
print("Item:")
print(c[i])
If you want to iterate on all elements in the ndarray, simply iterate over the flat attribute:
In [106]:
for i in c.flat:
print("Item:", i)
In [107]:
q1 = np.full((3,4), 1.0)
q1
Out[107]:
In [108]:
q2 = np.full((4,4), 2.0)
q2
Out[108]:
In [109]:
q3 = np.full((3,4), 3.0)
q3
Out[109]:
In [110]:
q4 = np.vstack((q1, q2, q3))
q4
Out[110]:
In [111]:
q4.shape
Out[111]:
In [112]:
q5 = np.hstack((q1, q3))
q5
Out[112]:
In [113]:
q5.shape
Out[113]:
This is possible because q1 and q3 both have 3 rows. But since q2 has 4 rows, it cannot be stacked horizontally with q1 and q3:
In [114]:
try:
q5 = np.hstack((q1, q2, q3))
except ValueError as e:
print(e)
In [115]:
q7 = np.concatenate((q1, q2, q3), axis=0) # Equivalent to vstack
q7
Out[115]:
In [116]:
q7.shape
Out[116]:
As you might guess, hstack is equivalent to calling concatenate with axis=1.
In [117]:
q8 = np.stack((q1, q3))
q8
Out[117]:
In [118]:
q8.shape
Out[118]:
In [119]:
r = np.arange(24).reshape(6,4)
r
Out[119]:
Now let's split it in three equal parts, vertically:
In [120]:
r1, r2, r3 = np.vsplit(r, 3)
r1
Out[120]:
In [121]:
r2
Out[121]:
In [122]:
r3
Out[122]:
There is also a split function which splits an array along any given axis. Calling vsplit is equivalent to calling split with axis=0. There is also an hsplit function, equivalent to calling split with axis=1:
In [123]:
r4, r5 = np.hsplit(r, 2)
r4
Out[123]:
In [124]:
r5
Out[124]:
In [125]:
t = np.arange(24).reshape(4,2,3)
t
Out[125]:
Now let's create an ndarray such that the axes 0, 1, 2 (depth, height, width) are re-ordered to 1, 2, 0 (depth→width, height→depth, width→height):
In [126]:
t1 = t.transpose((1,2,0))
t1
Out[126]:
In [127]:
t1.shape
Out[127]:
By default, transpose reverses the order of the dimensions:
In [128]:
t2 = t.transpose() # equivalent to t.transpose((2, 1, 0))
t2
Out[128]:
In [129]:
t2.shape
Out[129]:
NumPy provides a convenience function swapaxes to swap two axes. For example, let's create a new view of t with depth and height swapped:
In [130]:
t3 = t.swapaxes(0,1) # equivalent to t.transpose((1, 0, 2))
t3
Out[130]:
In [131]:
t3.shape
Out[131]:
NumPy 2D arrays can be used to represent matrices efficiently in python. We will just quickly go through some of the main matrix operations available. For more details about Linear Algebra, vectors and matrics, go through the Linear Algebra tutorial.
The T attribute is equivalent to calling transpose() when the rank is ≥2:
In [132]:
m1 = np.arange(10).reshape(2,5)
m1
Out[132]:
In [133]:
m1.T
Out[133]:
The T attribute has no effect on rank 0 (empty) or rank 1 arrays:
In [134]:
m2 = np.arange(5)
m2
Out[134]:
In [135]:
m2.T
Out[135]:
We can get the desired transposition by first reshaping the 1D array to a single-row matrix (2D):
In [136]:
m2r = m2.reshape(1,5)
m2r
Out[136]:
In [137]:
m2r.T
Out[137]:
Let's create two matrices and execute a matrix dot product using the dot method.
In [138]:
n1 = np.arange(10).reshape(2, 5)
n1
Out[138]:
In [139]:
n2 = np.arange(15).reshape(5,3)
n2
Out[139]:
In [140]:
n1.dot(n2)
Out[140]:
Caution: as mentionned previously, n1*n2 is not a dot product, it is an elementwise product.
In [141]:
import numpy.linalg as linalg
m3 = np.array([[1,2,3],[5,7,11],[21,29,31]])
m3
Out[141]:
In [142]:
linalg.inv(m3)
Out[142]:
You can also compute the pseudoinverse using pinv:
In [143]:
linalg.pinv(m3)
Out[143]:
In [144]:
m3.dot(linalg.inv(m3))
Out[144]:
You can create an identity matrix of size NxN by calling eye:
In [145]:
np.eye(3)
Out[145]:
The qr function computes the QR decomposition of a matrix:
In [146]:
q, r = linalg.qr(m3)
q
Out[146]:
In [147]:
r
Out[147]:
In [148]:
q.dot(r) # q.r equals m3
Out[148]:
The det function computes the matrix determinant:
In [149]:
linalg.det(m3) # Computes the matrix determinant
Out[149]:
The eig function computes the eigenvalues and eigenvectors of a square matrix:
In [150]:
eigenvalues, eigenvectors = linalg.eig(m3)
eigenvalues # λ
Out[150]:
In [151]:
eigenvectors # v
Out[151]:
In [152]:
m3.dot(eigenvectors) - eigenvalues * eigenvectors # m3.v - λ*v = 0
Out[152]:
The svd function takes a matrix and returns its singular value decomposition:
In [153]:
m4 = np.array([[1,0,0,0,2], [0,0,3,0,0], [0,0,0,0,0], [0,2,0,0,0]])
m4
Out[153]:
In [154]:
U, S_diag, V = linalg.svd(m4)
U
Out[154]:
In [155]:
S_diag
Out[155]:
The svd function just returns the values in the diagonal of Σ, but we want the full Σ matrix, so let's create it:
In [156]:
S = np.zeros((4, 5))
S[np.diag_indices(4)] = S_diag
S # Σ
Out[156]:
In [157]:
V
Out[157]:
In [158]:
U.dot(S).dot(V) # U.Σ.V == m4
Out[158]:
In [159]:
np.diag(m3) # the values in the diagonal of m3 (top left to bottom right)
Out[159]:
In [160]:
np.trace(m3) # equivalent to np.diag(m3).sum()
Out[160]:
The solve function solves a system of linear scalar equations, such as:
In [161]:
coeffs = np.array([[2, 6], [5, 3]])
depvars = np.array([6, -9])
solution = linalg.solve(coeffs, depvars)
solution
Out[161]:
Let's check the solution:
In [162]:
coeffs.dot(solution), depvars # yep, it's the same
Out[162]:
Looks good! Another way to check the solution:
In [163]:
np.allclose(coeffs.dot(solution), depvars)
Out[163]:
Instead of executing operations on individual array items, one at a time, your code is much more efficient if you try to stick to array operations. This is called vectorization. This way, you can benefit from NumPy's many optimizations.
For example, let's say we want to generate a 768x1024 array based on the formula $sin(xy/40.5)$. A bad option would be to do the math in python using nested loops:
In [164]:
import math
data = np.empty((768, 1024))
for y in range(768):
for x in range(1024):
data[y, x] = math.sin(x*y/40.5) # BAD! Very inefficient.
Sure, this works, but it's terribly inefficient since the loops are taking place in pure python. Let's vectorize this algorithm. First, we will use NumPy's meshgrid function which generates coordinate matrices from coordinate vectors.
In [165]:
x_coords = np.arange(0, 1024) # [0, 1, 2, ..., 1023]
y_coords = np.arange(0, 768) # [0, 1, 2, ..., 767]
X, Y = np.meshgrid(x_coords, y_coords)
X
Out[165]:
In [166]:
Y
Out[166]:
As you can see, both X and Y are 768x1024 arrays, and all values in X correspond to the horizontal coordinate, while all values in Y correspond to the the vertical coordinate.
Now we can simply compute the result using array operations:
In [167]:
data = np.sin(X*Y/40.5)
Now we can plot this data using matplotlib's imshow function (see the matplotlib tutorial).
In [168]:
import matplotlib.pyplot as plt
import matplotlib.cm as cm
fig = plt.figure(1, figsize=(7, 6))
plt.imshow(data, cmap=cm.hot, interpolation="bicubic")
plt.show()
In [169]:
a = np.random.rand(2,3)
a
Out[169]:
In [170]:
np.save("my_array", a)
Done! Since the file name contains no file extension was provided, NumPy automatically added .npy. Let's take a peek at the file content:
In [171]:
with open("my_array.npy", "rb") as f:
content = f.read()
content
Out[171]:
To load this file into a NumPy array, simply call load:
In [172]:
a_loaded = np.load("my_array.npy")
a_loaded
Out[172]:
In [173]:
np.savetxt("my_array.csv", a)
Now let's look at the file content:
In [174]:
with open("my_array.csv", "rt") as f:
print(f.read())
This is a CSV file with tabs as delimiters. You can set a different delimiter:
In [175]:
np.savetxt("my_array.csv", a, delimiter=",")
To load this file, just use loadtxt:
In [176]:
a_loaded = np.loadtxt("my_array.csv", delimiter=",")
a_loaded
Out[176]:
In [177]:
b = np.arange(24, dtype=np.uint8).reshape(2, 3, 4)
b
Out[177]:
In [178]:
np.savez("my_arrays", my_a=a, my_b=b)
Again, let's take a peek at the file content. Note that the .npz file extension was automatically added.
In [179]:
with open("my_arrays.npz", "rb") as f:
content = f.read()
repr(content)[:180] + "[...]"
Out[179]:
You then load this file like so:
In [180]:
my_arrays = np.load("my_arrays.npz")
my_arrays
Out[180]:
This is a dict-like object which loads the arrays lazily:
In [181]:
my_arrays.keys()
Out[181]:
In [182]:
my_arrays["my_a"]
Out[182]:
Now you know all the fundamentals of NumPy, but there are many more options available. The best way to learn more is to experiment with NumPy, and go through the excellent reference documentation to find more functions and features you may be interested in.