# 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 (arrays with dimension > 2)
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 (when trying to perform same operation on two arrays that are not same size)
8. Masking (tell numpy to ignore certain values)
9. Combining and splitting arrays
10. Vectorizing loops (using vector / matrix operations instead of loops)

Note: There are numpy alternatives, but they tend to use the same or similar syntax


```python
import numpy as np
```

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

dtype, shape, strides

### Vector


```python
# 1D array (which is like a vector)
x = np.array([1,2,3])
x
```




    array([1, 2, 3])




```python
# -1 means numpy should guess the correct number. In this case, we specify one column, and Python chooses 3 rows
# Creates a column vector
np.array([1, 2, 3]).reshape([-1,1])
```




    array([[1],
           [2],
           [3]])




```python
type(x)
```




    numpy.ndarray




```python
x.dtype
```




    dtype('int64')




```python
x.shape
```




    (3,)




```python
# number of bytes
x.strides
```




    (8,)



### Matrix


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




    array([[1, 2, 3],
           [4, 5, 6]], dtype=int32)




```python
x.dtype
```




    dtype('int32')




```python
# 2 rows, 3 columns
x.shape
```




    (2, 3)




```python
x.strides
```




    (12, 4)



### Tensor


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


```python
x
```




    array([[[ 0,  1,  2,  3],
            [ 4,  5,  6,  7],
            [ 8,  9, 10, 11]],
    
           [[12, 13, 14, 15],
            [16, 17, 18, 19],
            [20, 21, 22, 23]]])



## Creating `ndarray`s

### From a file


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

    Writing numbers.txt



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




    array([[1, 2, 3],
           [4, 5, 6]])



### From Python lists or tuples


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




    array([[1, 2, 3],
           [4, 5, 6]])



### From ranges

arange, linspace, logspace


```python
# np.arange is same as range for numpy. Fills rows first by default
np.arange(1, 7).reshape((2,3))
```




    array([[1, 2, 3],
           [4, 5, 6]])




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




    array([ 1.,  4.,  7., 10.])




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




    array([    1,    10,   100,  1000, 10000])



### From a function

`fromfunciton`


```python
# Function says what value each cell should have based on position
np.fromfunction(lambda i, j: i*3 + j + 1, (2,3))
```




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




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




    array([[8, 5, 4, 5, 8],
           [5, 2, 1, 2, 5],
           [4, 1, 0, 1, 4],
           [5, 2, 1, 2, 5],
           [8, 5, 4, 5, 8]])



#### How to visualize `fromfunction` 


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


```python
i
```




    array([[0, 0, 0, 0, 0],
           [1, 1, 1, 1, 1],
           [2, 2, 2, 2, 2],
           [3, 3, 3, 3, 3],
           [4, 4, 4, 4, 4]])




```python
j
```




    array([[0, 1, 2, 3, 4],
           [0, 1, 2, 3, 4],
           [0, 1, 2, 3, 4],
           [0, 1, 2, 3, 4],
           [0, 1, 2, 3, 4]])




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




    array([[8, 5, 4, 5, 8],
           [5, 2, 1, 2, 5],
           [4, 1, 0, 1, 4],
           [5, 2, 1, 2, 5],
           [8, 5, 4, 5, 8]])



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


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




    array([[ 0, -1, -1, -1, -1],
           [-1,  0, -1, -1, -1],
           [-1, -1,  0, -1, -1],
           [-1, -1, -1,  0, -1],
           [-1, -1, -1, -1,  0]])




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




    array([[ 0,  1,  1,  1,  1],
           [-1,  0,  1,  1,  1],
           [-1, -1,  0,  1,  1],
           [-1, -1, -1,  0,  1],
           [-1, -1, -1, -1,  0]])




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




    array([[0, 0, 0, 0, 0],
           [0, 1, 1, 1, 1],
           [0, 1, 2, 2, 2],
           [0, 1, 2, 3, 3],
           [0, 1, 2, 3, 4]])




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




    array([[0, 1, 2, 3, 4],
           [1, 1, 2, 3, 4],
           [2, 2, 2, 3, 4],
           [3, 3, 3, 3, 4],
           [4, 4, 4, 4, 4]])



### From special constructors

