# Numbers

The `numpy` array is the foundation of essentially all numerical computing in Python, so it is important to understand the array and how to use it well.

## Learning objectives

1. Attributes of an array
2. How to create vectors, matrices, tensors
3. How to index and slice arrays
4. Generating random arrays and sampling
5. Universal functions, vectorization and matrix multiplication
6. Array axes and marginal calculations
7. Broadcasting
8. Masking
9. Combining and splitting arrays
10. Vectorizing loops


```python
import numpy as np
```

## The `ndarray`: Vectors, matrices and tenosrs

dtype, shape, strides

### Vector


```python
x = np.array([1,2,3])
x
```


```python
type(x)
```


```python
x.dtype
```


```python
x.shape
```


```python
x.strides
```

### Matrix


```python
x = np.array([[1,2,3], [4,5,6]], dtype=np.int32)
x
```


```python
x.dtype
```


```python
x.shape
```


```python
x.strides
```

### Tensor


```python
x = np.arange(24).reshape((2,3,4))
```


```python
x
```

## Creating `ndarray`s

### From a file


```python
%%file numbers.txt
a,b,c # can also skip headers
1,2,3
4,5,6
```


```python
np.loadtxt('numbers.txt', dtype='int', delimiter=',',
           skiprows=1, comments='#')
```

### From Python lists or tuples


```python
np.array([
    [1,2,3],
    [4,5,6]
])
```

### From ranges

arange, linspace, logspace


```python
np.arange(1, 7).reshape((2,3))
```


```python
np.linspace(1, 10, 4)
```


```python
np.logspace(0, 4, 5, dtype='int')
```

### From a function

`fromfunciton`


```python
np.fromfunction(lambda i, j: i*3 + j + 1, (2,3))
```


```python
np.fromfunction(lambda i, j: (i-2)**2 + (j-2)**2, (5,5), dtype='int')
```

#### How to visualize `fromfunction` 


```python
j = np.repeat([np.arange(5)], 5, axis=0)
i = j.T
```


```python
i
```


```python
j
```


```python
(i-2)**2 + (j-2)**2
```

#### Using element-wise functions in `fromfunction`


```python
np.fromfunction(lambda i, j: np.where(i==j,0, -1), (5,5))
```


```python
np.fromfunction(lambda i, j: np.where(i<j, 1, np.where(i==j,0, -1)), (5,5))
```


```python
np.fromfunction(lambda i, j: np.minimum(i,j), (5,5), dtype='int')
```


```python
np.fromfunction(lambda i, j: np.maximum(i,j), (5,5), dtype='int')
```

### From special constructors

zeros, ones, eye, diag


```python
np.zeros((2,3))
```


```python
np.ones((2,3))
```


```python
np.eye(3)
```


```python
np.eye(3, 4)
```


```python
np.eye(4, k=-1)
```


```python
np.diag([1,2,3,4])
```


```python
np.diag([1,2,3,4], k=1)
```

### From random variables

#### Convenience functions

rand, randn


```python
np.random.rand(2,3)
```


```python
np.random.randn(2,3)
```

#### Distributions

uniform, normal, randint, poisson, multinomial, multivariate_ normal


```python
np.random.uniform(0, 1, (2,3))
```


```python
np.random.normal(0, 1, (2,3))
```


```python
np.random.randint(0, 10, (4,5))
```


```python
np.random.poisson(10, (4,5))
```


```python
np.random.multinomial(n=5, pvals=np.ones(5)/5, size=8)
```


```python
np.random.multivariate_normal(mean=[10,20,30], cov=np.eye(3), size=4)
```

### Sampling using `choice`

Works much like the R `sample` function.


```python
x = np.random.permutation(list('ABCDEF'))
```


```python
x
```


```python
np.random.choice(x, 3)
```


```python
np.random.choice(x, 10)
```


```python
try:
    np.random.choice(x, 10, replace=False)
except ValueError as e:
    print(e)
```

## Indexing 


```python
x = np.arange(20).reshape((4,5))
x
```

### Extracing a scalar


```python
x[1,1]
```

### Extracting a vector


```python
x[1]
```

### Using slices


```python
x[1,:]
```


```python
x[:,1]
```


```python
x[1:3,1:3]
```

### Using slices with strides


```python
x[::2,::2]
```

### Extrcting blocks with arbitrary row and column lists (fancy indexing)

