# coding: utf-8

# # CSC338. Homework 4
# 
# Due Date: Wednesday Feburary 5, 9pm
# 
# Please see the guidelines at https://www.cs.toronto.edu/~lczhang/338/homework.html
# 
# ### What to Hand In
# 
# Please hand in 2 files:
# 
# - Python File containing all your code, named `hw4.py`.
# - PDF file named `hw4_written.pdf` containing your solutions to the written parts of the assignment. Your solution can be hand-written, but must be legible. Graders may deduct marks for illegible or poorly presented solutions.
# 
# If you are using Jupyter Notebook to complete the work, your notebook can be exported as a .py file (File -> Download As -> Python). Your code will be auto-graded using Python 3.6, so please make sure that your code runs. There will be a 20% penalty if you need a remark due to small issues that renders your code untestable.
# 
# **Make sure to remove or comment out all matplotlib or other expensive code before submitting your homework!**
# 
# Submit the assignment on **MarkUs** by 9pm on the due date. See the syllabus for the course policy regarding late assignments. All assignments must be done individually.

# In[ ]:


import math
import numpy as np


# ## Question 1
# 
# For this question, we will again start from code from tutorial 3.

# In[ ]:


# Code from tutorial 3

def backward_substitution(A, b):
    """Return a vector x with np.matmul(A, x) == b, where 
        * A is an nxn numpy matrix that is upper-triangular and non-singular
        * b is an nx1 numpy vector
    """
    n = A.shape[0]
    x = np.zeros_like(b, dtype=np.float)
    for i in range(n-1, -1, -1):
        s = 0
        for j in range(n-1, i, -1):
            s += A[i,j] * x[j]
        x[i] = (b[i] - s) / A[i,i]
    return x

def eliminate(A, b, k):
    """Eliminate the k-th row of A, in the system np.matmul(A, x) == b,
    so that A[i, k] = 0 for i < k. The elimination is done in place."""
    n = A.shape[0]
    for i in range(k + 1, n):
        m = A[i, k] / A[k, k]
        for j in range(k, n):
            A[i, j] = A[i, j] - m * A[k, j]
        b[i] = b[i] - m * b[k]

def gauss_elimination(A, b):
    """Return a vector x with np.matmul(A, x) == b using
    the Gauss Elimination algorithm, without partial pivoting."""
    for k in range(A.shape[0] - 1):
        eliminate(A, b, k)
    x = backward_substitution(A, b)
    return x


# ### Part (a) -- 1 pt
# 
# Solve the system $Ax = b$ for the values of $A$ and $b$ below using
# the function `gauss_elimination` from last time.
# Save the solution you obtain in the variable `soln_nopivot`.

# In[ ]:


e = pow(2, -100)
A = np.array([[e, 1],
              [1, 1]])
b = np.array([1 + e, 2])

soln_nopivot = None


# ### Part (b) -- 2 pt
# 
# Write a helper function `partial_pivot` that performs partial pivoting
# on `A` at column `k`, so that the function 
# `gauss_elimination_partial_pivot` performs Gauss Elimination with Partial Pivoting.

# In[ ]:


def partial_pivot(A, b, k):
    """Perform partial pivoting for column k. That is, swap row k
    with row j > k so that the new element at A[k,k] is the largest
    amongst all other values in column k below the diagonal.
    
    This function should modify A and b in place.
    """
    # TODO

def gauss_elimination_partial_pivot(A, b):
    """Return a vector x with np.matmul(A, x) == b using
    the Gauss Elimination algorithm, with partial pivoting."""
    for k in range(A.shape[0] - 1):
        partial_pivot(A, b, k)
        eliminate(A, b, k)
    x = backward_substitution(A, b)
    return x


# ### Part (c) -- 1 pt
# 
# Solve the system $Ax = b$ for the values of $A$ and $b$ below using `gauss_elimination_partial_pivot`.  Save the solution you obtain in the variable `soln_pivot`.

# In[ ]:


e = pow(2, -100)
A = np.array([[e, 1],
              [1, 1]])
