#!/usr/bin/env python
# coding: utf-8

# ChEn-3170: Computational Methods in Chemical Engineering Spring 2020 UMass Lowell; Prof. V. F. de Almeida **12Feb20**
# 
# # Laboratory Work 05 (18Feb20)
# $
#   \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}}
#   \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}}
#   \newcommand{\Cmtrx}{\boldsymbol{\mathsf{C}}}
#   \newcommand{\Mmtrx}{\boldsymbol{\mathsf{M}}}
#   \newcommand{\Imtrx}{\boldsymbol{\mathsf{I}}}
#   \newcommand{\Pmtrx}{\boldsymbol{\mathsf{P}}}
#   \newcommand{\Lmtrx}{\boldsymbol{\mathsf{L}}}
#   \newcommand{\Umtrx}{\boldsymbol{\mathsf{U}}}
#   \newcommand{\xvec}{\boldsymbol{\mathsf{x}}}
#   \newcommand{\yvec}{\boldsymbol{\mathsf{y}}}
#   \newcommand{\zvec}{\boldsymbol{\mathsf{z}}}
#   \newcommand{\avec}{\boldsymbol{\mathsf{a}}}
#   \newcommand{\bvec}{\boldsymbol{\mathsf{b}}}
#   \newcommand{\cvec}{\boldsymbol{\mathsf{c}}}
#   \newcommand{\rvec}{\boldsymbol{\mathsf{r}}}
#   \newcommand{\norm}[1]{\bigl\lVert{#1}\bigr\rVert}
#   \DeclareMathOperator{\rank}{rank}
#   \DeclareMathOperator{\abs}{abs}
# $

# ### Name: `your name`

# ### Rubric for each assignment: 
# 
# |       Context          |  Points |
# | -----------------------     | ------- |
# | Precision of the answer     |   80%   |
# | Answer Markdown readability |   10%   |
# | Code readability            |   10%   |
#     

# ### <span style="color:red">Guidance:</span>
#  +  <span style="color:red"> 
#     Save your work frequently to a file locally to your computer.
#     </span>
#  +  <span style="color:red">
#     Before submitting, `Kernel` -> `Restart & Run All`, to verify your notebook runs correctly.
#     </span>
#  +  <span style="color:red">
#     Save your file again.
#     </span>

# ---
# ### Table of Contents
# * [Assignment 1 (35pts) ](#a1): $\Lmtrx$ forward solve.
#  - [1.1)](#a11) Code the forward solve in a `Python` function.
#  - [1.2)](#a12) Verify your forward solve code against `NumPy` (provide your own test).
#  - [1.3)](#a13) Test your forward solve code with a given matrix and right side vector.
# * [Assignment 2 (35pts) ](#a2): $\Umtrx$ forward solve.
#  - [2.1)](#a21) Code the backward solve in a `Python` function.
#  - [2.2)](#a22) Verify your backward solve code against `NumPy` (provide your own test).
#  - [2.3)](#a23) Test your backward solve code with a given matrix and right side vector.
# * [Assignment 3 (30pts)](#a3) $\Lmtrx\Umtrx$ factorization code and $\Lmtrx\Umtrx\xvec = \bvec$ solution.
#  - [3.1)](#a31) Code the LU factorization in a `Python` function and test against a given matrix.
#  - [3.2)](#a32) Verify your factorization code against `NumPy`.
#  - [3.3)](#a33) Test your $\Lmtrx\Umtrx\,\xvec = \bvec$ solution with a given matrix and right side vector.
#  ---

# ## <span style="color:blue">Assignment 1 (35 pts): For each item below respond in a separate notebook cell.</span><a id="a1"></a>

# <span style="color:blue"> 
# 1.1) Program a solution algorithm, as a `Python` function, for a system of equations with a lower triangular matrix of coefficients $\Lmtrx\,\xvec=\bvec$ (suggestion: use a test matrix and right side vector of *variable size* to test your code). The algorithm for solving $\Lmtrx\,\xvec=\bvec$ is as follows: 
# \begin{equation*}
#   x_i = \Bigl(b_i - \sum\limits_{j=1}^{i-1} L_{i,j}\,x_j \Bigr)\,L^{-1}_{i,i} \ \ \forall \ \ i=1,\ldots,m
# \end{equation*}
# for $i$ and $j$ with offset 1. Recall that `NumPy` and `Python` have offset 0 for their sequence data types. Provide an output example for $\Lmtrx$, $\bvec$, and $\xvec$. 
# </span><a id="a11"></a>

# In[1]:


'''1.1) Forward solve function'''


# In[ ]:


'''1.1) Test L x = b'''


# <span style="color:blue"> 1.2) Using the `numpy.linalg.solve()` function solve the same system, $\Lmtrx\,\yvec = \bvec$, print $\yvec$, and compute the norm of the
# difference of the solutions, i.e., $\norm{\xvec-\yvec}$.
#    </span><a id="a12"></a>

# In[4]:


'''1.2) Compare solutions'''


