<jupyter_start><jupyter_text>Using `numpy`====The foundation for numerical computation in Python is the `numpy` package, and essentially all scientific libraries in Python build on this - e.g. `scipy`, `pandas`, `statsmodels`, `scikit-learn`, `cv2` etc. The basic data structure in `numpy` is the NDArray, and it is essential to become familiar with how to slice and dice this object.Numpy also has the `random`, and `linalg` modules that we will discuss in later lectures. Resources----- [Numpy for R users](http://mathesaurus.sourceforge.net/r-numpy.html)- [NumPy: creating and manipulating numerical data](http://www.scipy-lectures.org/intro/numpy/index.html)- [Advanced Numpy](http://www.scipy-lectures.org/advanced/advanced_numpy/index.html)- [100 Numpy Exercises](http://www.labri.fr/perso/nrougier/teaching/numpy.100/)  NDArrayThe base structure in `numpy` is `ndarray`, used to represent vectors, matrices and higher-dimensional arrays. Each `ndarray` has the following attributes:- dtype = corresponds to data types in C- shape = dimensions of array- strides = number of bytes to step in each direction when traversing the array<jupyter_code>x = np.array([1,2,3,4,5,6])
print(x)
print('dytpe', x.dtype)
print('shape', x.shape)
print('strides', x.strides)
x.shape = (2,3)
print(x)
print('dytpe', x.dtype)
print('shape', x.shape)
print('strides', x.strides)
x = x.astype('complex')
print(x)
print('dytpe', x.dtype)
print('shape', x.shape)
print('strides', x.strides)<jupyter_output>[[ 1.+0.j  2.+0.j  3.+0.j]
 [ 4.+0.j  5.+0.j  6.+0.j]]
dytpe complex128
shape (2, 3)
strides (48, 16)<jupyter_text>Array creation----<jupyter_code>np.array([1,2,3])
np.array([1,2,3], np.float64)
np.arange(3)
np.arange(3, 6, 0.5)
np.array([[1,2,3],[4,5,6]])
np.ones(3)
np.zeros((3,4))
np.eye(4)
np.diag([1,2,3,4])
np.fromfunction(lambda i, j: i**2+j**2, (4,5))<jupyter_output><empty_output><jupyter_text>Array manipulation----<jupyter_code>x = np.fromfunction(lambda i, j: i**2+j**2, (4,5))
x
x.shape
x.size
x.dtype
x.astype(np.int64)
x.T
x.reshape(2,-1)<jupyter_output><empty_output><jupyter_text>Array indexing----<jupyter_code>x
x[0]
x[0,:]
x[:,0]
x[-1]
x[1,1]
x[:, 1:3]<jupyter_output><empty_output><jupyter_text>Boolean indexing<jupyter_code>x >= 2
x[x > 2]<jupyter_output><empty_output><jupyter_text>Fancy indexing<jupyter_code>x[0, [1,2]]<jupyter_output><empty_output><jupyter_text>Calculations and broadcasting----Broadcasting refers to the set of rules that numpy uses to perfrom operations on arrays with different shapes. See official [documentation](http://docs.scipy.org/doc/numpy-1.10.1/user/basics.broadcasting.html) for a clear explanation of the rules. Array shapes can be manipulated using the `reshape` method or by inserting a new axis with `np.newaxis`. Note that `np.newaxis` is an alias for `None`, which I sometimes use in my examples.<jupyter_code>x = np.fromfunction(lambda i, j: i**2+j**2, (2,3))
x
x * 5
x + x
x @ x.T
x.T @ x
np.log1p(x)
np.exp(x)<jupyter_output><empty_output><jupyter_text>Combining and splitting arrays----<jupyter_code>x
np.r_[x, x]
np.vstack([x, x])
np.concatenate([x, x], axis=0)
np.c_[x,x]
np.hstack([x, x])
np.concatenate([x,x], axis=1)
y = np.r_[x, x]
y
a, b, c = np.hsplit(y, 3)
a
b
c
np.vsplit(y, [3])
np.split(y, [3], axis=0)
np.hstack(np.hsplit(y, 3))<jupyter_output><empty_output><jupyter_text>Reductions----<jupyter_code>y
y.sum()
y.sum(0) # column sum
y.sum(1) # row sum<jupyter_output><empty_output><jupyter_text>Standardize by column mean and standard deviation<jupyter_code>z = (y - y.mean(0))/y.std(0)
z
z.mean(0), z.std(0)<jupyter_output><empty_output><jupyter_text>Standardize by row mean and standard deviation<jupyter_code>z = (y - y.mean(1)[:,None])/y.std(1)[:,None]
z
z.mean(1), z.std(1)<jupyter_output><empty_output><jupyter_text>Example: Calculating pairwise distance matrix using broadcasting and vectorization----Calculate the pairwise distance matrix between the following points- (0,0)- (4,0)- (4,3)- (0,3)<jupyter_code>def distance_matrix_py(pts):
    """Returns matrix of pairwise Euclidean distances. Pure Python version."""
    n = len(pts)
    p = len(pts[0])
    m = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            s = 0
            for k in range(p):
                s += (pts[i,k] - pts[j,k])**2
            m[i, j] = s**0.5
    return m
def distance_matrix_np(pts):
    """Returns matrix of pairwise Euclidean distances. Vectorized numpy version."""
    return np.sum((pts[None,:] - pts[:, None])**2, -1)**0.5
pts = np.array([(0,0), (4,0), (4,3), (0,3)])
pts
pts.shape
n = pts.shape[0]
p = pts.shape[1]
dist = np.zeros((n, n))
for i in range(n):
    for j in range(n):
        s = 0
        for k in range(p):
            s += (pts[i, k] - pts[j, k])**2
        dist[i, j] = np.sqrt(s)
dist<jupyter_output><empty_output><jupyter_text>Using broadcasting<jupyter_code>pts[None, :].shape
pts[:, None].shape
m = pts[None, :] - pts[:, None]
m
m**2
(m**2).shape<jupyter_output><empty_output><jupyter_text>We want to end up with a 4 by 4 matrix, so sum over the axis with dimension 2. This is axis=2, or axis=-1 since it is the first axis from the end.<jupyter_code>np.sum((pts[None, :] - pts[:, None])**2, -1)<jupyter_output><empty_output><jupyter_text>Basically, the distance matrix can be calculated in one line of numpy code<jupyter_code>np.sqrt(np.sum((pts[None, :] - pts[:, None])**2, -1))<jupyter_output><empty_output><jupyter_text>Let's put them in functions and compare the time.<jupyter_code>def pdist1(pts):
    n = pts.shape[0]
    p = pts.shape[1]
    dist = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            s = 0
            for k in range(p):
                s += (pts[i, k] - pts[j, k])**2
            dist[i, j] = s
    return np.sqrt(dist)
def pdist2(pts):
    return np.sqrt(np.sum((pts[None, :] - pts[:, None])**2, -1))<jupyter_output><empty_output><jupyter_text>Check that the outputs are the same<jupyter_code>np.alltrue(pdist1(pts) == pdist2(pts))
pts = np.random.random((1000, 2))
%timeit pdist1(pts)
%timeit pdist2(pts)<jupyter_output>10 loops, best of 3: 77.3 ms per loop<jupyter_text>But don't give up on loops yet<jupyter_code>from numba import njit
@njit
def pdist3(pts):
    n = pts.shape[0]
    p = pts.shape[1]
    dist = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            s = 0
            for k in range(p):
                s += (pts[i, k] - pts[j, k])**2
            dist[i, j] = s
    return np.sqrt(dist)
%timeit pdist3(pts)<jupyter_output>The slowest run took 27.17 times longer than the fastest. This could mean that an intermediate result is being cached 
1 loops, best of 3: 16.1 ms per loop<jupyter_text>What is going on?This is 3-5 times faster than the broadcasting version! We have just performed Just In Time (JIT) compilation of a function, which will be discussed in a later lecture. Example: Consructing leave-one-out arrays----Another example of numpy trickery is to construct a leave-one-out matrix of a vector of length k. In the matrix, each row is a vector of length k-1, with a different vector component dropped each time. This can be used for LOOCV to evalaute the out-of-sample accuracy of a predictive model.For example, suppose you have data points [(1,4), (2,7), (3,11), (4,9), (5,15)] that you want to perfrom LOOCV on for a simple regression model. For each cross-validation, you use one point for testing, and the remaining 4 points for training. In other words, you want the training set to be:```[(2,7), (3,11), (4,9), (5,15)][(1,4), (3,11), (4,9), (5,15)][(1,4), (2,7),  (4,9), (5,15)][(1,4), (2,7), (3,11), (5,15)][(1,4), (2,7), (3,11), (4,9)]```Here is one way to do create the training set using numpy tricks. Create a triangular matrix with N rows, N-1 columns and offset from diagnonal by -1<jupyter_code>N = 5
np.tri(N)
np.tri(N, N-1)
np.tri(N, N-1, -1)<jupyter_output><empty_output><jupyter_text>Use broadcasting to create a new index matrix<jupyter_code>np.arange(1, N)
np.arange(1, N) - np.tri(N, N-1, -1)
idx = np.arange(1, N) - np.tri(N, N-1, -1).astype('int')
data = np.array([(1,4), (2,7), (3,11), (4,9), (5,15)])
data
data[idx]<jupyter_output><empty_output><jupyter_text>All but oneR uses negative indexing to mean delete the component at that index. Because Python uses negative indexing to mean count from the end, we have to do a little more work to get the same effect. Here are two ways of deleting one item from a vector.<jupyter_code>def f1(a, k):
    idx = np.ones_like(a).astype('bool')
    idx[k] = 0
    return a[idx]
def f2(a, k):
    return np.r_[a[:k], a[k+1:]]
a = np.arange(100)
k = 50
%timeit f1(a, k)
%timeit f2(a, k)<jupyter_output>10000 loops, best of 3: 47.6 µs per loop<jupyter_text>Universal functions (Ufuncs)Functions that work on both scalars and arrays are known as ufuncs. For arrays, ufuncs apply the function in an element-wise fashion. Use of ufuncs is an esssential aspect of vectorization and typically much more computationally efficient than using an explicit loop over each element.<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline

xs = np.linspace(0, 2*np.pi, 100)
ys = np.sin(xs) # np.sin is a universal function
plt.plot(xs, ys);<jupyter_output><empty_output><jupyter_text>Generalized ufuncsA universal function performs vectorized looping over scalars. A generalized ufunc performs looping over vectors or arrays. Currently, numpy only ships with a single generalized ufunc. However, they play an important role for JIT compilation with `numba`, a topic we will cover in future lectures.<jupyter_code>from numpy.core.umath_tests import matrix_multiply

print(matrix_multiply.signature)
us = np.random.random((5, 2, 3)) # 5 2x3 matrics
vs = np.random.random((5, 3, 4)) # 5 3x4 matrices
us
vs
# perform matrix multiplication for each of the 5 sets of matrices
ws = matrix_multiply(us, vs) 
ws.shape
ws<jupyter_output><empty_output><jupyter_text>Saving and loading NDArrays----  Saving to and loading from text files<jupyter_code>x1 = np.arange(1,10).reshape(3,3)
x1
np.savetxt('../data/x1.txt', x1)
!cat ../data/x1.txt
x2 = np.loadtxt('../data/x1.txt')
x2<jupyter_output><empty_output><jupyter_text>Saving to and loading from binary files (much faster and also preserves dtype)<jupyter_code>np.save('../data/x1.npy', x1)
!cat ../data/x1.npy
x3 = np.load('../data/x1.npy')
x3<jupyter_output><empty_output><jupyter_text>Version information----<jupyter_code>%load_ext version_information
%version_information numpy, numba, matplotlib<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1523222748815166, 'fraction_numerical': 0.03213270142180095, 'mean_word_length': 4.529874213836478, 'pattern_counts': {'":': 0, '<': 101, '<?xml version=': 0, '>': 103, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '240050', 'license': 'MIT', 'n_tokens_mistral': 4129, 'n_tokens_neox': 3970, 'n_words': 1145, 'path': 'notebook/05_Numbers.ipynb', 'repo_name': 'cliburn/sta-663-2017', 'size': '10550'}