zeros, ones, eye, diag


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




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




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




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




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




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




```python
# Doesn't have to be square
np.eye(3, 4)
```




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




```python
# Can offset 1's from diagonal
np.eye(4, k=-1)
```




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




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




    array([[1, 0, 0, 0],
           [0, 2, 0, 0],
           [0, 0, 3, 0],
           [0, 0, 0, 4]])




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




    array([[0, 1, 0, 0, 0],
           [0, 0, 2, 0, 0],
           [0, 0, 0, 3, 0],
           [0, 0, 0, 0, 4],
           [0, 0, 0, 0, 0]])



### From random variables

#### Convenience functions

rand, randn


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




    array([[0.62714498, 0.76026772, 0.77202773],
           [0.02715102, 0.84982795, 0.05837353]])




```python
# Standard normal
np.random.randn(2,3)
```




    array([[-0.97359113, -0.95378846, -1.16320165],
           [ 1.76855437, -0.2479148 ,  0.06691788]])



#### Distributions

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


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




    array([[0.08631122, 0.27358747, 0.53793108],
           [0.3190033 , 0.82258753, 0.2233583 ]])




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




    array([[-1.45531815, -1.20090964, -2.20962917],
           [-1.64417944,  0.42055602, -0.40536003]])




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




    array([[0, 7, 6, 9, 0],
           [4, 6, 3, 1, 9],
           [9, 9, 1, 8, 2],
           [3, 6, 9, 0, 7]])




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




    array([[11, 16,  9, 12,  7],
           [ 9, 12,  6,  9, 15],
           [12, 14,  8,  8, 14],
           [25,  6,  9, 19, 17]])




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




    array([[1, 3, 1, 0, 0],
           [1, 1, 2, 1, 0],
           [0, 1, 1, 3, 0],
           [1, 0, 1, 1, 2],
           [1, 1, 0, 2, 1],
           [1, 1, 0, 1, 2],
           [1, 3, 0, 0, 1],
           [1, 2, 1, 1, 0]])




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




    array([[10.28065478, 20.52686803, 29.66922868],
           [ 9.12998199, 19.84150138, 31.07355396],
           [ 9.35603618, 19.4041181 , 28.91581478],
           [10.32892276, 18.5119342 , 30.65627675]])



### Sampling using `choice`

Works much like the R `sample` function.


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


```python
# Breaks string into list with each letter being an element
list('ABCDEF')
```




    ['A', 'B', 'C', 'D', 'E', 'F']




```python
x
```




    array(['A', 'F', 'D', 'E', 'C', 'B'], dtype='<U1')




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




    array(['E', 'F', 'F'], dtype='<U1')




```python
# Need to specify sampling without replacement if desired
np.random.choice(x, 10)
```




    array(['C', 'E', 'C', 'C', 'D', 'F', 'C', 'F', 'C', 'F'], dtype='<U1')




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

    Cannot take a larger sample than population when 'replace=False'


## Indexing 


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




    array([[ 0,  1,  2,  3,  4],
           [ 5,  6,  7,  8,  9],
           [10, 11, 12, 13, 14],
           [15, 16, 17, 18, 19]])



### Extracing a scalar


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




    6



### Extracting a vector


```python
# Single number returns a row vector
x[1]
```




    array([5, 6, 7, 8, 9])



### Using slices


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




    array([5, 6, 7, 8, 9])




```python
# All rows, column 1
x[:,1]
```




    array([ 1,  6, 11, 16])




```python
# Recall, end of range is exclusive
x[1:3,1:3]
```




    array([[ 6,  7],
           [11, 12]])



### Using slices with strides


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




    array([[ 0,  2,  4],
           [10, 12, 14]])



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

`np.ix_`


```python
# Only words when one of the dimensions returns all (i.e., a colon)
x[:, [0,3]]
```




    array([[ 0,  3],
           [ 5,  8],
           [10, 13],
           [15, 18]])



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]]
```




    array([ 0, 13])



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


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




    array([[ 0,  3],
           [10, 13]])



### 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
```




    [0, 1, 2]




```python
# Original array gets altered!
b = a2[:]
b[1] = 9
a2
```




    array([0, 9, 2])




