import numpy as np
np.set_printoptions(precision=3) # Keeps output from getting cluttered.


x = np.random.rand(10)

sum_x = 0

for i in range(x.shape[0]):
    sum_x += x[i]

print(sum_x)

# Write the dot_loop function

def dot_loop(x1, x2):
    """ Returns the dot product of the 1-d numpy arrays x1 and x2 """

    # YOUR CODE HERE
    raise NotImplementedError()

# Tests for the dot_loop function.

import numpy as np

a = np.array([1.0, 2.0, 3.0])
b = np.array([-1.0, 0.5, 4.0])

np.testing.assert_allclose(dot_loop(a, b), 12.0)
np.testing.assert_allclose(dot_loop(b, a), 12.0)   # the dot product is symmetric
np.testing.assert_allclose(dot_loop(a, a), 14.0)
np.testing.assert_allclose(dot_loop(np.array([1.0, 0.0]),
                                   np.array([0.0, 1.0])), 0.0)  # orthogonal

# The shape property returns a tuple containing the shape of a numpy array,
# and tuples can be unpacked into separate variables:

A = np.random.random((5, 3))
print(A.shape)

rows, cols = A.shape
print("rows: {}, cols: {}".format(rows, cols))

# np.zeros creates an array of a given shape filled with zeros.  This is the
# usual way to allocate an array that you are about to fill in element by
# element.  Note the DOUBLE parentheses: the shape is a single tuple argument.

C = np.zeros((2, 3))
print(C)

# Individual entries of a 2-d array are accessed with comma notation.
# Use C[i, j].  Do NOT use C[i][j].

C[1, 2] = 7.0
print(C)

# ##
# Exercise: Complete matmul_loop
# ##

def matmul_loop(A, B):
    """ Return the matrix product AB using nested for loops.
         Arguments:
          A - An m x d numpy array
          B - A  d x k numpy array
         Returns:
          An m x k numpy array
    """
    # YOUR CODE HERE
    raise NotImplementedError()

# Tests for matmul_loop

B = np.array([[1.0, 2.0, -2.0],
              [3.0, 0.0, 2.0]])

A = np.array([[2.0, -1.0, 0.0],
              [1.0, 3.0, -2.0],
              [0.0, 1.0, 2.0]])

# This is the product you computed by hand above:
expected_BA = np.array([[4.0, 3.0, -8.0],
                        [6.0, -1.0, 4.0]])

np.testing.assert_allclose(matmul_loop(B, A), expected_BA)
assert matmul_loop(B, A).shape == (2, 3)

# Non-square shapes need to work too:
expected_ABt = np.array([[0.0, 6.0],
                         [11.0, -1.0],
                         [-2.0, 4.0]])

np.testing.assert_allclose(matmul_loop(A, B.T), expected_ABt)

# A matrix times a single column vector:
y = np.array([[2.0], [-4.0], [1.0]])
np.testing.assert_allclose(matmul_loop(A, y), np.array([[8.0], [-12.0], [-2.0]]))

# PREDICT BEFORE YOU RUN THIS CELL.  What will the result be?

P = np.array([[1.0, 2.0],
              [3.0, 4.0]])

Q = np.array([[0.0, 1.0],
              [1.0, 0.0]])

# The * operator multiplies corresponding entries.  It does not compute
# a matrix product.
print("P * Q (elementwise):\n{}".format(P * Q))

# np.sum reduces an array to a single number by adding up all of its entries:

print(np.sum(np.array([1.0, 2.0, 3.0])))

def dot(x1, x2):
    """ Returns the dot product of the 1-d numpy arrays x1 and x2 """
    # Note: This should be a one-liner!
    # YOUR CODE HERE
    raise NotImplementedError()

# Tests for dot

a = np.array([1.0, 2.0, 3.0])
b = np.array([-1.0, 0.5, 4.0])

np.testing.assert_allclose(dot(a, b), 12.0)
np.testing.assert_allclose(dot(b, a), 12.0)
np.testing.assert_allclose(dot(a, a), 14.0)
np.testing.assert_allclose(dot(np.array([1.0, 0.0]),
                              np.array([0.0, 1.0])), 0.0)

# TEST SPEED
dim = 10000 # Try changing this value to see how it impacts the running times.
a = np.random.random(dim)
b = np.random.random(dim)

print('time for dot_loop')
%timeit dot_loop(a, b)

print('--------\n\nTime for dot without loops')
%timeit dot(a, b)

a = np.array([1.0, 2.0, 3.0])
b = np.array([-1.0, 0.5, 4.0])

