Jupyter Snippet NP ch02-code-listing
Jupyter Snippet NP ch02-code-listing
Chapter 2: Vectors, matrices and multidimensional arrays
Robert Johansson
Source code listings for Numerical Python - Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib (ISBN 978-1-484242-45-2).
import numpy as np
The NumPy array object
data = np.array([[1, 2], [3, 4], [5, 6]])
type(data)
numpy.ndarray
data
array([[1, 2],
[3, 4],
[5, 6]])
data.ndim
2
data.shape
(3, 2)
data.size
6
data.dtype
dtype('int64')
data.nbytes
48
Data types
np.array([1, 2, 3], dtype=np.int)
array([1, 2, 3])
np.array([1, 2, 3], dtype=np.float)
array([ 1., 2., 3.])
np.array([1, 2, 3], dtype=np.complex)
array([ 1.+0.j, 2.+0.j, 3.+0.j])
data = np.array([1, 2, 3], dtype=np.float)
data
array([ 1., 2., 3.])
data.dtype
dtype('float64')
data = np.array([1, 2, 3], dtype=np.int)
data.dtype
dtype('int64')
data
array([1, 2, 3])
data = np.array([1, 2, 3], dtype=np.float)
data
array([ 1., 2., 3.])
data.astype(np.int)
array([1, 2, 3])
d1 = np.array([1, 2, 3], dtype=float)
d2 = np.array([1, 2, 3], dtype=complex)
d1 + d2
array([ 2.+0.j, 4.+0.j, 6.+0.j])
(d1 + d2).dtype
dtype('complex128')
np.sqrt(np.array([-1, 0, 1]))
/Users/rob/miniconda/envs/py3.6-npm-e2/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in sqrt
"""Entry point for launching an IPython kernel.
array([ nan, 0., 1.])
np.sqrt(np.array([-1, 0, 1], dtype=complex))
array([ 0.+1.j, 0.+0.j, 1.+0.j])
Real and imaginary parts
data = np.array([1, 2, 3], dtype=complex)
data
array([ 1.+0.j, 2.+0.j, 3.+0.j])
data.real
array([ 1., 2., 3.])
data.imag
array([ 0., 0., 0.])
Creating arrays
Arrays created from lists and other array-like objects
np.array([1, 2, 3, 4])
array([1, 2, 3, 4])
data.ndim
1
data.shape
(3,)
np.array([[1, 2], [3, 4]])
array([[1, 2],
[3, 4]])
data.ndim
1
data.shape
(3,)
Arrays filled with constant values
np.zeros((2, 3))
array([[ 0., 0., 0.],
[ 0., 0., 0.]])
np.ones(4)
array([ 1., 1., 1., 1.])
data = np.ones(4)
data.dtype
dtype('float64')
data = np.ones(4, dtype=np.int64)
data.dtype
dtype('int64')
x1 = 5.4 * np.ones(10)
x2 = np.full(10, 5.4)
x1 = np.empty(5)
x1.fill(3.0)
x1
array([ 3., 3., 3., 3., 3.])
x2 = np.full(5, 3.0)
x2
array([ 3., 3., 3., 3., 3.])
Arrays filled with incremental sequences
np.arange(0.0, 10, 1)
array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])
np.linspace(0, 10, 11)
array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])
Arrays filled with logarithmic sequences
np.logspace(0, 2, 5) # 5 data points between 10**0=1 to 10**2=100
array([ 1. , 3.16227766, 10. , 31.6227766 , 100. ])
Mesh-grid arrays
x = np.array([-1, 0, 1])
y = np.array([-2, 0, 2])
X, Y = np.meshgrid(x, y)
X
array([[-1, 0, 1],
[-1, 0, 1],
[-1, 0, 1]])
Y
array([[-2, -2, -2],
[ 0, 0, 0],
[ 2, 2, 2]])
Z = (X + Y) ** 2
Z
array([[9, 4, 1],
[1, 0, 1],
[1, 4, 9]])
Creating uninitialized arrays
np.empty(3, dtype=np.float)
array([ 1., 2., 3.])
Creating arrays with properties of other arrays
def f(x):
y = np.ones_like(x)
# compute with x and y
return y
Creating matrix arrays
np.identity(4)
array([[ 1., 0., 0., 0.],
[ 0., 1., 0., 0.],
[ 0., 0., 1., 0.],
[ 0., 0., 0., 1.]])
np.eye(3, k=1)
array([[ 0., 1., 0.],
[ 0., 0., 1.],
[ 0., 0., 0.]])
np.eye(3, k=-1)
array([[ 0., 0., 0.],
[ 1., 0., 0.],
[ 0., 1., 0.]])
np.diag(np.arange(0, 20, 5))
array([[ 0, 0, 0, 0],
[ 0, 5, 0, 0],
[ 0, 0, 10, 0],
[ 0, 0, 0, 15]])
Index and slicing
One-dimensional arrays
a = np.arange(0, 11)
a
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
a[0] # the first element
0
a[-1] # the last element
10
a[4] # the fifth element, at index 4
4
a[1:-1]
array([1, 2, 3, 4, 5, 6, 7, 8, 9])
a[1:-1:2]
array([1, 3, 5, 7, 9])
a[:5]
array([0, 1, 2, 3, 4])
a[-5:]
array([ 6, 7, 8, 9, 10])
a[::-2]
array([10, 8, 6, 4, 2, 0])
Multidimensional arrays
f = lambda m, n: n + 10 * m
A = np.fromfunction(f, (6, 6), dtype=int)
A
array([[ 0, 1, 2, 3, 4, 5],
[10, 11, 12, 13, 14, 15],
[20, 21, 22, 23, 24, 25],
[30, 31, 32, 33, 34, 35],
[40, 41, 42, 43, 44, 45],
[50, 51, 52, 53, 54, 55]])
A[:, 1] # the second column
array([ 1, 11, 21, 31, 41, 51])
A[1, :] # the second row
array([10, 11, 12, 13, 14, 15])
A[:3, :3] # upper half diagonal block matrix
array([[ 0, 1, 2],
[10, 11, 12],
[20, 21, 22]])
A[3:, :3] # lower left off-diagonal block matrix
array([[30, 31, 32],
[40, 41, 42],
[50, 51, 52]])
A[::2, ::2] # every second element starting from 0, 0
array([[ 0, 2, 4],
[20, 22, 24],
[40, 42, 44]])
A[1::2, 1::3] # every second element starting from 1, 1
array([[11, 14],
[31, 34],
[51, 54]])
Views
B = A[1:5, 1:5]
B
array([[11, 12, 13, 14],
[21, 22, 23, 24],
[31, 32, 33, 34],
[41, 42, 43, 44]])
B[:, :] = 0
A
array([[ 0, 1, 2, 3, 4, 5],
[10, 0, 0, 0, 0, 15],
[20, 0, 0, 0, 0, 25],
[30, 0, 0, 0, 0, 35],
[40, 0, 0, 0, 0, 45],
[50, 51, 52, 53, 54, 55]])
C = B[1:3, 1:3].copy()
C
array([[0, 0],
[0, 0]])
C[:, :] = 1 # this does not affect B since C is a copy of the view B[1:3, 1:3]
C
array([[1, 1],
[1, 1]])
B
array([[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]])
Fancy indexing and Boolean-valued indexing
A = np.linspace(0, 1, 11)
A[np.array([0, 2, 4])]
array([ 0. , 0.2, 0.4])
A[[0, 2, 4]]
array([ 0. , 0.2, 0.4])
A > 0.5
array([False, False, False, False, False, False, True, True, True,
True, True], dtype=bool)
A[A > 0.5]
array([ 0.6, 0.7, 0.8, 0.9, 1. ])
A = np.arange(10)
A
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
indices = [2, 4, 6]
B = A[indices]
B[0] = -1 # this does not affect A
A
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
A[indices] = -1
A
array([ 0, 1, -1, 3, -1, 5, -1, 7, 8, 9])
A = np.arange(10)
B = A[A > 5]
B[0] = -1 # this does not affect A
A
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
A[A > 5] = -1
A
array([ 0, 1, 2, 3, 4, 5, -1, -1, -1, -1])
Reshaping and resizing
data = np.array([[1, 2], [3, 4]])
np.reshape(data, (1, 4))
array([[1, 2, 3, 4]])
data.reshape(4)
array([1, 2, 3, 4])
data = np.array([[1, 2], [3, 4]])
data
array([[1, 2],
[3, 4]])
data.flatten()
array([1, 2, 3, 4])
data.flatten().shape
(4,)
data = np.arange(0, 5)
column = data[:, np.newaxis]
column
array([[0],
[1],
[2],
[3],
[4]])
row = data[np.newaxis, :]
row
array([[0, 1, 2, 3, 4]])
data = np.arange(5)
data
array([0, 1, 2, 3, 4])
np.vstack((data, data, data))
array([[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4]])
data = np.arange(5)
data
array([0, 1, 2, 3, 4])
np.hstack((data, data, data))
array([0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4])
data = data[:, np.newaxis]
np.hstack((data, data, data))
array([[0, 0, 0],
[1, 1, 1],
[2, 2, 2],
[3, 3, 3],
[4, 4, 4]])
Vectorized expressions
Arithmetic operations
x = np.array([[1, 2], [3, 4]])
y = np.array([[5, 6], [7, 8]])
x + y
array([[ 6, 8],
[10, 12]])
y - x
array([[4, 4],
[4, 4]])
x * y
array([[ 5, 12],
[21, 32]])
y / x
array([[ 5. , 3. ],
[ 2.33333333, 2. ]])
x * 2
array([[2, 4],
[6, 8]])
2 ** x
array([[ 2, 4],
[ 8, 16]])
y / 2
array([[ 2.5, 3. ],
[ 3.5, 4. ]])
(y / 2).dtype
dtype('float64')
x = np.array([1, 2, 3, 4]).reshape(2,2)
z = np.array([1, 2, 3, 4])
x / z
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-146-c325a0617380> in <module>()
----> 1 x / z
ValueError: operands could not be broadcast together with shapes (2,2) (4,)
z = np.array([[2, 4]])
z
array([[2, 4]])
z.shape
(1, 2)
x / z
array([[ 0.5, 0.5],
[ 1.5, 1. ]])
zz = np.concatenate([z, z], axis=0)
zz
array([[2, 4],
[2, 4]])
x / zz
array([[ 0.5, 0.5],
[ 1.5, 1. ]])
z = np.array([[2], [4]])
z.shape
(2, 1)
x / z
array([[ 0.5 , 1. ],
[ 0.75, 1. ]])
zz = np.concatenate([z, z], axis=1)
zz
array([[2, 2],
[4, 4]])
x / zz
array([[ 0.5 , 1. ],
[ 0.75, 1. ]])
x = np.array([[1, 3], [2, 4]])
x = x + y
x
array([[ 6, 9],
[ 9, 12]])
x = np.array([[1, 3], [2, 4]])
x += y
x
array([[ 6, 9],
[ 9, 12]])
Elementwise functions
x = np.linspace(-1, 1, 11)
x
array([-1. , -0.8, -0.6, -0.4, -0.2, 0. , 0.2, 0.4, 0.6, 0.8, 1. ])
y = np.sin(np.pi * x)
np.round(y, decimals=4)
array([-0. , -0.5878, -0.9511, -0.9511, -0.5878, 0. , 0.5878,
0.9511, 0.9511, 0.5878, 0. ])
np.add(np.sin(x) ** 2, np.cos(x) ** 2)
array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
np.sin(x) ** 2 + np.cos(x) ** 2
array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
def heaviside(x):
return 1 if x > 0 else 0
heaviside(-1)
0
heaviside(1.5)
1
x = np.linspace(-5, 5, 11)
heaviside(x)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-172-52f34380192d> in <module>()
----> 1 heaviside(x)
<ipython-input-168-68cdb49d4533> in heaviside(x)
1 def heaviside(x):
----> 2 return 1 if x > 0 else 0
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
heaviside = np.vectorize(heaviside)
heaviside(x)
array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1])
def heaviside(x):
return 1.0 * (x > 0)
Aggregate functions
data = np.random.normal(size=(15,15))
np.mean(data)
-0.090472669644641304
data.mean()
-0.090472669644641304
data = np.random.normal(size=(5, 10, 15))
data.sum(axis=0).shape
(10, 15)
data.sum(axis=(0, 2)).shape
(10,)
data.sum()
25.869033203350682
data = np.arange(1,10).reshape(3,3)
data
array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
data.sum()
45
data.sum(axis=0)
array([12, 15, 18])
data.sum(axis=1)
array([ 6, 15, 24])
Boolean arrays and conditional expressions
a = np.array([1, 2, 3, 4])
b = np.array([4, 3, 2, 1])
a < b
array([ True, True, False, False], dtype=bool)
np.all(a < b)
False
np.any(a < b)
True
if np.all(a < b):
print("All elements in a are smaller than their corresponding element in b")
elif np.any(a < b):
print("Some elements in a are smaller than their corresponding elemment in b")
else:
print("All elements in b are smaller than their corresponding element in a")
Some elements in a are smaller than their corresponding elemment in b
x = np.array([-2, -1, 0, 1, 2])
x > 0
array([False, False, False, True, True], dtype=bool)
1 * (x > 0)
array([0, 0, 0, 1, 1])
x * (x > 0)
array([0, 0, 0, 1, 2])
def pulse(x, position, height, width):
return height * (x >= position) * (x <= (position + width))
x = np.linspace(-5, 5, 11)
pulse(x, position=-2, height=1, width=5)
array([0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0])
pulse(x, position=1, height=1, width=5)
array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1])
def pulse(x, position, height, width):
return height * np.logical_and(x >= position, x <= (position + width))
x = np.linspace(-4, 4, 9)
np.where(x < 0, x**2, x**3)
array([ 16., 9., 4., 1., 0., 1., 8., 27., 64.])
np.select([x < -1, x < 2, x >= 2],
[x**2 , x**3 , x**4])
array([ 16., 9., 4., -1., 0., 1., 16., 81., 256.])
np.choose([0, 0, 0, 1, 1, 1, 2, 2, 2],
[x**2, x**3, x**4])
array([ 16., 9., 4., -1., 0., 1., 16., 81., 256.])
x[abs(x) > 2]
array([-4., -3., 3., 4.])
np.nonzero(abs(x) > 2)
(array([0, 1, 7, 8]),)
x[np.nonzero(abs(x) > 2)]
array([-4., -3., 3., 4.])
Set operations
a = np.unique([1,2,3,3])
b = np.unique([2,3,4,4,5,6,5])
np.in1d(a, b)
array([False, True, True], dtype=bool)
1 in a
True
1 in b
False
np.all(np.in1d(a, b))
False
np.union1d(a, b)
array([1, 2, 3, 4, 5, 6])
np.intersect1d(a, b)
array([2, 3])
np.setdiff1d(a, b)
array([1])
np.setdiff1d(b, a)
array([4, 5, 6])
Operations on arrays
data = np.arange(9).reshape(3, 3)
data
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
np.transpose(data)
array([[0, 3, 6],
[1, 4, 7],
[2, 5, 8]])
data = np.random.randn(1, 2, 3, 4, 5)
data.shape
(1, 2, 3, 4, 5)
data.T.shape
(5, 4, 3, 2, 1)
Matrix and vector operations
A = np.arange(1, 7).reshape(2, 3)
A
array([[1, 2, 3],
[4, 5, 6]])
B = np.arange(1, 7).reshape(3, 2)
B
array([[1, 2],
[3, 4],
[5, 6]])
np.dot(A, B)
array([[22, 28],
[49, 64]])
np.dot(B, A)
array([[ 9, 12, 15],
[19, 26, 33],
[29, 40, 51]])
A = np.arange(9).reshape(3, 3)
A
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
x = np.arange(3)
x
array([0, 1, 2])
np.dot(A, x)
array([ 5, 14, 23])
A.dot(x)
array([ 5, 14, 23])
A = np.random.rand(3,3)
B = np.random.rand(3,3)
Ap = B @ A @ np.linalg.inv(B)
Ap
array([[-0.15628289, 1.01950542, 0.54368897],
[ 0.1421288 , 0.60925824, 0.90360847],
[ 0.03321325, 0.63335306, 0.47040536]])
Ap = np.dot(B, np.dot(A, np.linalg.inv(B)))
Ap
array([[-0.15628289, 1.01950542, 0.54368897],
[ 0.1421288 , 0.60925824, 0.90360847],
[ 0.03321325, 0.63335306, 0.47040536]])
Ap = B.dot(A.dot(np.linalg.inv(B)))
Ap
array([[-0.15628289, 1.01950542, 0.54368897],
[ 0.1421288 , 0.60925824, 0.90360847],
[ 0.03321325, 0.63335306, 0.47040536]])
A = np.matrix(A)
B = np.matrix(B)
Ap = B * A * B.I
A = np.asmatrix(A)
B = np.asmatrix(B)
Ap = B * A * B.I
Ap = np.asarray(Ap)
Ap
array([[-0.15628289, 1.01950542, 0.54368897],
[ 0.1421288 , 0.60925824, 0.90360847],
[ 0.03321325, 0.63335306, 0.47040536]])
np.inner(x, x)
5
np.dot(x, x)
5
y = x[:, np.newaxis]
y
array([[0],
[1],
[2]])
np.dot(y.T, y)
array([[5]])
x = np.array([1, 2, 3])
np.outer(x, x)
array([[1, 2, 3],
[2, 4, 6],
[3, 6, 9]])
np.kron(x, x)
array([1, 2, 3, 2, 4, 6, 3, 6, 9])
np.kron(x[:, np.newaxis], x[np.newaxis, :])
array([[1, 2, 3],
[2, 4, 6],
[3, 6, 9]])
np.kron(np.ones((2,2)), np.identity(2))
array([[ 1., 0., 1., 0.],
[ 0., 1., 0., 1.],
[ 1., 0., 1., 0.],
[ 0., 1., 0., 1.]])
np.kron(np.identity(2), np.ones((2,2)))
array([[ 1., 1., 0., 0.],
[ 1., 1., 0., 0.],
[ 0., 0., 1., 1.],
[ 0., 0., 1., 1.]])
x = np.array([1, 2, 3, 4])
y = np.array([5, 6, 7, 8])
np.einsum("n,n", x, y)
70
np.inner(x, y)
70
A = np.arange(9).reshape(3, 3)
B = A.T
np.einsum("mk,kn", A, B)
array([[ 5, 14, 23],
[ 14, 50, 86],
[ 23, 86, 149]])
np.alltrue(np.einsum("mk,kn", A, B) == np.dot(A, B))
True
Versions
%reload_ext version_information
%version_information numpy