Links to materials the tutorial is based on can be found in the end of the notebook, in the "References" section. And sure, I'm not going to retell NumPy quickstart tutorial.

Let's first compare amounts of bytes taken by standard Python list and identical NumPy array.

In [ ]:

```
import sys
import numpy as np
```

In [ ]:

```
mb = 1024 * 1024
python_list = list(range(0, 1000000))
numpy_array = np.array(range(0, 1000000))
print("Python list size: {0:.3f} MB".format(sys.getsizeof(python_list) / mb))
print("Numpy array size: {0:.3f} MB".format(numpy_array.nbytes / mb, "MB"))
print("ArraySize/ListSize ratio: {0:.2f}".format(numpy_array.nbytes / sys.getsizeof(python_list)))
```

In [ ]:

```
mixed_list = [1, 2, 3, 4.0, 5.0, 'abc', 'def', True, False]
```

Creating mixed int and string array:

In [ ]:

```
np_array_mixed = np.array([i if i % 2 == 0 else str(i) for i in range(0, 1000000)])
np_array_mixed_size = np_array_mixed.nbytes
print("Mixed numpy array size:", np_array_mixed_size / mb, "MB")
```

And when storing string only:

In [ ]:

```
np_array_str = np.array([str(i) for i in range(0, 1000000)])
np_array_str_size = np_array_str.nbytes
print("Strings-only numpy array size:", np_array_str_size / mb, "MB")
```

In [ ]:

```
print('String array contains type:', np_array_str.dtype)
print('"Mixed" array contains type:', np_array_mixed.dtype)
print('"Mixed" array elements:"', np_array_mixed[:6])
```

In [ ]:

```
nums_list = list(range(1, 10000000))
numpy_nums = np.array(nums_list)
```

In [ ]:

```
%%time
reciprocal_list = []
for num in nums_list:
reciprocal_list.append(1 / num)
```

Python loops shows as quite slow. Maybe list comprehension are quicker?

In [ ]:

```
%%time
reciprocal_list = [1 / num for num in nums_list]
```

A little bit. What about NumPy?

In [ ]:

```
%%time
reciprocal_array = 1 / numpy_nums
```

Much more better!

In [ ]:

```
%%time
nums_list = list(range(1, 10000000))
```

In [ ]:

```
%%time
numpy_nums = np.array(nums_list)
```

In [ ]:

```
%%time
numpy_nums = np.array(range(1, 10000000))
```

In [ ]:

```
%%time
numpy_nums = np.arange(1, 10000000)
```

In [ ]:

```
%%time
reciprocal_array = 1 / numpy_nums
```

In [ ]:

```
%%time
point = [1.0, 1.0]
other_points = list([float(x), x + 1.0] for x in range(10000000))
dists = [((point[0] - op[0])**2 + (point[1] - op[1])**2)**0.5 for op in other_points]
```

In [ ]:

```
dists[:10]
```

In [ ]:

```
%%time
np_point = np.ones(2)
np_other_points = np.vstack([np.arange(0.0, 10000000.0), np.arange(1.0, 10000001.0)]).T
np_dists = np.sum((np_point - np_other_points)**2, axis=1)**0.5
```

In [ ]:

```
np_dists[:10]
```

The results are the same, but the performance advantage is significant.

One of the things that should be noticed about the example above is that it handles matrixes with different dimensions. The np_point has shape (2,) and np_other_points - (10000000, 2). Nevertheless, element-wise operations are performed with them, producing a matrix with shape (10000000,). That is possible by the virtue of ** matrix broadcasting** mechanism.

Matrix broadcasting in NumPy is a set of rules that allowes to expand two matrixes with different dimentions to match shapes of each other, in order to perform element-by-element operations. What is important, this set of rules is a "virtual" mechanism, that just allows to understand how matrixes will interact. No real expansion and memory allocation is performed.

Its simplest case is summing a matrix with a scalar number, that will be added to each element of the matrix:

In [ ]:

```
M = np.ones((3, 3))
print("M:")
print(M)
scalar = 2
M_added = M + scalar
print()
print("M +", scalar, ":")
print(M_added)
```

A matrix can be summed with a vector in the similar way:

In [ ]:

```
M = np.ones((3, 3))
print(M)
v = np.array([0, 1, 2])
print()
print(v)
M_added = M + v
print()
print(M_added)
```

In a more complex case the both matrixes are broadcasted:

In [ ]:

```
a = np.arange(3)
b = np.arange(3)[:, np.newaxis]
print("a:", a)
print("b:\n", b)
summed = a + b
print("\nSummed:")
print(summed)
```

**Rule 1:** If the two arrays differ in their number of dimensions, the shape of the one with fewer dimensions is padded with ones on its leading (left) side.

**Rule 2:** If the shape of the two arrays does not match in any dimension, the array with shape equal to 1 in that dimension is stretched to match the other shape.

**Rule 3:** If in any dimension the sizes disagree and neither is equal to 1, an error is raised.

So, adding a two-dimensional array to a one-dimensional array is performed this way:

In [ ]:

```
M = np.ones((2, 3))
a = np.arange(3)
print("M:", M.shape)
print(M)
print()
print("a:", a.shape)
print(a)
```

The shape of a is pad on the left, as it has fewer dimensions:

M.shape -> (2, 3)

a.shape -> (1, 3)

The the 1-dimension of a is stretched to match M:

M.shape -> (2, 3)

a.shape -> (2, 3)

Now the shapes match and matrixes can be summed:

In [ ]:

```
M + a
```

That's what happens when both arrays need to be broadcasted:

In [ ]:

```
a = np.arange(3).reshape((3, 1))
b = np.arange(3)
print("a:", a.shape)
print(a)
print()
print("b:", b.shape)
print(b)
```

The shape of b is pad left with ones:

a.shape -> (3, 1)

b.shape -> (1, 3)

Then both matrixes has dimension to be expanded:

a.shape -> (3, 3)

b.shape -> (3, 3)

Then the matrixes can be easily summed:

In [ ]:

```
a + b
```

However, broadcasting is not always possible:

In [ ]:

```
M = np.ones((3, 2))
print("M:", M.shape)
print(M)
print()
a = np.arange(3)
print("a:", a.shape)
print(a)
print()
```

M.shape -> (3, 2)

a.shape -> (1, 3)

M.shape -> (3, 2)

a.shape -> (3, 3)

Matrixes shape does not match, the error is raised when trying to perform operations with them:

In [ ]:

```
M + a
```

*right* manually. This way the Rule 1 will be skipped and according to the Rule 2 NumPy will just expand matrix to the needed size:

In [ ]:

```
a = a[:, np.newaxis]
print("a:", a.shape)
print(a)
print()
```

In [ ]:

```
M + a
```

In [ ]:

```
np.newaxis is None
```

*np.newaxis* constant is useful when converting a 1D array into a row vector or a column vector, by adding new dimensions from left or right side:

In [ ]:

```
arr = np.arange(3)
arr.shape
```

In [ ]:

```
row_vec = arr[np.newaxis, :]
print(row_vec.shape)
print(row_vec)
```

In [ ]:

```
col_vec = arr[:, np.newaxis]
print(col_vec.shape)
print(col_vec)
```

*np.reshape(-1, 1)* and *np.reshape(1, -1)* down to minor implementation details. But the np.newaxis allows to stack dimentions using slice syntax, without specifying original shape:

In [ ]:

```
M = np.ones((5, 5))
```

In [ ]:

```
M[np.newaxis, :, :, np.newaxis, np.newaxis].shape
```

In [ ]:

```
M[np.newaxis, ..., np.newaxis, np.newaxis].shape
```

In [ ]:

```
M.reshape(1, -1, -1, 1, 1 ).shape
```

This will work:

In [ ]:

```
M.reshape(1, *M.shape, 1, 1 ).shape
```

but doesn't it look a little bit clumsy?

**Pytorch** allows to initialize its "Tensors" from numpy arrays, but often requires input in the form "minibatch × in_channels × iW" or "minibatch × in_channels × iH × iW" (torch.nn.functional). There, minibatch and in_channels can be equal to 1, but they must be present.

In [ ]:

```
A = np.full((2, 2), 2)
print("A:")
print(A)
print()
B = np.full((2, 2), 3)
print("B:")
print(B)
print()
print("A*B:")
print(A*B)
```

