Python-Workshop-3-NumPy.md

@@ -70,7 +70,8 @@ The first few examples will be short and focus on arrays and statistics. These
 examples are inspired by [*Coffee Break NumPy* by Christian
 Mayer](https://www.amazon.com/Coffee-Break-NumPy-Science-Mastery/dp/1076932614).
 
-**Creating Arrays**
+
+## Creating Arrays
 
 ```python
     array_size = 8
@@ -108,6 +109,45 @@ Mayer](https://www.amazon.com/Coffee-Break-NumPy-Science-Mastery/dp/1076932614).
     print()
 ```
 
+
+## Broadcasting
+
+```python
+    prices = [1.00, 2.95, 8.40, 3.50, 3.30, 16.91]
+    prices = [0.9 * price for price in prices]
+
+    print(prices)
+```
+
+```python
+    prices = np.array([1.00, 2.95, 8.40, 3.50, 3.30, 16.91])
+    prices *= 0.9
+
+    print(prices)
+
+    print()
+    print("*" * 80)
+    print()
+```
+
+```python
+    def op_wrapper_py():
+        prices = range(1, 1000000, 1)
+        prices = [0.9 * price for price in prices]
+
+    py_list = timeit.timeit(op_wrapper_py, number=100)
+
+    def op_wrapper_np():
+        prices = np.arange(0, 1000000, 1, dtype=np.float64)
+        prices[:] *= 0.9
+
+    np_array = timeit.timeit(op_wrapper_np, number=100)
+
+    print(f"Python Time: {py_list:.4f}")
+    print(f"NumPy Time : {np_array:.4f}")
+```
+
+
 **I/O**
 
 **Indexing**
@@ -116,11 +156,126 @@ Mayer](https://www.amazon.com/Coffee-Break-NumPy-Science-Mastery/dp/1076932614).
 
 **Boolean (Mask) Index Arrays**
 
-**Broadcasting**
 
-**Linear Algebra**
+# Fun Examples!
+
+The examples in this section are fun!
+
+
+## Linear Algebra
+
+In CS 417/517 Computational Methods... I require students to [implement a
+Matrix Solver... from
+scratch](https://www.cs.odu.edu/~tkennedy/cs417/s21/Assts/matrixSolverExercise/).
+NumPy provides implementations of matrix operations (e.g., multiplication).
+
+Let us implement a quick NumPy based Matrix solver, for this [Discrete Case
+Least Squares Approximation
+Problem](https://www.cs.odu.edu/~tkennedy/cs417/s21/Public/approximationExample2/).
+
+```python
+import numpy as np
+
+def print_matrices(matrix_XTX, matrix_XTY):
+    """
+    Print the XTX and XTY matrices
+    """
+
+    print("{:*^40}".format("XTX"))
+    print(matrix_XTX)
+
+    print()
+    print("{:*^40}".format("XTY"))
+    print(matrix_XTY)
+
+
+def _backsolve(matrix_XTX, matrix_XTY):
+
+    num_rows, _ = matrix_XTX.shape
+
+    for i in reversed(range(1, num_rows)):
+        for j in reversed(range(0, i)):
+            s = matrix_XTX[j, i]
+
+            matrix_XTX[j, i] -= (s * matrix_XTX[i, i])
+            matrix_XTY[j] -= (s * matrix_XTY[i])
+
+
+def solve_matrix(matrix_XTX, matrix_XTY):
+    """
+    Solve a matrix and return the resulting solution vector
+    """
+
+    # Get the dimensions (shape) of the XTX matrix
+    num_rows, num_columns = matrix_XTX.shape
+
+    for i in range(0, num_rows):
+        # Find column with largest entry
+        largest_idx = i
+        current_col = i
+        for j in range(i + 1, num_rows):
+
+            if matrix_XTX[largest_idx, i] < matrix_XTX[j, current_col]:
+                largest_idx = j
+
+        # Swap
+        if largest_idx != current_col:
+            matrix_XTX[[i, largest_idx], :] = matrix_XTX[[largest_idx, i], :]
+            matrix_XTY[[i, largest_idx]] = matrix_XTY[[largest_idx, i]]
+
+        # Scale
+        scaling_factor = matrix_XTX[i, i]
+        matrix_XTX[i, :] /= scaling_factor
+        matrix_XTY[i] /= scaling_factor
+
+        # Eliminate
+        for row_i in range(i + 1, num_rows):
+            s = matrix_XTX[row_i][i]
+
+            matrix_XTX[row_i] = matrix_XTX[row_i] - s * matrix_XTX[i]
+            matrix_XTY[row_i] = matrix_XTY[row_i] - s * matrix_XTY[i]
+
+        #  print("{:-^80}".format(f"Iteration #{i:}"))
+        #  print_matrices(matrix_XTX, matrix_XTY)
+
+    _backsolve(matrix_XTX, matrix_XTY)
+
+    return matrix_XTY
+
+
+def main():
+
+    # Set up input data points, X, Y, and XT
+    points = [(0., 0.), (1., 1.), (2., 4.)]
+
+    matrix_X = np.array([[1., 0., 0.],
+                         [1., 1., 1.],
+                         [1., 2., 4.]])
+
+    matrix_Y = np.array([0,
+                         1,
+                         4])
+
+    matrix_XT = matrix_X.transpose()
+
+    # Compute XTX and XTY
+    matrix_XTX = np.matmul(matrix_XT, matrix_X)
+    matrix_XTY = np.matmul(matrix_XT, matrix_Y)
+
+    print_matrices(matrix_XTX, matrix_XTY)
+
+    print()
+    print("{:-^40}".format("Solution"))
+    solution = solve_matrix(matrix_XTX, matrix_XTY)
+    print(solution)
+
+
+if __name__ == "__main__":
+    main()
+```
+
 
-**Statistics**
+## Statistics