`np.ix_`


```python
x[:, [0,3]]
```

Warning: Fancy indexing can only be used for 1 dimension at a time.

In the example below, `numpy` treats the arguments as *paired* coordinates, and returns the values at (0,0) and (2,3).


```python
x[[0,2],[0,3]]
```

Use the helper `np.ix_` to extract arbitrary blocks.


```python
x[np.ix_([0,2], [0,3])]
```

### A slice is a view, not a copy

**Warning**

```python
b = a[:]
```

makes a copy if `a` is a list but not if `a` is a numpy array


```python
a1 = list(range(3))
a2 = np.arange(3)
```


```python
b = a1[:]
b[1] = 9
a1
```


```python
b = a2[:]
b[1] = 9
a2
```


```python
x
```


```python
y = x[1:-1, 1:-1]
y
```


```python
y *= 10
```


```python
y
```


```python
x
```

Use the copy method to convert a view to a copy


```python
z = x[1:-1, 1:-1].copy()
```


```python
z
```


```python
z[:] = 0
```


```python
z
```


```python
x
```

### Boolean indexing


```python
x[x % 2 == 0]
```


```python
x [x > 3]
```

### Functions that return indexes


```python
idx = np.nonzero(x)
idx
```


```python
x[idx]
```


```python
idx = np.where(x > 3)
idx
```


```python
x[idx]
```

## Universal functions


```python
x
```

Operations


```python
x + x
```

Element-wise functions


```python
np.log1p(x)
```


```python
x.clip(10, 100)
```

Scans


```python
np.cumsum(x, axis=1)
```

Reductions


```python
np.sum(x)
```


```python
x.prod()
```

## Margins and the `axis` argument


```python
x
```

The 0th axis has 4 items, the 1st axis has 5 items.


```python
x.shape
```


```python
x.mean()
```

### Marginalizing out the 0th axis = column summaries


```python
x.mean(axis=0)
```

### Marginalizing out the 1st axis = row summaries


```python
x.mean(axis=1)
```

Note marginalizing out the last axis is a common default.


```python
x.mean(axis=-1)
```

### Marginalization works for higher dimensions in the same way


```python
x = np.random.random((2,3,4))
x
```


```python
x.shape
```


```python
x.mean(axis=0).shape
```


```python
x.mean(axis=1).shape
```


```python
x.mean(axis=2).shape
```


```python
x.mean(axis=(0,1)).shape
```


```python

x.mean(axis=(0,2)).shape
```


```python
x.mean(axis=(1,2)).shape
```

## Broadcasting

Broadcasting is what happens when `numpy` tries to perform binary operations on two arrays with different shapes. In general, shapes are *promoted* to make the arrays compatible using the following rule

- For each axis from highest to lowest
    - If both dimensions are the same, do nothing
    - If one of the dimensions is 1 or None and the other is $k$, promote to $k$
    - Otherwise print error message


```python
x = np.zeros((3,2))
x.shape
```




    (3, 2)




```python
x
```




    array([[0., 0.],
           [0., 0.],
           [0., 0.]])



Shapes are compatible


```python
y = np.ones(2)
y.shape
```




    (2,)




```python
x + y
```




    array([[1., 1.],
           [1., 1.],
           [1., 1.]])



Shapes are compatible


```python
y = np.ones((1,2))
y.shape
```


```python
x + y
```

Shapes are incompatible but can be made compaible by adding empty dimension


```python
y = np.ones(3)
y.shape
```




    (3,)




```python
try:
    x + y
except ValueError as e:
    print(e)
```

    operands could not be broadcast together with shapes (3,2) (3,) 



```python
y[:, None].shape
```




    (3, 1)




```python
x + y[:, None]
```




    array([[1., 1.],
           [1., 1.],
           [1., 1.]])



Shapes are incompatible


```python
y = np.ones((2,2))
y.shape
```


```python
try:
    x + y
except ValueError as e:
    print(e)
```

### More examples of broadcasting


```python
x1 = np.arange(12)
```


```python
x1
```


```python
x1 * 10
```


```python
x2 = np.random.randint(0,10,(3,4))
```


```python
x2
```


```python
x2 * 10
```


```python
x2.shape
```

### Column-wise broadcasting


```python
mu = np.mean(x2, axis=0)
mu.shape
```


```python
x2 - mu
```


```python
(x2 - mu).mean(axis=0)
```