print(dot(a, b))
print(np.dot(a, b))

# Slicing can be used to extract a single row or column from a numpy array.
# A[i, :] is a row and A[:, j] is a column.  Both come back as 1-d arrays:

A = np.random.random((5, 3))
print(A)
print("\nrow 0:            {}".format(A[0, :]))
print("middle column:    {}".format(A[:, 1]))

def matmul_dot(A, B):
    """ Return the matrix product AB using two loops and your dot function.
         Arguments:
          A - An m x d numpy array
          B - A  d x k numpy array
         Returns:
          An m x k numpy array
    """
    # YOUR CODE HERE
    raise NotImplementedError()

# Tests for matmul_dot

B = np.array([[1.0, 2.0, -2.0],
              [3.0, 0.0, 2.0]])

A = np.array([[2.0, -1.0, 0.0],
              [1.0, 3.0, -2.0],
              [0.0, 1.0, 2.0]])

expected_BA = np.array([[4.0, 3.0, -8.0],
                        [6.0, -1.0, 4.0]])

np.testing.assert_allclose(matmul_dot(B, A), expected_BA)
assert matmul_dot(B, A).shape == (2, 3)

expected_ABt = np.array([[0.0, 6.0],
                         [11.0, -1.0],
                         [-2.0, 4.0]])

np.testing.assert_allclose(matmul_dot(A, B.T), expected_ABt)

y = np.array([[2.0], [-4.0], [1.0]])
np.testing.assert_allclose(matmul_dot(A, y), np.array([[8.0], [-12.0], [-2.0]]))

# Your two implementations should agree on arbitrary inputs:
R1 = np.random.random((7, 5))
R2 = np.random.random((5, 4))
np.testing.assert_allclose(matmul_dot(R1, R2), matmul_loop(R1, R2))

A_small = np.array([[1.0, 2.0, 0.0],
                    [4.0, -1.0, 3.0]])
B_small = np.array([[2.0, 1.0],
                    [0.0, -3.0],
                    [5.0, 2.0]])

print("matmul_dot(A_small, B_small):\n{}".format(matmul_dot(A_small, B_small)))
print("\nA_small @ B_small:\n{}".format(A_small @ B_small))
np.testing.assert_allclose(matmul_dot(A_small, B_small), A_small @ B_small)

# Read @ as "matrix product."  It is different from *, which multiplies
# corresponding entries.  For these arrays, both operations are defined,
# but they produce different answers:

P = np.array([[1.0, 2.0],
              [3.0, 4.0]])
Q = np.array([[0.0, 1.0],
              [1.0, 0.0]])

print("P * Q (elementwise):\n{}".format(P * Q))
print("\nP @ Q (matrix product):\n{}".format(P @ Q))

# TEST SPEED
#
# Both of your implementations are complete, so time them against numpy's
# built-in @ on the same matrices.  matmul_loop uses %time, which runs it
# ONCE.  Using %timeit there would run it repeatedly and make you wait
# much longer.

n = 200
A = np.random.random((n, n))
B = np.random.random((n, n))

print('time for matmul_loop (three loops, one run)')
%time C_slow = matmul_loop(A, B)

print('--------\n\ntime for matmul_dot (two loops)')
%timeit matmul_dot(A, B)

A @ B  # A throwaway call, so we do not time one-time library startup costs.

print('--------\n\ntime for A @ B (no loops)')
%timeit A @ B

from sklearn.datasets import load_digits
digits, labels = load_digits(return_X_y=True)


print(digits.shape)
print(labels.shape)

# Show the fist and second digits in the data set:
import matplotlib.pyplot as plt
plt.gray()
plt.matshow(digits[0, :].reshape(8, 8)) ## make into a square matrix for matshow
plt.show()
plt.matshow(digits[20, :].reshape(8, 8)) ## make into a square matrix for matshow
plt.show()

# "Broadcasting" can be used to apply the same operation across all rows of
# an array.  Here a length-4 array is subtracted from every row of a 3 x 4
# data matrix.

data = np.array([[3.0, 1.0, 7.0, 4.0],
                 [5.0, 9.0, 7.0, 0.0],
                 [1.0, 5.0, 7.0, 2.0]])

col_means = np.mean(data, axis=0)   # one mean per COLUMN, so length 4

print("data shape:      {}".format(data.shape))
print("col_means shape: {}".format(col_means.shape))
print("\ncolumn means: {}".format(col_means))

centered = data - col_means   # Broadcasting!  Shapes are 3x4 and 4

print("\ndata with each column's mean subtracted:\n{}".format(centered))

# The axis argument controls which direction an aggregation collapses.