```python
x
```




    array([[ 0,  1,  2,  3,  4],
           [ 5,  6,  7,  8,  9],
           [10, 11, 12, 13, 14],
           [15, 16, 17, 18, 19]])




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




    array([[ 6,  7,  8],
           [11, 12, 13]])




```python
y *= 10
```


```python
y
```




    array([[ 60,  70,  80],
           [110, 120, 130]])




```python
x
```




    array([[  0,   1,   2,   3,   4],
           [  5,  60,  70,  80,   9],
           [ 10, 110, 120, 130,  14],
           [ 15,  16,  17,  18,  19]])



**Use the copy method to convert a view to a copy**


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


```python
z
```




    array([[ 60,  70,  80],
           [110, 120, 130]])




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


```python
z
```




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




```python
x
```




    array([[  0,   1,   2,   3,   4],
           [  5,  60,  70,  80,   9],
           [ 10, 110, 120, 130,  14],
           [ 15,  16,  17,  18,  19]])



### Boolean indexing


```python
# Always returns 1D array
x[x % 2 == 0]
```




    array([  0,   2,   4,  60,  70,  80,  10, 110, 120, 130,  14,  16,  18])




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




    array([  4,   5,  60,  70,  80,   9,  10, 110, 120, 130,  14,  15,  16,
            17,  18,  19])



### Functions that return indexes


```python
# Gives indices for non-zero entries
idx = np.nonzero(x)
idx
```




    (array([0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]),
     array([1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4]))




```python
x[idx]
```




    array([  1,   2,   3,   4,   5,  60,  70,  80,   9,  10, 110, 120, 130,
            14,  15,  16,  17,  18,  19])




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




    (array([0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]),
     array([4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4]))




```python
x[idx]
```




    array([  4,   5,  60,  70,  80,   9,  10, 110, 120, 130,  14,  15,  16,
            17,  18,  19])



## Universal functions


```python
x
```




    array([[  0,   1,   2,   3,   4],
           [  5,  60,  70,  80,   9],
           [ 10, 110, 120, 130,  14],
           [ 15,  16,  17,  18,  19]])



Operations


```python
x + x
```




    array([[  0,   2,   4,   6,   8],
           [ 10, 120, 140, 160,  18],
           [ 20, 220, 240, 260,  28],
           [ 30,  32,  34,  36,  38]])



Element-wise functions


```python
np.log1p(x) # log(x) + 1
```




    array([[0.        , 0.69314718, 1.09861229, 1.38629436, 1.60943791],
           [1.79175947, 4.11087386, 4.26267988, 4.39444915, 2.30258509],
           [2.39789527, 4.7095302 , 4.79579055, 4.87519732, 2.7080502 ],
           [2.77258872, 2.83321334, 2.89037176, 2.94443898, 2.99573227]])




```python
# Gives lower and upper bound
x.clip(10, 100)
```




    array([[ 10,  10,  10,  10,  10],
           [ 10,  60,  70,  80,  10],
           [ 10, 100, 100, 100,  14],
           [ 15,  16,  17,  18,  19]])



Scans


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




    array([[  0,   1,   3,   6,  10],
           [  5,  65, 135, 215, 224],
           [ 10, 120, 240, 370, 384],
           [ 15,  31,  48,  66,  85]])



Reductions


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




    703




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




    0



## Margins and the `axis` argument


```python
x
```




    array([[  0,   1,   2,   3,   4],
           [  5,  60,  70,  80,   9],
           [ 10, 110, 120, 130,  14],
           [ 15,  16,  17,  18,  19]])



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


```python
x.shape
```




    (4, 5)




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




    35.15



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


```python
# Marginalizes 0th axis, which refers to rows... means rows will "disappear"
x.mean(axis=0)
```




    array([ 7.5 , 46.75, 52.25, 57.75, 11.5 ])



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


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




    array([ 2. , 44.8, 76.8, 17. ])



Note marginalizing out the last axis is a common default.


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




    array([ 2. , 44.8, 76.8, 17. ])



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


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




    array([[[0.50405392, 0.57784173, 0.41159044, 0.72692582],
            [0.47458899, 0.02542905, 0.80778833, 0.87714966],
            [0.58456692, 0.35611667, 0.70781391, 0.44196183]],
    
           [[0.08716115, 0.32468677, 0.82071695, 0.52673336],
            [0.29664671, 0.3067223 , 0.6104232 , 0.61916896],
            [0.20120738, 0.40392386, 0.57092189, 0.29781173]]])