### Row wise broadcasting


```python
mu = np.mean(x2, axis=1)
mu.shape
```


```python
try:
    x2 - mu
except ValueError as e:
    print(e)
```

### We can add a "dummy" axis using None or `np.newaxis`


```python
mu[:, None].shape
```


```python
x2 - mu[:, None]
```


```python
x2 - mu[:, np.newaxis]
```


```python
np.mean(x2 - mu[:, None], axis=1)
```

#### Reshaping works too


```python
x2 - mu.reshape((-1,1))
```

#### Exercise in broadcasting

Creating a 12 by 12 multiplication table


```python
x = np.arange(1, 13)
x[:,None] * x[None,:]
```

Scaling to have zero mean and unit standard devation for each feature.


```python
x = np.random.normal(10, 5,(3,4))
x
```

Scaling column-wise


```python
(x - x.mean(axis=0))/x.std(axis=0)
```

Scaling row-wise


```python
(x - x.mean(axis=1)[:, None])/x.std(axis=1)[:, None]
```

## Masking

- [Ref](https://docs.scipy.org/doc/numpy/reference/maskedarray.generic.html)


```python
import numpy.ma as ma
```


```python
x = np.arange(20).reshape(4,5)
```


```python
x
```


```python
mask = x % 2 == 0
mask
```

- Note that values that are True in the mask are not used in the array
- Values that are False are *not* masked and so remain
- So the above mask keeps only the *odd* numbers in the array `x`


```python
m = ma.masked_array(x, mask)
```


```python
m
```


```python
m.data
```


```python
m.mask
```


```python
m.sum(axis=0).data
```


```python
m.sum(axis=1).data
```

### Often used with missing value sentinels


```python
import warnings

with warnings.catch_warnings():
    warnings.simplefilter('ignore', RuntimeWarning)
    x1 = x / mask
```


```python
x1
```


```python
x1.sum()
```


```python
x2 = ma.masked_invalid(x1)
x2
```


```python
x2.data
```


```python
x2.mask
```


```python
x2.sum()
```


```python
x2.filled(0)
```

## Combining `ndarray`s


```python
x1 = np.zeros((3,4))
x2 = np.ones((3,5))
x3 = np.eye(4)
```


```python
x1
```


```python
x2
```


```python
x3
```

### Binding rows when number of columns is the same


```python
np.r_[x1, x3]
```

### Binding columns when number of rows is the same


```python
np.c_[x1, x2]
```

### You can combine more than 2 at a time


```python
np.c_[x1, x2, x1]
```

### Stacking


```python
np.vstack([x1, x3])
```


```python
np.hstack([x1, x2])
```


```python
np.dstack([x2, 2*x2, 3*x2])
```

### Generic stack with axis argument


```python
np.stack([x2, 2*x2, 3*x2], axis=0)
```


```python
np.stack([x2, 2*x2, 3*x2], axis=1)
```


```python
np.stack([x2, 2*x2, 3*x2], axis=2)
```

### Repetition and tiling

#### For a vector


```python
x = np.array([1,2,3])
```


```python
np.repeat(x, 3)
```


```python
np.tile(x, 3)
```

#### For a matrix


```python
x = np.arange(6).reshape((2,3))
x
```


```python
np.repeat(x, 3)
```


```python
np.repeat(x, 3, axis=0)
```


```python
np.repeat(x, 3, axis=1)
```


```python
np.tile(x, (3,2))
```

## Splitting `ndarray`s


```python
x = np.arange(32).reshape((4,8))
```


```python
x
```


```python
np.split(x, 4)
```


```python
np.split(x, 4, axis=1)
```

## Saving and loading arrays


```python
x = np.arange(16).reshape(4,4)
y = np.arange(20).reshape(-1,4)
```


```python
np.save('x.npy', x)
```


```python
np.load('x.npy')
```


```python
np.savez('xy.npz', x=x, y=y)
```


```python
arr = np.load('xy.npz')
```


```python
arr['x']
```


```python
arr['y']
```


```python
import os
```


```python
os.remove('x.npy')
```


```python
os.remove('xy.npz')
```

### Einstein summation notation

See https://obilaniu6266h16.wordpress.com/2016/02/04/einstein-summation-in-numpy/ for a comprehensible explanation before reading the `numpy.einsum` help documentation. Then try out several examples.

On the first try, ignore implicit outputs and broadcasting `...`.


```python
x = np.arange(1,10).reshape(3,3)
x
```


```python
np.einsum('ii', x)
```


```python
np.einsum('ji', x)
```


```python
np.einsum('ii->i', x)
```


```python
np.einsum('ij->i', x)
```


```python
np.einsum('ij->j', x)
```


```python
np.einsum('mn, np -> mp', x, x)
```


```python
np.einsum('...i,...i', x, x)
```


```python
np.einsum('j...,j...', x,x)
```

## Vectorization

### Example 1

The operators and functions (ufuncs) in Python are vectorized, and will work element-wise over all entries in an `ndarray`.


```python
xs = np.zeros(10, dtype='int')
for i in range(10):
    xs[i] = i**2
xs
```


```python
xs = np.arange(10)**2
xs
```

Using ufuncs


```python
np.sqrt(xs)
```


```python
np.log1p(xs)
```

### Example 2

Scalar product.


```python
n = 10

xs = np.random.rand(n)
ys = np.random.rand(n)

s = 0
for i in range(n):
    s += xs[i] * ys[i]
s
```


```python
np.dot(xs, ys)
```


```python
xs @ ys
```


```python
np.einsum('i,i->', xs, ys)
```

### Example 3

\begin{align}
y_0 &= \alpha + \beta_1 x_1 + \beta_2 x_2 \\
y_1 &= \alpha + \beta_1 x_1 + \beta_2 x_2 \\
y_2 &= \alpha + \beta_1 x_1 + \beta_2 x_2 \\
\end{align}





```python
m = 3
n = 2

alpha = np.random.rand(1)
betas = np.random.rand(n,1)
xs = np.random.rand(m,n)
```


```python
alpha
```


```python
betas
```


```python
xs
```

### Using loops


```python
ys = np.zeros((m,1))
for i in range(m):
    ys[i] = alpha
    for j in range(n):
        ys[i] += betas[j] * xs[i,j]
ys
```

### Removing inner loop


```python
ys = np.zeros((m,1))
for i in range(m):
    ys[i] = alpha + xs[i,:].T @ betas
ys
```

### Removing all loops


```python
ys = alpha + xs @ betas
ys
```

### Alternative approach

The calculaiton with explicit intercepts and coefficients is common in deep learning, where $\alpha$ is called the bias ($b$) and $\beta$ are called the weights ($w$), and each equation is $y[i] = b + w[i]*x[i]$.

It is common in statisiics to use an augmented matrix in which the first column is all ones, so that all that is needed is a single matrix multiplicaiotn.


```python
X = np.c_[np.ones(m), xs]
X
```


```python
alpha
```


```python
betas
```


```python
betas_ = np.concatenate([[alpha], betas])
betas_
```


```python
ys = X @ betas_
ys
```


```python
X.shape, betas_.shape
```


```python
np.einsum('ij,jk->ik', X, betas_)
```

### Simulating diffusion

See how we use views to avoid making (slow) copies of the matrix at each step.


```python
%matplotlib inline
import matplotlib.pyplot as plt
```


```python
w = 100
h = 100
x = np.zeros((w+2,h+2), dtype='float')
x[(w//2-1):(w//2+2), (h//2-1):(h//2+2)] = 1

wts = np.ones(5)/5

for i in range(41):
    if i % 10 == 0:    
        plt.figure()
        plt.imshow(x[1:-1, 1:-1], interpolation='nearest')
        
    center = x[1:-1, 1:-1]
    left = x[:-2, 1:-1]
    right = x[2:, 1:-1]
    bottom = x[1:-1, :-2]
    top = x[1:-1, 2:]
    nbrs = np.dstack([center, left, right, bottom, top])
    x = np.sum(wts * nbrs, axis=-1)
```


```python

```

0a3f52669567ecb999c883f15313810b8704179b
ipynb
Jupyter Notebook
notebooks/S03_Numpy.ipynb
itachi4869/sta-663-2021
6f7a50f4a95207a26b01afa4439d992d767a1387
MIT
notebooks/S03_Numpy.ipynb
itachi4869/sta-663-2021
6f7a50f4a95207a26b01afa4439d992d767a1387
MIT
notebooks/S03_Numpy.ipynb
itachi4869/sta-663-2021
6f7a50f4a95207a26b01afa4439d992d767a1387
MIT
Qwen/Qwen-72B
1. YES
2. YES