# <span style="color:blue">
# 1.3) Import the following image URL: 
#     
#  + https://raw.githubusercontent.com/dpploy/chen-3170/master/notebooks/images/cermet.png
# 
# as a matrix $\Amtrx$ (need internet connection) and extract from it a lower triangular matrix of the same dimensions, show the image, and solve for a right side vector $\bvec$ with all ones as elements (*i.e.* $b_i = 1\ \forall \ i=1,\ldots,m$), using **your** forward solve function. Compute the difference of your solution to the solution using `numpy.linalg.solve()` function and show the norm of the difference.
# </span></span><a id="a13"></a>

# In[6]:


'''1.3) A matrix and its L part'''


# In[7]:


'''1.3) Solve L y = b for b equal to a vector of ones'''


# ## <span style="color:blue">Assignment 2 (35 pts): For each item below respond in a separate notebook cell.</span><a id="a2"></a>

# <span style="color:blue">
# 2.1) Program a solution algorithm, as a `Python` function, for a system of equations with an upper triangular matrix of coefficients $\Umtrx\,\xvec=\bvec$ (suggestion: use a test matrix and right side vector of *variable size* to test your code). The algorithm for solving $\Umtrx\,\xvec=\bvec$  is as follows: 
# \begin{equation*}
# x_i = \Bigl(b_i - \sum\limits_{j=i+1}^{m} U_{i,j}\,x_j \Bigr)\,U^{-1}_{i,i} \ \ \forall \ \ i=m,\ldots,1
# \end{equation*}
# for $i$ and $j$ with offset 1. Recall that `NumPy` and `Python` have offset 0 for their sequence data types. Provide an output example for $\Umtrx$, $\bvec$, and $\xvec$.
# </span><a id="a21"></a>

# In[8]:


'''2.1) Backward solve function'''
 


# In[ ]:


'''2.1) Example U x = b'''


# <span style="color:blue">
# 2.2) Using the `numpy.linalg.solve()` function solve the same system, $\Umtrx\,\yvec = \bvec$, print $\yvec$, and compute the norm of the difference of the solutions, i.e., $\norm{\xvec-\yvec}$.
# </span><a id="a22"></a>

# In[11]:


'''2.2) Compare solutions'''


# <span style="color:blue">
# 2.3) Extract from $\Amtrx$ an upper triangular matrix of the same dimensions, show the image, and solve for a right side vector $\bvec$ with all ones as elements (<i>i.e.</i> $b_i = 1\ \forall \ i=1,\ldots,m$), using <b>your</b> backward solve function. Compute the difference of your solution to the solution using the `numpy.linalg.solve()` function and show the norm of the difference.
# </span></span><a id="a23"></a>

# In[12]:


'''2.3) A matrix and its U part'''


# In[13]:


'''2.3) Solve U y = b for b equal to a vector of ones'''


# ## <span style="color:blue">Assignment 3 (30 pts): For each item below respond in a separate notebook cell.</span><a id="a3"></a>

# <span style="color:blue">
# 3.1) Program (in a function) a $\Lmtrx\,\Umtrx$ factorization algorithm <b>(without using pivoting)</b> for the matrix $\overset{(120x120)}{\Amtrx}$ and compute the $\Lmtrx\,\Umtrx$ factors. Compute the $\max\bigl(\abs(\Amtrx-\Lmtrx\,\Umtrx)\bigr)$. The factorization is obtained by elimination steps $k = 1,\ldots,m-1$ so that 
#     
# \begin{equation*}
#     A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\, m_{i,k} \quad\quad \forall \quad\quad i=k+1,\ldots,m \quad\quad \text{and} \quad\quad j=k+1,\ldots,m
# \end{equation*}
# 
# where the multipliers $m_{i,k}$ are given by $m_{i,k} = \frac{A^{(k)}_{i,k}}{A^{(k)}_{k,k}}$. When $k = m-1$, $A^{(m)}_{i,j}$, is upper triangular, that is, $U_{i,j} = A^{(m)}_{i,j}$ . The lower triangular matrix is obtained using the multipliers $m_{i,k}$, that is $L_{i,j} = m_{i,j} \ \forall \ i>j$,  $L_{i,i}=1$, and $L_{i,j}=0 \ \forall \ i<j$.
# </span><a id="a31"></a>

# In[14]:


'''3.1) LU factorization function'''


# In[17]:


'''3.1) Compute LU factors of A 120x120'''


# <span style="color:blue">
# 3.2) Print the smallest magnitude diagonal entry of the $\Umtrx$ factor, compute the $\rank(\Amtrx)$ and the $\det(\Amtrx)$ using the $\Lmtrx\,\Umtrx$ factors.
# </span><a id="a32"></a>

# In[19]:


'''3.2 Compute: min(diag(U)), rank(A), det(A) using the LU factors'''


# <span style="color:blue">
# 3.3) Using your codes (Assignment 1 and 2) solve the linear system of algebraic equations $\Amtrx\,\xvec=\bvec$ for $\bvec$ with all ones as elements (*i.e.* $b_i = 1\ \forall \ i=1,\ldots,m$). Solve the same system with `numpy.linalg.solve` and compute the norm of the solutions difference.
# </span><a id="a33"></a>

# In[20]:


'''3.3 Solve A x = b for b equal to a vector of ones '''