```python
x.shape
```




    (2, 3, 4)




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




    (3, 4)




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




    (2, 4)




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




    (2, 3)




```python
# Can marginalize out multiple axes
x.mean(axis=(0,1)).shape
```




    (4,)




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




    (3,)




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




    (2,)



## 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
# When comparing dimensions, start at right-most dimension and work to the left (checking rules above)
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
```




    (1, 2)




```python
x + y
```




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



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
# Inserts a dummy dimension
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
```




    (2, 2)




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

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


### More examples of broadcasting


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


```python
x1
```




    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])




```python
x1 * 10
```




    array([  0,  10,  20,  30,  40,  50,  60,  70,  80,  90, 100, 110])




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


```python
x2
```




    array([[0, 4, 1, 1],
           [3, 7, 0, 5],
           [5, 2, 0, 0]])




```python
x2 * 10
```




    array([[ 0, 40, 10, 10],
           [30, 70,  0, 50],
           [50, 20,  0,  0]])




```python
x2.shape
```




    (3, 4)



### Column-wise broadcasting


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




    (4,)




```python
x2 - mu
```




    array([[-2.66666667, -0.33333333,  0.66666667, -1.        ],
           [ 0.33333333,  2.66666667, -0.33333333,  3.        ],
           [ 2.33333333, -2.33333333, -0.33333333, -2.        ]])




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




    array([1.48029737e-16, 2.96059473e-16, 3.70074342e-17, 0.00000000e+00])



### Row wise broadcasting


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




    (3,)




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

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


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


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




    (3, 1)




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




    array([[-1.5 ,  2.5 , -0.5 , -0.5 ],
           [-0.75,  3.25, -3.75,  1.25],
           [ 3.25,  0.25, -1.75, -1.75]])




```python
# Does same thing as the "None" syntax above
x2 - mu[:, np.newaxis]
```




    array([[-1.5 ,  2.5 , -0.5 , -0.5 ],
           [-0.75,  3.25, -3.75,  1.25],
           [ 3.25,  0.25, -1.75, -1.75]])




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




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



#### Reshaping works too


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




    array([[-1.5 ,  2.5 , -0.5 , -0.5 ],
           [-0.75,  3.25, -3.75,  1.25],
           [ 3.25,  0.25, -1.75, -1.75]])



#### Exercise in broadcasting

Creating a 12 by 12 multiplication table


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




    array([[  1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12],
           [  2,   4,   6,   8,  10,  12,  14,  16,  18,  20,  22,  24],
           [  3,   6,   9,  12,  15,  18,  21,  24,  27,  30,  33,  36],
           [  4,   8,  12,  16,  20,  24,  28,  32,  36,  40,  44,  48],
           [  5,  10,  15,  20,  25,  30,  35,  40,  45,  50,  55,  60],
           [  6,  12,  18,  24,  30,  36,  42,  48,  54,  60,  66,  72],
           [  7,  14,  21,  28,  35,  42,  49,  56,  63,  70,  77,  84],
           [  8,  16,  24,  32,  40,  48,  56,  64,  72,  80,  88,  96],
           [  9,  18,  27,  36,  45,  54,  63,  72,  81,  90,  99, 108],
           [ 10,  20,  30,  40,  50,  60,  70,  80,  90, 100, 110, 120],
           [ 11,  22,  33,  44,  55,  66,  77,  88,  99, 110, 121, 132],
           [ 12,  24,  36,  48,  60,  72,  84,  96, 108, 120, 132, 144]])



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


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




    array([[14.17567293,  6.50796953, 11.12009193,  2.62072251],
           [ 8.59736853,  5.25756209, 15.14943199, 11.67007694],
           [ 0.65966542,  9.4391506 ,  7.21314915, 10.57761823]])



Scaling column-wise


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




    array([[ 1.14766576, -0.31969229, -0.01259191, -1.40554016],
           [ 0.14181192, -1.03319519,  1.23099228,  0.8382051 ],
           [-1.28947768,  1.35288748, -1.21840037,  0.56733506]])