# axis=0 collapses the rows, leaving one number per COLUMN.
print("sum down the columns: {}".format(np.sum(data, axis=0)))

# axis=1 collapses the columns, leaving one number per ROW.
print("sum across the rows:  {}".format(np.sum(data, axis=1)))

# Arithmetic is elementwise, so squaring first and then summing down the
# columns gives the sum of squared deviations for each column:
print("\nsum of squared deviations per column: {}".format(
    np.sum(centered ** 2, axis=0)))

# np.argmax returns the INDEX of the largest entry, not the entry itself.
# When there is a tie it returns the first index that holds the maximum.

example = np.array([12.0, 45.0, 3.0, 45.0])

print("largest value: {}".format(np.max(example)))
print("its index:     {}".format(np.argmax(example)))

def most_variable_feature(X):
    """ Return the index of the feature with the greatest variance.
         Arguments:
          X - An m x d array of m examples with d features each
         Returns:
          The column index of the feature with the largest variance
    """
    # YOUR CODE HERE
    raise NotImplementedError()

# Tests for most_variable_feature

# Column 1 spreads the most; column 2 never changes at all.
small = np.array([[2.0, 10.0, 5.0, 1.0],
                  [4.0, 14.0, 5.0, 3.0],
                  [6.0, 18.0, 5.0, 2.0]])
assert most_variable_feature(small) == 1

# Large values are not the same thing as large variance.  Column 0 holds the
# biggest numbers, but column 1 is the one that actually varies.
shifted = np.array([[100.0, 0.0],
                    [101.0, 8.0],
                    [99.0, -8.0]])
assert most_variable_feature(shifted) == 1

# A tie goes to the first such column, which is what np.argmax does.
tied = np.array([[0.0, 0.0],
                 [1.0, 1.0]])
assert most_variable_feature(tied) == 0

# And on the real data set:
assert most_variable_feature(digits) == 42

# Of course, numpy already has a name for "average squared distance from the
# mean".  Your function should agree with it:

print("your answer: {}".format(most_variable_feature(digits)))
print("np.var:      {}".format(np.argmax(np.var(digits, axis=0))))

# The per-pixel variance is itself an 8x8 image.  It shows where the pen
# actually goes.  The border pixels barely change, and a few never change at
# all, which means they can tell a classifier nothing.

pixel_variance = np.var(digits, axis=0)
print("\npixels that never vary: {}".format(np.flatnonzero(pixel_variance == 0)))

plt.matshow(pixel_variance.reshape(8, 8))
plt.title("variance of each pixel")
plt.show()

# Use numpy operators to pull out all columns
# with indices that are multiples of 3 or 7

X = np.random.random((3, 10))
print("X:\n{}".format(X))

all_cols = np.arange(0, X.shape[1])

print("\nAll column indices: {}".format(all_cols))

# Which are multiples of three?
cols_mult_3 = all_cols % 3 == 0

print("\nMultiples of three:")
print(cols_mult_3)

# Which are multiples of seven?
cols_mult_7 = all_cols % 7 == 0

print("\nMultiples of seven:")
print(cols_mult_7)

# We can use logical operators to combine boolean numpy arrays:
col_mask = cols_mult_3 | cols_mult_7

print("\nFinal mask:")
print(col_mask)

X_reduced = X[:, col_mask]

print("\nX_reduced:\n{}".format(X_reduced))

# Write code here to initialize the variables digits_subset and labels_subset.

# YOUR CODE HERE
raise NotImplementedError()

print(digits_subset.shape) # should be (361, 64)
print(labels_subset.shape) # should be (361,)

assert(digits_subset.shape == (361, 64))
assert(labels_subset.shape == (361,))
assert(labels_subset[0] == 1)
assert(labels_subset[-1] == 7)

# Look at some of the digits, assuming the variables have been initialized correctly.
plt.gray()
plt.matshow(digits_subset[0, :].reshape(8, 8)) # should be a "1"
plt.figure()
plt.matshow(digits_subset[180, :].reshape(8, 8))  # should be "7"
plt.show()

Dot Products, Matrix Products, and Numpy¶

Data Organization¶

How A Matrix Product Works¶

QUESTION¶

PART 1 - Processing Numpy Arrays Using Loops¶

Activity 1¶

Activity 2¶

PART 2 - Processing Numpy Arrays WITHOUT Loops¶

Activity 3¶

Activity 4¶

Matrix Products in Numpy¶

QUESTION¶

Part 3 - Working With Real Data¶

Activity 5 - Broadcasting and Feature Variance¶

Activity 6¶