ChEn-3170: Computational Methods in Chemical Engineering Spring 2020 UMass Lowell; Prof. V. F. de Almeida 28Jan20
06. Computational Linear Algebra Fundamentals¶
$ \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Cmtrx}{\boldsymbol{\mathsf{C}}} \newcommand{\Dmtrx}{\boldsymbol{\mathsf{D}}} \newcommand{\Mmtrx}{\boldsymbol{\mathsf{M}}} \newcommand{\Imtrx}{\boldsymbol{\mathsf{I}}} \newcommand{\Pmtrx}{\boldsymbol{\mathsf{P}}} \newcommand{\Qmtrx}{\boldsymbol{\mathsf{Q}}} \newcommand{\Lmtrx}{\boldsymbol{\mathsf{L}}} \newcommand{\Umtrx}{\boldsymbol{\mathsf{U}}} \newcommand{\xvec}{\boldsymbol{\mathsf{x}}} \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} $
Table of Contents¶
- Matrix solve
- $\Lmtrx$ forward solve
- $\Umtrx$ backward solve
- $\Amtrx = \Pmtrx\,\Lmtrx\,\Umtrx$ factorization
- $\Lmtrx$ forward solve
- $\Umtrx$ backward solve
- $\Amtrx = \Lmtrx\,\Umtrx$ factorization
- $\Pmtrx\,\Amtrx = \Lmtrx\,\Umtrx$ factorization
- $\Pmtrx\,\Amtrx\,\Qmtrx = \Lmtrx\,\Umtrx$ factorization
Objectives¶
- Introduce the elements of computational linear algebra needed in this course to analyse and solve system of linear algebraic equations.
- Implement a direct method of solution of linear algebraic equations (also known as matrix factorization).
Theory¶
The course notes (OneNote ChEn-3170-linalg cover basic elements of linear system of algebraic equations as applied to computational stoichiometry. Particular attention is given to conditions for the existance and uniqueness of solutions of general algebraic systems.
Basic theoretical aspects of solving for $\overset{(n)}{\xvec}$ in the matrix equation $\overset{(m\times n)}{\Amtrx}\,\overset{(n)}{\xvec} = \overset{(m)}{\bvec}$ are covered. $\overset{(m\times n)}{\Amtrx}$ is a matrix, $\overset{(m)}{\bvec}$ and $\overset{(n)}{\xvec}$ are vectors where $m$ indicates the number of rows (or equations) and $n$ number of columns (or unknowns).
Matrix-vector and matrix-matrix product operations¶
The following operations between vectors and matrices are obtained directly from the buil-in functions in the numpy package.
'''Import the NumPy package as usual'''
import numpy as np
Inner product of two vectors: $\avec \cdot \bvec$.
'''Vector inner product or dot product of vectors'''
a_vec = np.array( np.random.random(3) )
b_vec = np.array( np.random.random(3) )
np.set_printoptions( precision=3, threshold=20, edgeitems=12, linewidth=100 )
a_vec_dot_b_vec = np.dot( a_vec, b_vec ) # clear linear algebra operation
print('a.b =', a_vec_dot_b_vec)
a_vec_x_b_vec = a_vec @ b_vec # consistent linear algebra multiplication
print('a@b =', a_vec_x_b_vec )
a.b = 0.6381292212080929 a@b = 0.6381292212080929
print( 'a.b = %10.4e'%a_vec_dot_b_vec ) # formatting with scientific notation
a.b = 6.3813e-01
Matrix vector product: $\Amtrx\,\bvec$.
'''Matrix-vector product'''
a_mtrx = np.array( [ [ 2., 1., 1.], # per course notes (NB 03/04)
[ 4., -6., 0.],
[-2., 7., 2.]
] )
b_vec = np.array( [5., -2., 9.] ) # per course notes
a_mtrx_x_b_vec = a_mtrx @ b_vec # linear algebra matrix-vector product
print('A b =', a_mtrx_x_b_vec)
A b = [17. 32. -6.]
Matrix-vector product: $\Imtrx\,\bvec = \bvec$. Note: $\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \, \begin{pmatrix} b_1\\b_2\\b_3 \end{pmatrix} = \begin{pmatrix} b_1\\b_2\\b_3 \end{pmatrix} $.
'''Identity-matrix vector product'''
i_mtrx = np.eye(3)
i_mtrx_x_b_vec = i_mtrx @ b_vec # linear algebra matrix-vector product
print('b =', b_vec)
print('I x b =', i_mtrx_x_b_vec)
b = [ 5. -2. 9.] I x b = [ 5. -2. 9.]
Matrix-matrix product: $\Imtrx\,\Amtrx = \Amtrx$. Note: $\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \, \begin{pmatrix} A_{1,1} & A_{1,2} & A_{1,3} \\ A_{2,1} & A_{2,2} & A_{2,3} \\ A_{3,1} & A_{3,2} & A_{3,3} \end{pmatrix} = \begin{pmatrix} A_{1,1} & A_{1,2} & A_{1,3} \\ A_{2,1} & A_{2,2} & A_{2,3} \\ A_{3,1} & A_{3,2} & A_{3,3} \end{pmatrix} $.
'''Matrix-matrix product IA = A'''
i_mtrx_x_a_mtrx = i_mtrx @ a_mtrx # linear algebra matrix-matrix product
print('I x A =\n', i_mtrx_x_a_mtrx)
print('A =\n', a_mtrx)
I x A = [[ 2. 1. 1.] [ 4. -6. 0.] [-2. 7. 2.]] A = [[ 2. 1. 1.] [ 4. -6. 0.] [-2. 7. 2.]]
Matrix-matrix product: $\Amtrx\,\Bmtrx = \Cmtrx$. Note: $\begin{pmatrix} A_{1,1} & A_{1,2} & A_{1,3} \\ A_{2,1} & A_{2,2} & A_{2,3} \\ A_{3,1} & A_{3,2} & A_{3,3} \end{pmatrix} \, \begin{pmatrix} B_{1,1} & B_{1,2} & B_{1,3} \\ B_{2,1} & B_{2,2} & B_{2,3} \\ B_{3,1} & B_{3,2} & B_{3,3} \end{pmatrix} = \begin{pmatrix} C_{1,1} & C_{1,2} & C_{1,3} \\ C_{2,1} & C_{2,2} & C_{2,3} \\ C_{3,1} & C_{3,2} & C_{3,3} \end{pmatrix} $ where each $C_{i,j}$ is a vector product of the $i$th row of $\Amtrx$ and the $j$th column of $\Bmtrx$, i.e. $C_{i,j} = \sum\limits_{k=1}^3 A_{i,k}\, B_{k,j}$.
'''Matrix-matrix product AB = C'''
b_mtrx = np.array( [ [ 5. , 5. , 5.],
[-2. , -2. , -2.],
[ 9. , 9. , 9.] ]
)
c_mtrx = a_mtrx @ b_mtrx # linear algebra matrix-matrix product
print('A =\n', a_mtrx)
print('B =\n', b_mtrx)
print('C =\n', c_mtrx)
A = [[ 2. 1. 1.] [ 4. -6. 0.] [-2. 7. 2.]] B = [[ 5. 5. 5.] [-2. -2. -2.] [ 9. 9. 9.]] C = [[17. 17. 17.] [32. 32. 32.] [-6. -6. -6.]]
'''Matrix-matrix product BA = D'''
b_mtrx = np.array( [[5. , 5. , 5.],
[-2., -2., -2.],
[9. , 9. , 9.]]
)
d_mtrx = b_mtrx @ a_mtrx # linear algebra matrix-matrix product
print('A =\n', a_mtrx)
print('B =\n', b_mtrx)
print('D =\n', d_mtrx)
A = [[ 2. 1. 1.] [ 4. -6. 0.] [-2. 7. 2.]] B = [[ 5. 5. 5.] [-2. -2. -2.] [ 9. 9. 9.]] D = [[20. 10. 15.] [-8. -4. -6.] [36. 18. 27.]]
NumPy and SciPy Linear Algebra¶
NumPy has extensive support for linear algebra arrays. We collect here the relevant operations for this course. However additional resources are instead added to SciPy for general scientific computing including linear algebra.
Linear algebra operations are obtained from the linalg sub-package of the numpy package, and the linalg sub-package of scipy.
'''Import the NumPy linear algebra sub-package as usual'''
#import numpy.linalg as linalg
#from numpy import linalg # often used alternative
'''or leave it commented since the usage of np.linalg is self-documenting'''
'or leave it commented since the usage of np.linalg is self-documenting'
The 2-norm or norm (magnitude) of a vector $\bvec$ is indicated as $\norm{\bvec}$ and computed as follows:
'''Vector norm (or magnitude)'''
norm_b_vec = np.linalg.norm( b_vec ) # default norm is the 2-norm
print('||b|| =', norm_b_vec) # same as magnitude
||b|| = 10.488088481701515
Matrix solve ¶
Solve for $\xvec$ in the matrix equation $\Amtrx\,\xvec = \bvec$, where $\Amtrx = \begin{pmatrix} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2 \end{pmatrix} $ and $\bvec = \begin{pmatrix} 5\\ -2\\ 9 \end{pmatrix}$.
'''Matrix solver (this is short for solution of a linear algebraic system of equations)'''
x_vec = np.linalg.solve( a_mtrx, b_vec ) # solve linear system for A, b; per course notes 04
print('solution x =', x_vec)
solution x = [1. 1. 2.]
The residual vector defined as $\rvec = \bvec - \Amtrx\,\xvec$ is of importance. So is its norm $\norm{\rvec}$.
'''Verify the accuracy of the solution'''
res_vec = b_vec - a_mtrx @ x_vec
print('b - A x =',res_vec)
print('||b - A x|| =',np.linalg.norm( res_vec ))
b - A x = [0. 0. 0.] ||b - A x|| = 0.0
The rank of a matrix of coefficients, $\rank(\Amtrx)$, of a linear algebraic system of equations determines weather the solution is unique or singular.
'''Matrix rank'''
k = np.linalg.matrix_rank( a_mtrx ) # rank; per course notes 14
print('rank(A) =',k)
print('shape(A) =',a_mtrx.shape)
if k == a_mtrx.shape[0] and k == a_mtrx.shape[1]: # flow control
print('A is non-singular; solution is unique ')
else:
print('A is singular')
rank(A) = 3 shape(A) = (3, 3) A is non-singular; solution is unique
Why is this matrix $\Bmtrx$ singular?
b_mtrx = np.array( [ [ 2., 1., 3.], # singular
[ 4., -6., -2.],
[-2., 7., 5.]])
k = np.linalg.matrix_rank( b_mtrx ) # rank
print('rank(B) =',k)
print('shape(B) =',b_mtrx.shape)
if k == b_mtrx.shape[0] and k == b_mtrx.shape[1]: # flow control
print('B is non-singular; solution is unique ')
else:
print('B is singular')
rank(B) = 2 shape(B) = (3, 3) B is singular
'''Matrix determinant'''
det_a_mtrx = np.linalg.det( a_mtrx ) # determinant;
print('det(A) =', det_a_mtrx)
det(A) = -15.999999999999998
'''Matrix determinant'''
det_b_mtrx = np.linalg.det( b_mtrx ) # determinant of a singular matrix
print('det(B) =', det_b_mtrx)
det(B) = 0.0
The inverse matrix is denoted as $\Amtrx^{-1}$ and is computed as the matrix that multiplies $\bvec$ and produces the solution $\xvec$, that is, $\xvec = \Amtrx^{-1}\,\bvec$.
'''Matrix inverse'''
a_mtrx_inv = np.linalg.inv( a_mtrx ) # matrix inverse; per course notes 17
print('A^-1 =\n', a_mtrx_inv)
A^-1 = [[ 0.75 -0.312 -0.375] [ 0.5 -0.375 -0.25 ] [-1. 1. 1. ]]
Recall $\Amtrx^{-1}\,\Amtrx = \Imtrx$ where $\Imtrx$ is the identity matrix.
'''Identity matrix'''
i_mtrx = a_mtrx_inv @ a_mtrx # identity matrix; per course notes 17
print('A^-1 A =\n',i_mtrx)
A^-1 A = [[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]]
Using the inverse, the same solution will be found: $\xvec = \Amtrx^{-1}\,\bvec$.
'''Solution using the inverse'''
x_vec_again = a_mtrx_inv @ b_vec # matrix-vector multiply; per course notes 17
print('solution x =', x_vec_again)
solution x = [1. 1. 2.]
This is the element-by-element reciprocal of the matrix $(\Amtrx)^{-1}$, which is very different than the inverse.
'''Inverse power of a matrix'''
#a_mtrx_to_negative_1 = a_mtrx**(-1) # this will cause an error (division by zero)
'Inverse power of a matrix'
Let's look at the determinant of a larger matrix, say $\Mmtrx$.
'''Generate a larger matrix from an image'''
%matplotlib inline
from matplotlib import pyplot as plt # import the pyplot function of the matplotlib package
plt.rcParams['figure.figsize'] = [20, 4] # extend the figure size on screen output
mtrx = plt.imread('https://raw.githubusercontent.com/dpploy/chen-3170/master/notebooks/images/cermet.png',format='png')
m_mtrx = mtrx.astype('float64')
plt.figure(1) # create a figure placeholder
plt.imshow( m_mtrx,cmap='gray')
plt.title('Cermet',fontsize=14)
plt.xlabel('x pixels',fontsize=12)
plt.ylabel('y pixels',fontsize=12)
plt.show()
print('M shape =', m_mtrx.shape)
M shape = (500, 500)
'''Larger matrix determinant'''
det_m_mtrx = np.linalg.det( m_mtrx ) # determinant
print('max(M) =',m_mtrx.max())
print('min(M) =',m_mtrx.min())
print('det(M) = %10.3e (not an insightful number)'%det_m_mtrx) # formatting numeric output
print('rank(M) = ',np.linalg.matrix_rank( m_mtrx, tol=1e-5 ) )
max(M) = 0.8901960849761963 min(M) = 0.062745101749897 det(M) = 0.000e+00 (not an insightful number) rank(M) = 500
Let's solve for this matrix with $\cvec$ as the right side vector, that is, $\Mmtrx\,\xvec = \cvec$.
'''Solve M x = c and plot x'''
c_vec = np.random.random(mtrx.shape[0]) # any c will do it
sol = np.linalg.solve( m_mtrx, c_vec ) # solve linear system for A, b
plt.figure(2)
plt.plot(range(c_vec.size),sol,'k')
plt.title('M x = c',fontsize=20)
plt.xlabel('n',fontsize=18)
plt.ylabel('$c_j$',fontsize=18)
print('')
res_vec = c_vec - m_mtrx @ sol
#print('c - M x =',res_vec)
print('||c - M x|| =%12.4e'%np.linalg.norm( res_vec ))
||c - M x|| = 8.3733e-12
$\Pmtrx\,\Lmtrx\,\Umtrx$ factorization¶
The factors: $\Pmtrx$, $\Lmtrx$, and $\Umtrx$ where $\Pmtrx\,\Lmtrx\,\Umtrx = \Amtrx$ can be obtained from the SciPy linear algebra package. $\Pmtrx$ is a permutation matrix if the underlying Gaussian elimination is used to construct the $\Lmtrx$ and $\Umtrx$ factors.
'''Import only the linear algebra package'''
import scipy.linalg
import numpy as np
'''P L U factors of A'''
a_mtrx = np.array( [[1, 2, 3],
[4, 5, 6],
[7, 8, 10]] )
(p_mtrx, l_mtrx, u_mtrx) = scipy.linalg.lu( a_mtrx )
print('P =\n',p_mtrx)
print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('Checking...')
print('PLU - A =\n', p_mtrx @ l_mtrx @ u_mtrx - a_mtrx)
P = [[0. 1. 0.] [0. 0. 1.] [1. 0. 0.]] L = [[1. 0. 0. ] [0.143 1. 0. ] [0.571 0.5 1. ]] U = [[ 7. 8. 10. ] [ 0. 0.857 1.571] [ 0. 0. -0.5 ]] Checking... PLU - A = [[0. 0. 0.] [0. 0. 0.] [0. 0. 0.]]
'''P^-1 = P^T (i.e. the transpose of a permutation matrix is its inverse)'''
pinv_mtrx = np.linalg.inv(p_mtrx)
print('P^-1 =\n', pinv_mtrx)
print('Checking...')
print('P^-1 - P^T =\n', pinv_mtrx - p_mtrx.transpose())
P^-1 = [[0. 0. 1.] [1. 0. 0.] [0. 1. 0.]] Checking... P^-1 - P^T = [[0. 0. 0.] [0. 0. 0.] [0. 0. 0.]]
'''PLU x = b; that is: Forward: L y = P^-1 b, Backward: U x = y '''
b_vec = np.array([1.,2.,3.])
y_vec = scipy.linalg.solve(l_mtrx, p_mtrx.transpose() @ b_vec) # L y = P^T b
x_vec = scipy.linalg.solve(u_mtrx, y_vec) # U x = y
print('x =', x_vec)
x_vec_gold = scipy.linalg.solve( a_mtrx, b_vec ) # solution using A x = b
print('||x - x_gold|| =',scipy.linalg.norm(x_vec-x_vec_gold))
x = [-3.333e-01 6.667e-01 3.172e-17] ||x - x_gold|| = 0.0
'''Deterninant of U or L: product of the diagonal'''
det_u = np.linalg.det(u_mtrx)
print('det(U) = %8.3e'%det_u)
diag_vec = np.diagonal(u_mtrx)
prod = np.prod(diag_vec)
print('diag(U) product = %8.3e'%prod )
det(U) = -3.000e+00 diag(U) product = -3.000e+00
'''Determinant of P (always +1 or -1)'''
det_p = np.linalg.det(p_mtrx)
print('det(P) = %8.3e'%det_p)
det(P) = 1.000e+00
'''Determinant of A = det(PLU)'''
det_l = np.prod( np.diagonal(l_mtrx) )
det_plu = det_p * det_l * det_u # last term is det of L
print('det(PLU) = %8.3e'%det_plu)
print('det(A) = %8.3e'%np.linalg.det(a_mtrx))
det(PLU) = -3.000e+00 det(A) = -3.000e+00
$\Lmtrx$ forward solve¶
A lower triangular matrix like any matrix can be used in a matrix solve.
'''L forward solve'''
l_mtrx = np.array( [[1., 0., 0.], # per course notes
[2., 3., 0.],
[4., 5., 6.]] )
b_vec = np.array( [1.,2.,3.] )
x_vec = np.linalg.solve( l_mtrx, b_vec )
np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)
x = [ 1.000e+00 -1.110e-17 -1.667e-01]
$\Umtrx$ backward solve¶
An upper triangular matrix like any matrix can be used in a matrix solve.
'''U backward solve'''
u_mtrx = np.array( [[1., 2., 3.], # per course notes
[0, 4., 5.],
[0., 0., 6.]] )
b_vec = np.array( [1.,2.,3.] )
x_vec = np.linalg.solve( u_mtrx, b_vec )
np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)
x = [-0.25 -0.125 0.5 ]
ChEn-3170 Linear Algebra¶
In this course various algorithms need to be programmed. These should be compared to SciPy and/or NumPy.
$\Lmtrx$ forward solve¶
A lower triangular matrix allows for a forward solve. The algorithm for $\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} \quad\ \forall \quad\ i=1,\ldots,m \end{equation*}for $i$ and $j$ with offset 1. See Python implementation below.
'''L forward solve'''
l_mtrx = np.array( [[1., 0., 0.], # per course notes
[2., 3., 0.],
[4., 5., 6.]] )
b_vec = np.array( [1.,2.,3.] )
from chen_3170.help import forward_solve # using the forward solve from this course's help package
help(forward_solve)
Help on function forward_solve in module chen_3170.help:
forward_solve(l_mtrx, b_vec, loop_option='use-dot-product')
Performs a forward solve with a lower triangular matrix and right side vector.
Parameters
----------
l_mtrx: numpy.ndarray, required
Lower triangular matrix.
b_vec: numpy.ndarray, required
Right-side vector.
loop_option: string, optional
This is an internal option to demonstrate the usage of an explicit
double loop or an implicit loop using a dot product.
Default: 'use-dot-product'
Returns
-------
x_vec: numpy.narray
Solution vector returned.
Examples
--------
'''Usage example'''
x_vec = forward_solve( l_mtrx, b_vec )
np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)
x = [ 1. 0. -0.167]
'''View the source code in the notebook'''
!cat "chen_3170/help.py" # ugly but works for now
#!/usr/bin/env python
#--*-- coding: utf-8 -*-
# This file is part of the ChEn-3170 Computational Methods in Chemical Engineering
# course at https://github.com/dpploy/chen-3170
def get_triangular_matrix( mode='lower', ndim=None, mtrx=None ):
'''
Returns a triangular matrix in-place.
If a matrix is given, the function will modify the input, in place, into a
triangular matrix. The mtrx object will be modified and reflected on the callee side.
Otherwise, the function generates a random triangular matrix.
Parameters
----------
mode: string, optional
Type of triangular matrix: 'lower' or 'upper'. Defaults to lower
triangular.
ndim: int, optional
Dimension of the square matrix. If a matrix is not provided this
argument is required.
mtrx: numpy.ndarray, optional
square matrix to be turned into a triangular matrix.
Returns
-------
mtrx: numpy.ndarray
If a matrix was not passed the return is random array. If a matrix
was passed, its view is modified.
Examples
--------
>>> a_mtrx = ce.get_triangular_matrx('lower',3)
>>> a_mtrx
array([[0.38819556, 0. , 0. ],
[0.12304746, 0.07522054, 0. ],
[0.96357929, 0.69187941, 0.2878785 ]])
'''
assert ndim is None or mtrx is None, 'ndim or mtrx must be given; not both.'
assert not (ndim is None and mtrx is None), 'either ndim or mtrx must be given.'
assert mode =='lower' or mode =='upper', 'invalid mode %r.'%mode
if mtrx is None:
import numpy as np
mtrx = np.random.random((ndim,ndim))
else:
assert mtrx.shape[0] == mtrx.shape[1], 'matrix not square.'
# ready to return matrix
if mode == 'lower':
for i in range(mtrx.shape[0]):
mtrx[i,i+1:] = 0.0
elif mode == 'upper':
for j in range(mtrx.shape[1]):
mtrx[j+1:,j] = 0.0
else:
assert False, 'oops. something is very wrong.'
return mtrx
#*********************************************************************************
def forward_solve(l_mtrx, b_vec, loop_option='use-dot-product'):
'''
Performs a forward solve with a lower triangular matrix and right side vector.
Parameters
----------
l_mtrx: numpy.ndarray, required
Lower triangular matrix.
b_vec: numpy.ndarray, required
Right-side vector.
loop_option: string, optional
This is an internal option to demonstrate the usage of an explicit
double loop or an implicit loop using a dot product.
Default: 'use-dot-product'
Returns
-------
x_vec: numpy.narray
Solution vector returned.
Examples
--------
'''
import numpy as np
# sanity test
assert isinstance(l_mtrx,np.ndarray) # l_mtrx must be np.ndarray
assert l_mtrx.shape[0] == l_mtrx.shape[1],'non-square matrix.' # l_mtrx must be square
assert np.all(np.abs(np.diagonal(l_mtrx)) > 0.0),'zero value on diagonal.'
rows_ids, cols_ids = np.where(np.abs(l_mtrx) > 0) # get i, j of non zero entries
assert np.all(rows_ids >= cols_ids),'non-triangular matrix.' # test i >= j
assert b_vec.shape[0] == l_mtrx.shape[0],'incompatible l_mtrx @ b_vec dimensions' # b_vec must be compatible to l_mtrx
assert loop_option == 'use-dot-product' or loop_option == 'use-double-loop'
# end of sanity test
m_rows = l_mtrx.shape[0]
n_cols = m_rows
x_vec = np.zeros(n_cols)
if loop_option == 'use-dot-product':
for i in range(m_rows):
sum_lx = np.dot( l_mtrx[i,:i], x_vec[:i] )
#sum_lx = l_mtrx[i,:i] @ x_vec[:i] # matrix-vec mult. alternative to dot product
x_vec[i] = b_vec[i] - sum_lx
x_vec[i] /= l_mtrx[i,i]
elif loop_option == 'use-double-loop':
for i in range(m_rows):
sum_lx = 0.0
for j in range(i):
sum_lx += l_mtrx[i,j] * x_vec[j]
x_vec[i] = b_vec[i] - sum_lx
x_vec[i] /= l_mtrx[i,i]
else:
assert False, 'not allowed option: %r'%loop_option
return x_vec
#*********************************************************************************
def plot_matrix(mtrx, color_map='bw', title=None):
'''
Plot matrix as an image.
Parameters
----------
mtrx: numpy.ndarray, required
Matrix data.
color_map: str, optional
Color map for image: 'bw' black and white
title: str, optional
Title for plot.
Returns
-------
None:
Examples
--------
'''
# sanity check
import numpy as np
assert isinstance(mtrx,np.ndarray)
import numpy as np
from matplotlib import pyplot as plt # import the pyplot function of the matplotlib package
# end of sanity check
plt.rcParams['figure.figsize'] = [20, 4] # extend the figure size on screen output
plt.figure(1)
if color_map == 'bw':
plt.imshow(np.abs(mtrx),cmap='gray')
else:
plt.imshow(mtrx,cmap=color_map)
if title is not None:
plt.title(title,fontsize=14)
print('matrix shape =',mtrx.shape) # inspect the array shape
plt.show()
return
#*********************************************************************************
def print_reactions(reactions):
'''
Nice printout of a reactions list.
Parameters
----------
reactions: list(str)
Reactions in the form of a list.
Returns
-------
None:
Examples
--------
'''
# sanity check
assert isinstance(reactions,list)
# end of sanity check
for r in reactions:
i = reactions.index(r)
print('r%s'%i,': ',r)
n_reactions = len(reactions)
print('n_reactions =',n_reactions)
return
#*********************************************************************************
def print_reaction_sub_mechanisms( sub_mechanisms, mode=None, print_n_sub_mech=None ):
'''
Nice printout of a scored reaction sub-mechanism list
Parameters
----------
sub_mechanims: list(str), required
Sorted reaction mechanims in the form of a list.
mode: string, optional
Printing mode: all, top, None. Default: all
Returns
-------
None:
Examples
--------
'''
# sanity check
assert mode is None or print_n_sub_mech is None
assert mode =='top' or mode =='all' or mode==None
assert isinstance(print_n_sub_mech,int) or print_n_sub_mech is None
# end of sanity check
if mode is None and print_n_sub_mech is None:
mode = 'all'
if print_n_sub_mech is None:
if mode == 'all':
print_n_sub_mech = len(sub_mechanisms)
elif mode == 'top':
scores = [sm[3] for sm in sub_mechanisms]
max_score = max(scores)
tmp = list()
for s in scores:
if s == max_score:
tmp.append(s)
print_n_sub_mech = len(tmp)
else:
assert False, 'illegal mode %r'%mode
for rm in sub_mechanisms:
if sub_mechanisms.index(rm) > print_n_sub_mech-1: continue
print('Reaction Sub Mechanism: %s (score %4.2f)'%(sub_mechanisms.index(rm),rm[3]))
for (i,r) in zip( rm[0], rm[1] ):
print('r%s'%i,r)
return
#*********************************************************************************
def read_arrhenius_experimental_data(filename):
'''
Read k versus T data for fitting an Arrhenius rate constant expression.
Parameters
----------
filename: string, required
File name of data file including the path.
Returns
-------
r_cte: float
Universal gas constant.
r_cte: string
Universal gas constant unit.
n_pts: int
Number of data points
temp: np.ndarray, float
Temperature data.
k_cte: np.ndarray, float
Reaction rate constant data.
Examples
--------
'''
import io # import io module
finput = open(filename, 'rt') # create file object
import numpy as np
for line in finput:
line = line.strip()
if line[0] == '#': # skip comments in the file
continue
var_line = line.split(' = ')
if var_line[0] == 'r_cte':
r_cte = float(var_line[1].split(' ')[0])
r_cte_units = var_line[1].split(' ')[1]
elif var_line[0] == 'n_pts':
n_pts = int(var_line[1])
temp = np.zeros(n_pts)
k_cte = np.zeros(n_pts)
idx = 0 # counter
else:
data = line.split(' ')
temp[idx] = float(data[0])
k_cte[idx] = float(data[1])
idx += 1
return (r_cte, r_cte_units, n_pts, temp, k_cte)
#*********************************************************************************
def plot_arrhenius_experimental_data( temp, k_cte ):
'''
Plot T versus k data for fitting an Arrhenius rate constant expression.
Parameters
----------
temp: nd.array, required
Temperature data.
k_cte: nd.array, required
Reaction rate constant data.
Returns
-------
None: None
Examples
--------
'''
import matplotlib.pyplot as plt
plt.figure(1, figsize=(7, 7))
plt.plot(temp, k_cte,'r*',label='experimental')
plt.xlabel(r'$T$ [K]',fontsize=14)
plt.ylabel(r'$k$ [s$^{-1}$]',fontsize=14)
plt.title('Arrhenius Rxn Rate Constant Data',fontsize=20)
plt.legend(loc='best',fontsize=12)
plt.grid(True)
plt.show()
print('')
return
#*********************************************************************************
def color_map( num_colors ):
'''
Nice colormap for plotting.
Parameters
----------
num_colors: int, required
Number of colors.
Returns
-------
color_map: list(tuple(R,G,B,A))
List with colors interpolated from internal list of primary colors.
'''
assert num_colors >= 1
import numpy as np
# primary colors
# use the RGBA decimal code
red = np.array((1,0,0,1))
blue = np.array((0,0,1,1))
magenta = np.array((1,0,1,1))
green = np.array((0,1,0,1))
orange = np.array((1,0.5,0,1))
black = np.array((0,0,0,1))
yellow = np.array((1,1,0,1))
cyan = np.array((0,1,1,1))
# order the primary colors here
color_map = list()
color_map = [red, blue, orange, magenta, green, yellow, cyan, black]
num_primary_colors = len(color_map)
if num_colors <= num_primary_colors:
return color_map[:num_colors]
# interpolate primary colors
while len(color_map) < num_colors:
j = 0
for i in range(len(color_map)-1):
color_a = color_map[2*i]
color_b = color_map[2*i+1]
mid_color = (color_a+color_b)/2.0
j = 2*i+1
color_map.insert(j,mid_color) # insert before index
if len(color_map) == num_colors:
break
return color_map
#*********************************************************************************
def get_covid_19_us_data( type='deaths' ):
'''
Load COVID-19 pandemic cumulative data from:
https://github.com/CSSEGISandData/COVID-19.
Parameters
----------
type: str, optional
Type of data. Deaths ('deaths') and confirmed cases ('confirmed').
Default: 'deaths'.
Returns
-------
data: tuple(int, list(str), list(int))
(population, dates, cases)
'''
import pandas as pd
if type == 'deaths':
df = pd.read_csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv')
#df.to_html('covid_19_deaths.html')
elif type == 'confirmed':
df = pd.read_csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv')
df_pop = pd.read_csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv')
#df.to_html('covid_19_deaths.html')
#df.to_html('covid_19_confirmed.html')
else:
assert True, 'invalid query type: %r (valid: "deaths", "confirmed"'%(type)
df = df.drop(['UID','iso2','iso3','Combined_Key','code3','FIPS','Lat', 'Long_','Country_Region'],axis=1)
df = df.rename(columns={'Province_State':'state/province','Admin2':'city'})
import numpy as np
state_names = list()
state_names_tmp = list()
for (i,istate) in enumerate(df['state/province']):
if istate.strip() == 'Wyoming' and df.loc[i,'city']=='Weston':
break
state_names_tmp.append(istate)
state_names_set = set(state_names_tmp)
state_names = list(state_names_set)
state_names = sorted(state_names)
dates = np.array(list(df.columns[3:]))
population = [0]*len(state_names)
cases = np.zeros( (len(df.columns[3:]),len(state_names)), dtype=np.float64)
for (i,istate) in enumerate(df['state/province']):
if istate.strip() == 'Wyoming' and df.loc[i,'city']=='Weston':
break
state_id = state_names.index(istate)
if type == 'confirmed':
population[state_id] += int(df_pop.loc[i,'Population'])
else:
population[state_id] += int(df.loc[i,'Population'])
cases[:,state_id] += np.array(list(df.loc[i, df.columns[3:]]))
return ( state_names, population, dates, cases )
#*********************************************************************************
def get_covid_19_global_data( type='deaths', distribution=True, cumulative=False ):
'''
Load COVID-19 pandemic cumulative data from:
https://github.com/CSSEGISandData/COVID-19
Parameters
----------
type: str, optional
Type of data. Deaths ('deaths') and confirmed cases ('confirmed').
Default: 'deaths'.
distribution: bool, optional
Distribution of new cases over dates.
Default: True
cumulative: bool, optional
Cumulative number of cases over dates.
Default: False
Returns
-------
data: tuple(int, list(str), list(int))
(contry_names, dates, cases)
'''
if cumulative is True:
distribution = False
import pandas as pd
if type == 'deaths':
df = pd.read_csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv')
#df.to_html('covid_19_global_deaths.html')
else:
assert True, 'invalid query type: %r (valid: "deaths"'%(type)
df = df.drop(['Lat', 'Long'],axis=1)
df = df.rename(columns={'Province/State':'state/province','Country/Region':'country/region'})
import numpy as np
country_names = list()
country_names_tmp = list()
for (i,icountry) in enumerate(df['country/region']):
country_names_tmp.append(icountry)
country_names_set = set(country_names_tmp)
country_names = list(country_names_set)
country_names = sorted(country_names)
dates = np.array(list(df.columns[2:]))
cases = np.zeros( (len(df.columns[2:]),len(country_names)), dtype=np.float64)
for (i,icountry) in enumerate(df['country/region']):
country_id = country_names.index(icountry)
cases[:,country_id] += np.array(list(df.loc[i, df.columns[2:]]))
if distribution:
for j in range(cases.shape[1]):
cases[:,j] = np.round(np.gradient( cases[:,j] ),0)
return ( country_names, dates, cases )
#*********************************************************************************
def forward_solve(l_mtrx, b_vec, loop_option='use-dot-product'):
'''
Performs a forward solve with a lower triangular matrix and right side vector.
Parameters
----------
l_mtrx: numpy.ndarray, required
Lower triangular matrix.
b_vec: numpy.ndarray, required
Right-side vector.
loop_option: string, optional
This is an internal option to demonstrate the usage of an explicit
double loop or an implicit loop using a dot product.
Default: 'use-dot-product'
Returns
-------
x_vec: numpy.narray
Solution vector returned.
Examples
--------
'''
import numpy as np
# sanity test
assert isinstance(l_mtrx,np.ndarray) # l_mtrx must be np.ndarray
assert l_mtrx.shape[0] == l_mtrx.shape[1],'non-square matrix.' # l_mtrx must be square
assert np.all(np.abs(np.diagonal(l_mtrx)) > 0.0),'zero value on diagonal.'
rows_ids, cols_ids = np.where(np.abs(l_mtrx) > 0) # get i, j of non zero entries
assert np.all(rows_ids >= cols_ids),'non-triangular matrix.' # test i >= j
assert b_vec.shape[0] == l_mtrx.shape[0],'incompatible l_mtrx @ b_vec dimensions' # b_vec must be compatible to l_mtrx
assert loop_option == 'use-dot-product' or loop_option == 'use-double-loop'
# end of sanity test
m_rows = l_mtrx.shape[0]
n_cols = m_rows
x_vec = np.zeros(n_cols)
if loop_option == 'use-dot-product':
for i in range(m_rows):
sum_lx = np.dot( l_mtrx[i,:i], x_vec[:i] )
#sum_lx = l_mtrx[i,:i] @ x_vec[:i] # matrix-vec mult. alternative to dot product
x_vec[i] = b_vec[i] - sum_lx
x_vec[i] /= l_mtrx[i,i]
elif loop_option == 'use-double-loop':
for i in range(m_rows):
sum_lx = 0.0
for j in range(i):
sum_lx += l_mtrx[i,j] * x_vec[j]
x_vec[i] = b_vec[i] - sum_lx
x_vec[i] /= l_mtrx[i,i]
else:
assert False, 'not allowed option: %r'%loop_option
return x_vec
A lower triangular matrix allows for a forward solve. The algorithm for $\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} \quad\ \forall \quad\ 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.
$\Umtrx$ backward solve¶
A upper triangular matrix allows for a backward solve. The algorithm for $\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} \quad\ \forall \quad\ 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.
'''U backward solve'''
u_mtrx = np.array( [[1., 2., 3.], # per course notes
[0, 4., 5.],
[0., 0., 6.]] )
b_vec = np.array( [1.,2.,3.] )
try:
from chen_3170.toolkit import backward_solve
except ModuleNotFoundError:
assert False, 'You need to provide your own backward_solve function here. Bailing out.'
x_vec = backward_solve( u_mtrx, b_vec )
np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)
x = [-0.25 -0.125 0.5 ]
$\Amtrx = \Lmtrx\,\Umtrx$ factorization¶
$\Lmtrx\,\Umtrx$ factorization algorithm (without using pivoting) for a square matrix $\overset{(m \times m)}{\Amtrx}$ computes the $\Lmtrx\,\Umtrx$ factors. 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\ \forall\ i=k+1,\ldots,m \quad\ \text{and} \quad\ j=k,\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$.
'''L U factors of A'''
a_mtrx = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 10.]] )
try:
from chen_3170.toolkit import lu_factorization
except ModuleNotFoundError:
assert False, 'You need to provide your own lu_factorization function here. Bailing out.'
(l_mtrx, u_mtrx, rank) = lu_factorization( a_mtrx )
print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('rank =',rank)
print('')
print('Checking...')
print('LU - A =\n', l_mtrx @ u_mtrx - a_mtrx)
L = [[1. 0. 0.] [4. 1. 0.] [7. 2. 1.]] U = [[ 1. 2. 3.] [ 0. -3. -6.] [ 0. 0. 1.]] rank = 3 Checking... LU - A = [[0. 0. 0.] [0. 0. 0.] [0. 0. 0.]]
$\Pmtrx\,\Amtrx = \Lmtrx\,\Umtrx$ factorization¶
The factors: $\Pmtrx$, $\Lmtrx$, and $\Umtrx$ where $\Lmtrx\,\Umtrx = \Pmtrx\,\Amtrx$ can be obtained from partial pivoting strategy to the $\Lmtrx\,\Umtrx$ factorization algorithm shown above. $\Pmtrx$ is a row permutation matrix obtained by the underlying Gaussian elimination used to construct the $\Lmtrx$ and $\Umtrx$ factors.
Program a $\Lmtrx\,\Umtrx$ factorization algorithm (using partial pivoting) for a square matrix $\overset{(m \times m)}{\Amtrx}$ and compute the $\Pmtrx\,\Lmtrx\,\Umtrx$ factors.
The factorization is obtained by elimination steps $k = 1,\ldots,m-1$ so that $A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\, m_{i,k} \ \forall\ i=k+1,\ldots,m \ \text{and}\ j=k,\ldots,m$ 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$. However, every $k$-step selects a pivot $A^{(k)}_{k,k}$ of maximum absolute value via row exchanges recorded in the permutation matrix $\Pmtrx$.
'''P L U factors of A'''
a_mtrx = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 10.]] )
try:
from chen_3170.toolkit import lu_factorization
except ModuleNotFoundError:
assert False, 'You need to provide your own lu_factorization function here. Bailing out.'
(p_mtrx, l_mtrx, u_mtrx, rank) = lu_factorization( a_mtrx, pivoting_option='partial' )
print('P =\n',p_mtrx)
print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('rank =',rank)
print('')
print('Checking...')
print('LU - PA =\n', l_mtrx @ u_mtrx - p_mtrx @ a_mtrx)
P = [[0. 0. 1.] [1. 0. 0.] [0. 1. 0.]] L = [[1. 0. 0. ] [0.143 1. 0. ] [0.571 0.5 1. ]] U = [[ 7.000e+00 8.000e+00 1.000e+01] [ 0.000e+00 8.571e-01 1.571e+00] [ 0.000e+00 5.551e-17 -5.000e-01]] rank = 3 Checking... LU - PA = [[0. 0. 0.] [0. 0. 0.] [0. 0. 0.]]
$\Pmtrx\,\Amtrx\,\Qmtrx = \Lmtrx\,\Umtrx$ factorization¶
The factors: $\Pmtrx$, $\Qmtrx$, $\Lmtrx$, and $\Umtrx$ where $\Lmtrx\,\Umtrx = \Pmtrx\,\Amtrx\,\Qmtrx$ can be obtained from a user-developed algorithm with complete pivoting. $\Pmtrx$ is a row permutation matrix and $\Qmtrx$ is a column permutation matrix, obtained by the underlying Gaussian elimination used to construct the $\Lmtrx$ and $\Umtrx$ factors.
Program a $\Lmtrx\,\Umtrx$ factorization algorithm (using complete pivoting) for a square matrix $\overset{(m \times m)}{\Amtrx}$ and compute the $\Pmtrx\,,\Qmtrx\,,\Lmtrx\,,\Umtrx$ factors. The factorization is obtained by elimination steps $k = 1,\ldots,m-1$ so that $A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\, m_{i,k} \ \forall\ i=k+1,\ldots,m \ \text{and}\ j=k,\ldots,m$ 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$. However, every $k$-step selects a pivot $A^{(k)}_{k,k}$ of maximum absolute value via row exchanges recorded in the permutation matrix $\Pmtrx$ and column exchanges recorded in the permutation matrix $\Qmtrx$.
'''P Q L U factors of A'''
a_mtrx = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 10.]] )
try:
from chen_3170.toolkit import lu_factorization
except ModuleNotFoundError:
assert False, 'You need to provide your own lu_factorization function here. Bailing out.'
(p_mtrx, q_mtrx, l_mtrx, u_mtrx, rank) = lu_factorization( a_mtrx, pivoting_option='complete' )
print('P =\n',p_mtrx)
print('Q =\n',q_mtrx)
print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('rank =',rank)
print('')
print('Checking...')
print('LU - PAQ =\n', l_mtrx @ u_mtrx - p_mtrx @ a_mtrx @ q_mtrx)
P = [[0. 0. 1.] [1. 0. 0.] [0. 1. 0.]] Q = [[0. 1. 0.] [0. 0. 1.] [1. 0. 0.]] L = [[1. 0. 0. ] [0.3 1. 0. ] [0.6 0.182 1. ]] U = [[10. 7. 8. ] [ 0. -1.1 -0.4 ] [ 0. 0. 0.273]] rank = 3 Checking... LU - PAQ = [[0. 0. 0.] [0. 0. 0.] [0. 0. 0.]]