Scaling row-wise


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




    array([[ 1.26477055, -0.47645994,  0.57088993, -1.35920053],
           [-0.42902964, -1.34096844,  1.36002041,  0.40997767],
           [-1.64381989,  0.64233703,  0.06269158,  0.93879128]])



## 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
```




    array([[ 0,  1,  2,  3,  4],
           [ 5,  6,  7,  8,  9],
           [10, 11, 12, 13, 14],
           [15, 16, 17, 18, 19]])




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




    array([[ True, False,  True, False,  True],
           [False,  True, False,  True, False],
           [ True, False,  True, False,  True],
           [False,  True, False,  True, False]])



- 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
```




    masked_array(
      data=[[--, 1, --, 3, --],
            [5, --, 7, --, 9],
            [--, 11, --, 13, --],
            [15, --, 17, --, 19]],
      mask=[[ True, False,  True, False,  True],
            [False,  True, False,  True, False],
            [ True, False,  True, False,  True],
            [False,  True, False,  True, False]],
      fill_value=999999)




```python
m.data
```




    array([[ 0,  1,  2,  3,  4],
           [ 5,  6,  7,  8,  9],
           [10, 11, 12, 13, 14],
           [15, 16, 17, 18, 19]])




```python
m.mask
```




    array([[ True, False,  True, False,  True],
           [False,  True, False,  True, False],
           [ True, False,  True, False,  True],
           [False,  True, False,  True, False]])




```python
# Only operates on the un-masked values
m.sum(axis=0).data
```




    array([20, 12, 24, 16, 28])




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




    array([ 4, 21, 24, 51])



### Often used with missing value sentinels


```python
import warnings

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


```python
x1
```




    array([[ 0., inf,  2., inf,  4.],
           [inf,  6., inf,  8., inf],
           [10., inf, 12., inf, 14.],
           [inf, 16., inf, 18., inf]])




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




    inf




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




    masked_array(
      data=[[0.0, --, 2.0, --, 4.0],
            [--, 6.0, --, 8.0, --],
            [10.0, --, 12.0, --, 14.0],
            [--, 16.0, --, 18.0, --]],
      mask=[[False,  True, False,  True, False],
            [ True, False,  True, False,  True],
            [False,  True, False,  True, False],
            [ True, False,  True, False,  True]],
      fill_value=1e+20)




```python
x2.data
```




    array([[ 0., inf,  2., inf,  4.],
           [inf,  6., inf,  8., inf],
           [10., inf, 12., inf, 14.],
           [inf, 16., inf, 18., inf]])




```python
x2.mask
```




    array([[False,  True, False,  True, False],
           [ True, False,  True, False,  True],
           [False,  True, False,  True, False],
           [ True, False,  True, False,  True]])




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




    90.0




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




    array([[ 0.,  0.,  2.,  0.,  4.],
           [ 0.,  6.,  0.,  8.,  0.],
           [10.,  0., 12.,  0., 14.],
           [ 0., 16.,  0., 18.,  0.]])



## Combining `ndarray`s


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


```python
x1
```




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




```python
x2
```




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




```python
x3
```




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



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


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




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



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


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




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



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


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




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



### Stacking


```python
# Same as row-bind
np.vstack([x1, x3])
```




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




```python
# Same as column-bind
np.hstack([x1, x2])
```




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




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




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



### Generic stack with axis argument


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




    array([[[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.]],
    
           [[3., 3., 3., 3., 3.],
            [3., 3., 3., 3., 3.],
            [3., 3., 3., 3., 3.]]])




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




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




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




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



### Repetition and tiling

#### For a vector


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


```python
np.repeat(x, 3) # Element by element
```




    array([1, 1, 1, 2, 2, 2, 3, 3, 3])




```python
np.tile(x, 3) # Whole array at a time
```




    array([1, 2, 3, 1, 2, 3, 1, 2, 3])



#### For a matrix


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




    array([[0, 1, 2],
           [3, 4, 5]])




```python
np.repeat(x, 3) # With no axis specification, creates 1D array
```




    array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5])




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




    array([[0, 1, 2],
           [0, 1, 2],
           [0, 1, 2],
           [3, 4, 5],
           [3, 4, 5],
           [3, 4, 5]])




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




    array([[0, 0, 0, 1, 1, 1, 2, 2, 2],
           [3, 3, 3, 4, 4, 4, 5, 5, 5]])




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




    array([[0, 1, 2, 0, 1, 2],
           [3, 4, 5, 3, 4, 5],
           [0, 1, 2, 0, 1, 2],
           [3, 4, 5, 3, 4, 5],
           [0, 1, 2, 0, 1, 2],
           [3, 4, 5, 3, 4, 5]])



## Splitting `ndarray`s


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


```python
x
```




    array([[ 0,  1,  2,  3,  4,  5,  6,  7],
           [ 8,  9, 10, 11, 12, 13, 14, 15],
           [16, 17, 18, 19, 20, 21, 22, 23],
           [24, 25, 26, 27, 28, 29, 30, 31]])




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




    [array([[0, 1, 2, 3, 4, 5, 6, 7]]),
     array([[ 8,  9, 10, 11, 12, 13, 14, 15]]),
     array([[16, 17, 18, 19, 20, 21, 22, 23]]),
     array([[24, 25, 26, 27, 28, 29, 30, 31]])]




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




    [array([[ 0,  1],
            [ 8,  9],
            [16, 17],
            [24, 25]]), array([[ 2,  3],
            [10, 11],
            [18, 19],
            [26, 27]]), array([[ 4,  5],
            [12, 13],
            [20, 21],
            [28, 29]]), array([[ 6,  7],
            [14, 15],
            [22, 23],
            [30, 31]])]



## 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')
```




    array([[ 0,  1,  2,  3],
           [ 4,  5,  6,  7],
           [ 8,  9, 10, 11],
           [12, 13, 14, 15]])