b = np.array([1 + e, 2])

soln_pivot = None


# ### Part (d) -- 2 pt
# 
# Do your answers in parts (a) and (d) match? If not, which is the correct
# answer? Include your explanation in the PDF write-up.
# 
# ## Question 2
# 
# ### Part (a) -- 3 pt
# 
# Consider the following matrices `M1`, `M2`, and `M3`. Compute each of their $L_1$, $L_2$, and $L_\infty$ norms. Save the results in the variables below.
# 
# For the $L_2$ norm you may find the function `np.linalg.norm` helpful. You can compute the $L_1$ and $L_\infty$ norms either by hand or write a function.

# In[ ]:


M1 = np.array([[3., 0.],
               [-4., 2.]])

M2 = np.array([[2., -2., 0.3],
               [0.5, 1., 0.9],
               [-4., -2., 5]])

M3 = np.array([[0.2, -0.2],
               [1.0, 0.2]])

# fill in these answers
M1_l_1 = 0
M1_l_2 = 0
M1_l_infty = 0
M2_l_1 = 0
M2_l_2 = 0
M2_l_infty = 0
M3_l_1 = 0
M3_l_2 = 0
M3_l_infty = 0


# ### Part (b) -- 4 pt
# 
# Show that an induced matrix norm $|| \cdot ||$ has the property
# $$||A + B|| \leq ||A|| + ||B||$$
# 
# ### Part (c) -- 4 pt
# 
# Is it true that for a vector, $||v||_\infty \le ||v||_1$? What about for a matrix: is it true that for a matrix, $||M||_\infty \le ||M||_1$? Include your solution and justification in your pdf writeup.
# 
# ## Question 3
# 
# ### Part (a) [3 pt]
# 
# Write a function `matrix_condition_number` that computes the condition
# number of a $2 \times 2$ matrix. Use the $$L_1$$ matrix norm.

# In[ ]:


def matrix_condition_number(M):
    """
    Returns the condition number of the 2x2 matrix M.
    Use the $L_1$ matrix norm.

    Precondition: M.shape == [2, 2] 
                  M is non-singular

    >>> matrix_condition_number(np.array([[1., 0.], [0., 1.]]))
    1
    """


# ### Part (b) [2 pt]
# 
# Classify each of the following matrices `A1`, `A2`, `A3` and `A4`,
# as well-conditioned or ill-conditioned.
# 
# You may do this question either by hand, or by using the function above.
# 
# Save the classifications in a Python array called `conditioning`.
# Each item of the array
# should be either the string "well" or the string "ill".

# In[ ]:


A1 = np.array([[3, 0], 
               [0, pow(3, -30)]])
A2 = np.array([[pow(3, 20), 0],
               [0,          pow(5, 50)]])
A3 = np.array([[pow(4, -40), 0],
               [0,          pow(2, -80)]])
A4 = np.array([[pow(5, -17), pow(5, -16)],
               [pow(5, -18), pow(5, -17)]])

# whether A1, A2, A3, A4 is "well" or "ill" conditioned
conditioning = [None, None, None, None]


# ### Part (c) [2 pt]
# 
# It should be immediate obvious that the matrix
# 
# $$\begin{bmatrix}100 & 2 \\ 201 & 4\end{bmatrix}$$
# 
# is ill-conditioned. Explain why. Include your answer in your pdf write-up.
# 
# 
# ### Part (d) [2 pt]
# 
# Suppose that $A$ and $B$ are two $n \times n$ matrices, and both are well-conditioned. Is $A(B^{-1})$ also well-conditioned? Why or why not? Include your answer and justifications in your pdf write-up. Be specific.
# 
# ### Part (e) [4 pt]
# 
# Describe an efficient algorithm to compute $d^T B^T A^{-1} B d$
# 
# Where:
# 
# * $A$ is an invertible $n \times n$ matrix, 
# * $B$ is an $n \times n$ matrix, and 
# * $d$ is an $n \times 1$ vectors
# 
# Be clear and specific. Include your strategy in your pdf write-up.