*matrix product*, that is defined this way (formula from wikipedia):

$$ A = \begin{pmatrix} a_{11}, a_{12} & \cdots & a_{1m}\\a_{21}, a_{22} & \cdots & a_{2m} \\ \vdots & \ddots & \vdots \\ a_{n1}, a_{n2} & \cdots & a_{nm} \end{pmatrix}, B = \begin{pmatrix} b_{11}, b_{12} & \cdots & b_{1p}\\b_{21}, b_{22} & \cdots & b_{2p} \\ \vdots & \ddots & \vdots \\ b_{m1}, b_{m2} & \cdots & b_{mp} \end{pmatrix} $$

Matrix product **C** = **AB**:
$$C = \begin{pmatrix} c_{11}, c_{12} & \cdots & c_{1p} \\c_{21}, c_{22} & \cdots & c_{2p} \\ \vdots & \ddots & \vdots \\ c_{n1}, c_{n2} & \cdots & c_{np} \end{pmatrix} $$
$c_{ij} = a_{i1}b_{1j} + ... + a_{im}b_{mj} = \sum_{k=1}^{m} a_{ik}b_{kj}$

*numpy.dot()* function or "@" shortcut is used for this purpose in NumPy. That is unpleasant to confuse these two operations, especially when matrix broadcasting exists:

In [ ]:

```
A = np.arange(1, 10).reshape(3,3)
print("A:")
print(A)
print()
B = np.full((3,), 3)
print("B:")
print(B)
print()
print("A*B:")
print(A*B)
print()
print("A@B:")
print(A@B)
```

So, be attentive :)

We used np.meshgrid function somewhere in the course, but without any explanations. That's pity in my opinion, as it is not so easy to grasp what it does from the NumPy documentation:

numpy.meshgrid(*xi, **kwargs)

Return coordinate matrices from coordinate vectors.

Make N-D coordinate arrays for vectorized evaluations of N-D scalar/vector fields over N-D grids, given one-dimensional coordinate arrays x1, x2,…, xn.

Actually, this function just creates a set of grids with coordinates of x, y, etc. on the corresponding grid locations.

In [ ]:

```
import matplotlib.pyplot as plt
from matplotlib import cm
xvalues = np.arange(-1, 1.05, 0.5)
yvalues = np.arange(-1, 1.05, 0.5)
xx, yy = np.meshgrid(xvalues, yvalues)
print(xx)
print()
print(yy)
grid = plt.plot(xx, yy, marker='.', color='k', linestyle='none')
```

In [ ]:

```
x = np.arange(-8, 8, 0.01)
y = np.arange(-8, 8, 0.01)
xx, yy = np.meshgrid(x, y, sparse=True)
z = np.sin(xx**2 + yy**2) / (xx**2 + yy**2)
h = plt.contourf(x,y,z, cmap=cm.PuBu_r)
```

In [ ]:

```
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(10, 8))
ax = fig.gca(projection='3d')
ax.view_init(60, 35)
# Make data.
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)
# Plot the surface.
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm, linewidth=0, antialiased=False)
```

In [ ]:

```
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn; seaborn.set()
rand = np.random.RandomState(42)
X = rand.rand(10, 2)
plt.scatter(X[:, 0], X[:, 1], s=100);
```

In [ ]:

```
dist_sq = np.sum((X[:, np.newaxis, :] - X[np.newaxis, :, :]) ** 2, axis=-1)
```

The dist_qs matrix has zeros on its diagonal, that proves calculation correctness. Distance between ii element is distance between i-th point and itself; that is zero.

In [ ]:

```
dist_sq.diagonal()
```

*np.argsort* function allows to sort elements in the each row and print indexes of other points in the order of their remoteness from the i-th point:

In [ ]:

```
nearests = np.argsort(dist_sq, axis=1)
print(nearests)
```

*np.argpartition* function allows to take them only, without sorting the whole rows:

In [ ]:

```
K = 2
nearest_partition = np.argpartition(dist_sq, K + 1, axis=1)
print(nearest_partition[:, :K+1])
```

That can be visualized:

In [ ]:

```
cmap = plt.get_cmap('viridis')
colors = cmap(np.linspace(0, 1, 10))
plt.scatter(X[:, 0], X[:, 1], s=150, color=colors)
# draw lines from each point to its two nearest neighbors
K = 2
for i in range(X.shape[0]):
for j in nearest_partition[i, :K+1]:
# plot a line from X[i] to X[j]
# the lines colors correspond to outgoing point color, but some lines obviously overlap
# using some zip magic:
plt.plot(*zip(X[j], X[i]), c=colors[i])
```

The Conway's Game of Life is a classical model of cellular automaton model and a zero-player game. Given an initial state, it starts to live its own life by applying the following rules on each step:

- Each cell on a 2D grid is "alive"(1) or "dead"(0)
- Any living cell that has 2 or 3 neighbors survives, else it dies [0,1 or 4+ neighbors]
- Any cell with exactly 3 neighbors becomes alive (if it was dead)

What it takes to implement this game using NumPy arrays? It's tricky, but not so long:

In [ ]:

```
def iterate(Z):
# Count neighbours
N = (Z[0:-2,0:-2] + Z[0:-2,1:-1] + Z[0:-2,2:] +
Z[1:-1,0:-2] + Z[1:-1,2:] +
Z[2: ,0:-2] + Z[2: ,1:-1] + Z[2: ,2:])
# Apply rules
birth = (N==3) & (Z[1:-1,1:-1]==0)
survive = ((N==2) | (N==3)) & (Z[1:-1,1:-1]==1)
Z[...] = 0
Z[1:-1,1:-1][birth | survive] = 1
return Z
```

First, it slices original matrix with 0, 1 and 2 vertical strides, thus getting matrices with "middle" rows, rows shifted by -1 from the "middle" ones, and rows shifted by +1. The similar action is performed with the columns. Combining this steps in different directions creates 8 shifted matrices; their elementwise sum gives amount of "alive" neighbors for every element of N matrix.

Then, by applying game rules to the N and Z matrices the boolean mask matrices "birth" and "survive" can be calculated. After that, NumPy boolean indexing allows to set to "1" only those cells that were born or had survived.

In [ ]:

```
import matplotlib.pyplot as plt
import numpy as np
%matplotlib notebook
plt.ion()
#initialize game field
Z = np.random.choice([0,1],size=(100,100))
fig = plt.figure()
ax = fig.add_subplot(111)
fig.show()
for _ in range(1000):
#update
Z = iterate(Z)
#re-draw image
ax.clear()
ax.imshow(Z,cmap='gray')
fig.canvas.draw()
```

When convolution is applied, for example to an image in computer vision, it takes the matrix with image pixels and slides across it with small window called "filter" or "kernel". The kernel is just little (3x3, 5x5, 7x7 or smth. like this) weight matrix that contains numerical coefficients. Each of these coefficients has to be multiplied by the value of pixel that is currently under the corresponding kernel cell. Then, all such multiplied element are summed, giving the value of an element of the resulting map. After that, the kernel is moved by 1 or more strides horizontaly or verticaly, to the next group of pixels, and the next weighted sum is calculated.

This image illustrates the process of applying the filter and calculating an element of resulting map:

In [ ]:

```
filters = np.array([[1,1,1],
[1,0,1],
[1,1,1]])
```

Let's see if it indeed does the same thing:

In [ ]:

```
from scipy import signal
Z = np.random.choice([0,1],size=(10,10))
N1 = (Z[0:-2,0:-2] + Z[0:-2,1:-1] + Z[0:-2,2:] +
Z[1:-1,0:-2] + Z[1:-1,2:] +
Z[2: ,0:-2] + Z[2: ,1:-1] + Z[2: ,2:])
N2 = signal.convolve2d(Z, filters, mode='valid')
print("Are N1 and N2 identical?")
print(np.array_equal(N1, N2))
```

That's all by this moment. Hopping, it was informative :)

- "Python Data Science Handbook" by Jake VanderPlas
- What is numpy.newaxis and when to use it
- NumPy documentation about numpy.dot, numpy.newaxis, numpy.meshgrid
- Matplotlib 3d surface example
- Exercise 88 from the 100 NumPy exercises