```python
np.savez('xy.npz', x=x, y=y) # Save multiple numpy arrays into a single file (with .npz file name ending)
```


```python
arr = np.load('xy.npz') # Single structure, which behaves as a dictionary
```


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




    array([[ 0,  1,  2,  3],
           [ 4,  5,  6,  7],
           [ 8,  9, 10, 11],
           [12, 13, 14, 15]])




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




    array([[ 0,  1,  2,  3],
           [ 4,  5,  6,  7],
           [ 8,  9, 10, 11],
           [12, 13, 14, 15],
           [16, 17, 18, 19]])




```python
import os
```


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


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

### Einstein summation notation

Useful for higher order tensors. Check out blog post on Einstein summation notation in most up-to-date version of notebook


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




    array([[1, 2, 3],
           [4, 5, 6],
           [7, 8, 9]])




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




    15




```python
np.einsum('ji->ij', x) # input -> output
```




    array([[1, 4, 7],
           [2, 5, 8],
           [3, 6, 9]])




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




    array([1, 5, 9])




```python
np.einsum('ij->i', x) # Returns row sums
```




    array([ 6, 15, 24])




```python
np.einsum('ij->j', x) # Returns column sums
```




    array([12, 15, 18])




```python
np.einsum('mn, np -> mp', x, x) # Standard matrix multiplication
```




    array([[ 30,  36,  42],
           [ 66,  81,  96],
           [102, 126, 150]])




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


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

## Vectorization

When at all possible, try to avoid using loops

### 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
```




    array([ 0,  1,  4,  9, 16, 25, 36, 49, 64, 81])




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




    array([ 0,  1,  4,  9, 16, 25, 36, 49, 64, 81])



Using ufuncs


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




    array([0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])




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




    array([0.        , 0.69314718, 1.60943791, 2.30258509, 2.83321334,
           3.25809654, 3.61091791, 3.91202301, 4.17438727, 4.40671925])



### 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
```




    2.2348383265132834




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




    2.2348383265132834




```python
xs @ ys # Matrix multiplication
```




    2.2348383265132834




