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)

"""
Explanation: 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
NumPy: creating and manipulating numerical data
Advanced Numpy
100 Numpy Exercises

NDArray
The 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
End of explanation
"""


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

"""
Explanation: Array creation
End of explanation
"""


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)

"""
Explanation: Array manipulation
End of explanation
"""


x

x[0]

x[0,:]

x[:,0]

x[-1]

x[1,1]

x[:, 1:3]

"""
Explanation: Array indexing
End of explanation
"""


x >= 2

x[x > 2]

"""
Explanation: Boolean indexing
End of explanation
"""


x[0, [1,2]]

"""
Explanation: Fancy indexing
End of explanation
"""


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)

"""
Explanation: Calculations and broadcasting
Broadcasting refers to the set of rules that numpy uses to perfrom operations on arrays with different shapes. See official documentation 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.
End of explanation
"""


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

"""
Explanation: Combining and splitting arrays
End of explanation
"""


y

y.sum()

y.sum(0) # column sum

y.sum(1) # row sum

"""
Explanation: Reductions
End of explanation
"""


z = (y - y.mean(0))/y.std(0)

z

z.mean(0), z.std(0)

"""
Explanation: Standardize by column mean and standard deviation
End of explanation
"""


z = (y - y.mean(1)[:,None])/y.std(1)[:,None]

z

z.mean(1), z.std(1)

"""
Explanation: Standardize by row mean and standard deviation
End of explanation
"""


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

"""
Explanation: 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)
End of explanation
"""


pts[None, :].shape

pts[:, None].shape

m = pts[None, :] - pts[:, None]
m

m**2

(m**2).shape

"""
Explanation: Using broadcasting
End of explanation
"""


np.sum((pts[None, :] - pts[:, None])**2, -1)

"""
Explanation: 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.
End of explanation
"""


np.sqrt(np.sum((pts[None, :] - pts[:, None])**2, -1))

"""
Explanation: Basically, the distance matrix can be calculated in one line of numpy code
End of explanation
"""


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

"""
Explanation: Let's put them in functions and compare the time.
End of explanation
"""


np.alltrue(pdist1(pts) == pdist2(pts))

pts = np.random.random((1000, 2))

%timeit pdist1(pts)

%timeit pdist2(pts)

"""
Explanation: Check that the outputs are the same
End of explanation
"""


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)

"""
Explanation: But don't give up on loops yet
End of explanation
"""


N = 5

np.tri(N)

np.tri(N, N-1)

np.tri(N, N-1, -1)

"""
Explanation: 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
End of explanation
"""


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]

"""
Explanation: Use broadcasting to create a new index matrix
End of explanation
"""


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)

"""
Explanation: All but one
R 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.
End of explanation
"""


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);

"""
Explanation: 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.
End of explanation
"""


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

"""
Explanation: Generalized ufuncs
A 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.
End of explanation
"""


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

"""
Explanation: Saving and loading NDArrays
Saving to and loading from text files
End of explanation
"""


np.save('../data/x1.npy', x1)

!cat ../data/x1.npy

x3 = np.load('../data/x1.npy')
x3

"""
Explanation: Saving to and loading from binary files (much faster and also preserves dtype)
End of explanation
"""


%load_ext version_information
%version_information numpy, numba, matplotlib

"""
Explanation: Version information
End of explanation
"""