```python
# Outer product
xs.reshape(10, 1) @ xs.reshape(1, 10)
np.outer(xs, ys)
```




    array([[0.29821288, 0.15492211, 0.43784934, 0.28829146, 0.13613758,
            0.41864202, 0.23636405, 0.0583262 , 0.44300641, 0.21491643],
           [0.41046199, 0.21323571, 0.60265845, 0.39680609, 0.18738058,
            0.57622137, 0.32533289, 0.08028053, 0.60975667, 0.29581226],
           [0.28347518, 0.14726585, 0.41621081, 0.27404409, 0.12940966,
            0.39795271, 0.22468293, 0.05544371, 0.42111302, 0.20429525],
           [0.16905646, 0.08782513, 0.24821618, 0.16343203, 0.07717621,
            0.23732757, 0.13399445, 0.03306504, 0.25113972, 0.12183582],
           [0.21045376, 0.10933109, 0.30899752, 0.20345205, 0.09607454,
            0.29544259, 0.16680602, 0.04116176, 0.31263695, 0.15167008],
           [0.14763515, 0.07669672, 0.21676446, 0.14272339, 0.06739713,
            0.20725557, 0.11701588, 0.02887534, 0.21931756, 0.10639789],
           [0.47104689, 0.24470967, 0.69161188, 0.45537535, 0.21503828,
            0.66127265, 0.37335258, 0.09213007, 0.69975781, 0.33947466],
           [0.39955291, 0.20756843, 0.58664126, 0.38625996, 0.18240047,
            0.56090682, 0.31668634, 0.07814687, 0.59355084, 0.2879503 ],
           [0.10378559, 0.05391679, 0.15238259, 0.10033268, 0.04737931,
            0.14569796, 0.08226064, 0.02029898, 0.15417738, 0.07479633],
           [0.32571953, 0.16921186, 0.47823583, 0.31488298, 0.14869469,
            0.45725686, 0.25816587, 0.06370611, 0.48386858, 0.23473996]])




```python
sum(xs * ys)
```




    2.2348383265132834




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




    2.2348383265132834




```python
np.einsum('i,i->', xs, ys) # Have two vectors and ask for a scalar... only have one option (dot product)
```




    2.2348383265132834



### 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
```




    array([0.75908805])




```python
betas
```




    array([[0.99604517],
           [0.2925755 ]])




```python
xs
```




    array([[0.01821893, 0.50077046],
           [0.88464832, 0.12055759],
           [0.56973445, 0.18519832]])



### Using loops


```python
%%timeit

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
```

    30.8 µs ± 639 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)


### Removing inner loop


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

    9.85 µs ± 49.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)


### Removing all loops


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

    2.65 µs ± 32 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)


### 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
```




    array([[1.        , 0.01821893, 0.50077046],
           [1.        , 0.88464832, 0.12055759],
           [1.        , 0.56973445, 0.18519832]])




```python
alpha
```




    array([0.75908805])




```python
betas
```




    array([[0.99604517],
           [0.2925755 ]])




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




    array([[0.75908805],
           [0.99604517],
           [0.2925755 ]])




```python
%%timeit
ys = X @ betas_
ys
```

    1.4 µs ± 18.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)



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

    3.69 µs ± 23.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)


### 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') # Create matrix with border of one along each side
x[(w//2-1):(w//2+2), (h//2-1):(h//2+2)] = 1 # Exact center is one

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') # Plot everything except for boundaries
    
    # When we use the slice notation, we are using views, not making a copy!
    center = x[1:-1, 1:-1] # Everything except for border
    top = x[:-2, 1:-1] # Offset by one to the top
    bottom = x[2:, 1:-1] # Offset by one to the bottom
    left = x[1:-1, :-2] # Offset by to the left
    right = x[1:-1, 2:] # Offset by to the right
    nbrs = np.dstack([center, left, right, bottom, top]) 
    x = np.sum(wts * nbrs, axis=-1) # Sum across depth to return 102 x 102 matrix which is result of averaging with neighbors
```


```python

```

2f040e7cbbed06b6fa29ff4c526ae831289256cd
ipynb
Jupyter Notebook
notebooks/copies/lectures/S03_Numpy.ipynb
robkravec/sta-663-2021
4dc8018f7b172eaf81da9edc33174768ff157939
MIT
notebooks/copies/lectures/S03_Numpy.ipynb
robkravec/sta-663-2021
4dc8018f7b172eaf81da9edc33174768ff157939
MIT
notebooks/copies/lectures/S03_Numpy.ipynb
robkravec/sta-663-2021
4dc8018f7b172eaf81da9edc33174768ff157939
MIT
Qwen/Qwen-72B
1. YES
2. YES