ChEn-3170: Computational Methods in Chemical Engineering Spring 2024 UMass Lowell; Prof. V. F. de Almeida 03Feb24
$ \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{\vvec}{\boldsymbol{\mathsf{v}}} \newcommand{\cvec}{\boldsymbol{\mathsf{c}}} \newcommand{\rvec}{\boldsymbol{\mathsf{r}}} \newcommand{\norm}[1]{\bigl\lVert{#1}\bigr\rVert} \DeclareMathOperator{\rank}{rank} $
Numpy
set_printoptions()
methodNumpy
set_printoptions()
method'''Python packages are accessed with an import directive as such:'''
import numpy as np # import the package and create the alias: np
import math
pi = math.pi
print('pi = %15.5e'%pi) # formatting numeric output: scientific notation
pi = 3.14159e+00
print('pi = %10.5f'%pi) # formatting numeric output: float
pi = 3.14159
print('pi = %10.5e and e = %8.3f'%(pi, math.e)) # formatting numeric output: sci. notation and float
pi = 3.14159e+00 and e = 2.718
Numpy
arrays'''Use set print options in Numpy'''
np.set_printoptions(precision=4, threshold=800, edgeitems=6, linewidth=105)
mtrx = np.random.random((10,100))
print(mtrx)
[[8.4138e-02 3.9647e-01 4.6171e-01 9.4997e-01 2.8638e-01 8.7224e-02 ... 4.0240e-01 4.0094e-01 8.5229e-01 8.8223e-01 7.5218e-01 4.0944e-01] [1.4907e-01 4.3748e-01 7.9740e-01 3.0196e-01 7.5468e-01 8.1755e-01 ... 8.5433e-02 9.0723e-01 7.3662e-01 2.6948e-01 8.9656e-01 4.8746e-01] [7.7662e-01 7.9629e-01 8.4358e-01 5.0597e-01 1.4333e-01 5.6789e-01 ... 2.3390e-01 2.6261e-01 4.4524e-01 4.7857e-01 8.0243e-01 9.4277e-01] [6.7695e-01 9.7351e-01 6.1864e-01 5.6811e-02 4.6991e-01 8.9605e-01 ... 8.3869e-01 5.1893e-01 7.0219e-01 6.7778e-01 6.9872e-01 3.0444e-01] [6.3906e-01 4.4086e-01 2.3749e-01 9.3303e-01 8.8561e-01 4.4480e-01 ... 6.9081e-01 5.6960e-01 7.5011e-01 5.2257e-01 1.0237e-01 4.7495e-01] [4.6809e-02 6.9367e-02 1.7840e-01 1.0421e-01 7.6502e-01 5.6952e-01 ... 6.7518e-01 7.6232e-01 2.8356e-01 4.5567e-01 6.2981e-01 9.7600e-01] [5.0889e-01 7.8494e-02 8.9570e-01 1.7171e-01 9.8414e-01 1.3523e-01 ... 8.1665e-01 9.4393e-02 3.4677e-01 6.0295e-01 8.8307e-01 4.5862e-04] [9.2530e-01 5.5120e-01 6.3884e-01 3.1400e-01 4.6168e-01 3.4693e-01 ... 5.5343e-01 1.1361e-01 7.4256e-01 1.4876e-01 1.1083e-02 8.3941e-01] [6.9449e-01 3.5907e-02 8.1853e-01 6.4938e-01 6.3062e-01 3.9009e-01 ... 7.6055e-01 4.3796e-02 9.0506e-01 9.6696e-01 6.3922e-01 6.3121e-01] [9.4677e-01 2.7162e-01 5.5534e-01 4.7331e-01 5.6053e-01 1.4532e-01 ... 1.0760e-01 8.6322e-01 3.2204e-01 4.5782e-01 1.1872e-01 1.0333e-01]]
#help(np.set_printoptions)
In all of engineering calculations use double-precision floating point numeric
'''Set double precision at creation time'''
x_vec = np.empty(10, dtype=np.float64)
print(type(x_vec))
print(x_vec.dtype)
<class 'numpy.ndarray'> float64
'''Set double precision after creation'''
x_vec = x_vec.astype(float)
print(type(x_vec))
print(x_vec.dtype)
<class 'numpy.ndarray'> float64
'''Set single precision after creation; not to be used'''
x_vec = x_vec.astype(np.float32)
print(type(x_vec))
print(x_vec.dtype)
<class 'numpy.ndarray'> float32
'''Element-by-element addition or subtraction'''
vec1 = np.array(np.random.random(5))
print('vec1 =',vec1)
vec2 = np.array(np.random.random(5))
print('vec2 =',vec2)
result = vec1 + vec2 # element-by-element sum
print('addition =',result)
result = vec1 - vec2 # element-by-element subtraction
print('difference =',result)
vec1 = [0.5453 0.4583 0.2304 0.7032 0.37 ] vec2 = [0.506 0.9526 0.5948 0.4733 0.0357] addition = [1.0513 1.4109 0.8251 1.1765 0.4058] difference = [ 0.0393 -0.4943 -0.3644 0.2299 0.3343]
'''Element-by-element product or division'''
vec1 = np.array(np.random.random(5))
print('vec1 =',vec1)
vec2 = np.array(np.random.random(5))
print('vec2 =',vec2)
result = vec1 * vec2 # element-by-element product
print('product =',result)
result = vec1 / vec2 # element-by-element division
print('division =',result)
vec1 = [0.2579 0.1909 0.0031 0.342 0.1954] vec2 = [0.2264 0.3046 0.066 0.573 0.6918] product = [0.0584 0.0581 0.0002 0.196 0.1352] division = [1.1393 0.6266 0.0466 0.5968 0.2824]
'''Product of all elements of a vector'''
vec1_prod = np.prod(vec1)
print('vec1 =', vec1)
print('vec1 product =', vec1_prod)
vec1 = [0.2579 0.1909 0.0031 0.342 0.1954] vec1 product = 1.0105005747408923e-05
The result of the inner product of two vectors: $\vvec_1 \cdot \vvec_2$ is a scalar.
'''Vector inner product or dot product'''
vec1 = np.array(np.random.random(5))
print('vec1 =',vec1)
vec2 = np.array(np.random.random(5))
print('vec2 =',vec2)
result = np.dot(vec1, vec2) # inner or dot product
print('dot product =',result)
vec1 = [0.3084 0.995 0.775 0.5653 0.4646] vec2 = [0.2284 0.3363 0.3029 0.1459 0.4418] dot product = 0.9275588056094912
'''More on vector inner product or dot product'''
'''Another way to compute the inner product'''
ele_by_ele_product_vec = vec1 * vec2
inner_product = ele_by_ele_product_vec.sum()
print('vec1 . vec2 = ', inner_product)
vec1 . vec2 = 0.9275588056094912
The magnitude of a vector, $\norm{\vvec}$, is related to its inner product: $\vvec \cdot \vvec = \norm{\vvec}^2$ which is a scalar.
'''Norm of a vector or its magnitude'''
vec = np.array(np.random.random(5))
vec_norm = np.sqrt(np.dot(vec,vec))
print('v = ', vec)
print('||v|| = ', vec_norm)
print('np.linalg.norm ||v|| = ', np.linalg.norm(vec))
v = [0.343 0.1219 0.255 0.9625 0.8877] ||v|| = 1.3827065468627124 np.linalg.norm ||v|| = 1.3827065468627124
'''Scaling of a vector'''
vec = np.array(np.random.random(5))
print('vec =',vec)
factor = 0.345
scaled = factor * vec # scaling of vec element-by-element product
print('scaled =', scaled) # assigned to new variable `scaled`
vec *= factor # in-place scaling
print('vec =',vec)
vec = [0.0205 0.2177 0.3561 0.5198 0.9982] scaled = [0.0071 0.0751 0.1229 0.1793 0.3444] vec = [0.0071 0.0751 0.1229 0.1793 0.3444]
'''Mathematical Operations on a Vector'''
vec = np.array(np.random.random(5))
print('vec =',vec)
log_vec = np.log(vec) # natural log element-by-element
print('log(vec) =',log_vec)
exp_vec = np.exp(log_vec) # exponential
print('exp(vec) =',exp_vec)
sin_vec = np.sin(vec) # sine
print('sin(vec) =',sin_vec)
vec_cubed = vec**3 # powers
print('vec^3 =',vec_cubed)
vec_mean = vec.mean() # arithmetic mean
print('mean(vec) =',vec_mean)
vec_std = vec.std() # standard deviation
print('std(vec) =',vec_std)
vec = [0.6139 0.7863 0.434 0.8923 0.863 ] log(vec) = [-0.4879 -0.2404 -0.8346 -0.114 -0.1473] exp(vec) = [0.6139 0.7863 0.434 0.8923 0.863 ] sin(vec) = [0.5761 0.7078 0.4205 0.7785 0.7598] vec^3 = [0.2314 0.4862 0.0818 0.7104 0.6428] mean(vec) = 0.7179130660462255 std(vec) = 0.17178674000332048
'''Searching a vector for entries matching a test'''
# what are the indices of the values in "vec" that satisfy: vec[] >= 0.3
(idx_ids, ) = np.where(vec >= 0.3)
print('vec =', vec)
print('indices = ', idx_ids)
vec = [0.6139 0.7863 0.434 0.8923 0.863 ] indices = [0 1 2 3 4]
'''Searching a vector for entries matching a test'''
# what are the indices of the values in "vec" that satisfy: vec[] == 0.3
(idx_ids, ) = np.where(vec == 0.3)
print('vec =', vec)
print('indices = ', idx_ids)
vec = [0.6139 0.7863 0.434 0.8923 0.863 ] indices = []
'''Zip creates a list of tuples on the fly'''
print(list(zip(vec1, vec2)))
[(0.3083937089424802, 0.22841714407520286), (0.9949698101315914, 0.3363147949745209), (0.7749536526133711, 0.30294856337531473), (0.5652559409417851, 0.14588439346608584), (0.4645754380083772, 0.4418232111080259)]
In all of engineering calculations use double-precision floating point numeric
'''Set double precision at creation time'''
mtrx = np.empty((5,5), dtype=np.float64)
print(type(mtrx))
print(mtrx.dtype)
<class 'numpy.ndarray'> float64
'''Set double precision after creation'''
mtrx = mtrx.astype(float)
print(type(mtrx))
print(mtrx.dtype)
<class 'numpy.ndarray'> float64
'''Element-by-element addition or subtraction'''
mat1 = np.random.random((3,3))
print('mat1 =\n', mat1)
mat2 = np.random.random((3,3))
print('mat2 =\n', mat2)
result = mat1 + mat2 # element-by-element sum
print('addition =\n', result)
result = mat1 - mat2 # element-by-element subtraction
print('difference =\n', result)
mat1 = [[0.9423 0.4184 0.9301] [0.1865 0.1441 0.3504] [0.5401 0.9425 0.8513]] mat2 = [[0.6776 0.2725 0.6021] [0.884 0.1885 0.1843] [0.2853 0.787 0.0605]] addition = [[1.6199 0.6909 1.5321] [1.0705 0.3327 0.5346] [0.8254 1.7296 0.9118]] difference = [[ 0.2646 0.1458 0.328 ] [-0.6975 -0.0444 0.1661] [ 0.2548 0.1555 0.7908]]
'''Element-by-element product or division'''
mat1 = np.random.random((3,3))
print('mat1 =\n', mat1)
mat2 = np.random.random((3,3))
print('mat2 =\n', mat2)
result = mat1 * mat2 # element-by-element product
print('product =\n', result)
result = mat1 / mat2 # element-by-element division (cross your fingers)
print('division =\n', result)
mat1 = [[0.8177 0.1258 0.3875] [0.6726 0.7571 0.0418] [0.2121 0.9535 0.462 ]] mat2 = [[0.3448 0.7567 0.7727] [0.4003 0.7337 0.2174] [0.5446 0.386 0.8489]] product = [[0.282 0.0952 0.2994] [0.2693 0.5555 0.0091] [0.1155 0.368 0.3922]] division = [[2.3716 0.1663 0.5015] [1.6802 1.0318 0.1921] [0.3894 2.4705 0.5442]]
'''Produce Noise on a Matrix Image (brick data)'''
from matplotlib import pyplot as plt # import the pyplot function of the matplotlib package
#%matplotlib inline
plt.rcParams['figure.figsize'] = [20, 4] # extend the figure size on screen output
# Read image from the images/ directory in the chen-3170 repo
block = plt.imread('images/glacier.png', format='png')
#block = plt.imread('https://raw.githubusercontent.com/dpploy/chen-3170/master/notebooks/images/cp21.8-T.png', format='png')
plt.figure(1)
plt.imshow(block)
plt.title('Matrix Reloaded', fontsize=14)
plt.show()
print('block shape =', block.shape) # inspect the array shape
block shape = (300, 400, 3)
'''Use Matrix Element-by-Element Multiplication'''
mtrx_shape = block.shape[0:2] # use the shape to automate noise_mtrx generation
noise_mtrx = np.random.random(mtrx_shape) # generate random matrix
mtrx_noise = block[:,:,2] * noise_mtrx # apply noise to the blue channel
plt.figure(2)
plt.imshow(mtrx_noise, cmap='gray')
plt.title('Noisy Matrix',fontsize=14)
plt.show()
mtrx_noise.shape
(300, 400)
'''Matrix Scaling (matrix product or division by a scalar)'''
mat1 = np.random.random((3,3))
print('mat1 =\n',mat1)
factor = 3.21
result = factor * mat1 # scaling of mat1 element-by-element; product with factor
print('scaled =\n',result)
mat1 = [[0.9643 0.5686 0.6097] [0.1848 0.1046 0.026 ] [0.0069 0.0106 0.7162]] scaled = [[3.0953 1.8252 1.9571] [0.5932 0.3357 0.0836] [0.022 0.0341 2.2989]]
'''Matrix Scaling of an Image'''
color_channel = np.copy(block[:,:,0]) # copy the red channel
color_channel /= color_channel.max() # scale to gray, 0-255 values
color_channel *= 255
gray_channel = color_channel.astype(np.uint8) # truncate all float data type to unsigned int 8
plt.figure(3)
plt.imshow(gray_channel, cmap='gray')
#plt.imshow(gray_channel)
plt.title('Matrix Gray Scaling',fontsize=14)
plt.show()
'''Other Mathematical Operations on a Matrix'''
mtrx = np.copy(block[:,:,0]) # copy the red channel
plt.figure(4)
plt.imshow(mtrx, cmap='gray') # show channel as a flat image with default colormap
plt.title('Original', fontsize=14)
plt.show()
mtrx_mean = mtrx.mean() # arithmetic mean
print('mean(mtrx) =', mtrx_mean)
mtrx_std = mtrx.std() # standard deviation
print('std(mtrx) =', mtrx_std)
mean(mtrx) = 0.38594216 std(mtrx) = 0.37984067
'''Other Mathematical Operations on a Matrix'''
log_mtrx = np.log(mtrx + .001) # natural log element-by-element
plt.figure(5)
plt.imshow(log_mtrx, cmap='gray')
plt.title('Log Transform', fontsize=14)
plt.show()
mtrx_mean = log_mtrx.mean() # arithmetic mean
print('mean(mtrx) =', mtrx_mean)
mtrx_std = log_mtrx.std() # standard deviation
print('std(mtrx) =', mtrx_std)
mean(mtrx) = -2.6418166 std(mtrx) = 2.5803103
'''Other Mathematical Operations on a Matrix'''
exp_mtrx = np.exp(log_mtrx) # exponential
plt.figure(6)
plt.imshow(exp_mtrx, cmap='gray')
plt.title('Exp of Log Transform', fontsize=14)
plt.show()
mtrx_mean = exp_mtrx.mean() # arithmetic mean
print('mean(mtrx) =', mtrx_mean)
mtrx_std = exp_mtrx.std() # standard deviation
print('std(mtrx) =', mtrx_std)
mean(mtrx) = 0.38694212 std(mtrx) = 0.37984067
'''Other Mathematical Operations on a Matrix'''
sin_mtrx = np.sin(mtrx + np.pi/2) # sine
plt.figure(7)
plt.imshow(sin_mtrx)
plt.title('Sine Transform', fontsize=14)
plt.show()
mtrx_mean = sin_mtrx.mean() # arithmetic mean
print('mean(mtrx) =', mtrx_mean)
mtrx_std = sin_mtrx.std() # standard deviation
print('std(mtrx) =', mtrx_std)
mean(mtrx) = 0.86122394 std(mtrx) = 0.1534877
'''Other Mathematical Operations on a Matrix'''
mtrx_cubed = mtrx**3 # powers
plt.figure(8)
plt.imshow(mtrx_cubed, cmap='gray')
plt.title('Cube Transform', fontsize=14)
plt.show()
mtrx_mean = mtrx_cubed.mean() # arithmetic mean
print('mean(mtrx) =', mtrx_mean)
mtrx_std = mtrx_cubed.std() # standard deviation
print('std(mtrx) =', mtrx_std)
mean(mtrx) = 0.2348578 std(mtrx) = 0.28249925
'''Matrix Transposition'''
'''clockwise rotation followed by horizontal right to left flip'''
mtrx = np.random.random((5,7))
np.set_printoptions(precision=3,threshold=20,edgeitems=12,linewidth=100) # one way to control printing of numpy arrays
print('mtrx =\n',mtrx)
mtrx_T = mtrx.transpose() # transpose of a mtrx: M[i,j] -> M[j,i]
print('mtrx^T =\n', mtrx_T)
mtrx = [[0.741 0.772 0.778 0.153 0.986 0.249 0.908] [0.079 0.596 0.263 0.201 0.448 0.361 0.885] [0.636 0.236 0.825 0.408 0.424 0.899 0.616] [0.851 0.251 0.732 0.151 0.694 0.183 0.594] [0.728 0.082 0.81 0.719 0.853 0.57 0.087]] mtrx^T = [[0.741 0.079 0.636 0.851 0.728] [0.772 0.596 0.236 0.251 0.082] [0.778 0.263 0.825 0.732 0.81 ] [0.153 0.201 0.408 0.151 0.719] [0.986 0.448 0.424 0.694 0.853] [0.249 0.361 0.899 0.183 0.57 ] [0.908 0.885 0.616 0.594 0.087]]
'''Matrix Transposition'''
'''Example of adding a transformed matrix to another transform transposed'''
'''note: to add a matrix to its transpose, a matrix must be square'''
n_rows = block.shape[0]
n_columns = n_rows
mtrx = np.copy(block[:n_rows,:n_columns,0]) # select a square block; red channel
sin_mtrx = np.sin(mtrx + np.pi/2) # sine
sin_mtrx /= sin_mtrx.max()
plt.figure(9)
plt.imshow(sin_mtrx)
plt.title('Sine Transform', fontsize=14)
plt.show()
mtrx_cubed = mtrx**3 # powers
plt.figure(10)
plt.imshow(mtrx_cubed)
plt.title('Cube Transform', fontsize=14)
plt.show()
plt.figure(11)
plt.imshow(sin_mtrx + mtrx_cubed.transpose()) # sine + cubed transposed
plt.title('Sine + Cube Transpose Transform', fontsize=14)
plt.show()
'''Searching a matrix for entries matching a test'''
# what are the indices of the values in "mtrx" that satisfy: mtrx >= 0.3
(idx_ids, jdx_ids) = np.where(mtrx >= 0.3)
np.set_printoptions(precision=3, threshold=20, edgeitems=5, linewidth=100)
print('matrix =\n', mtrx)
print('ith indices = ',idx_ids)
print('jth indices = ',jdx_ids)
matrix = [[0.329 0.329 0.329 0.329 0.329 ... 0.384 0.384 0.384 0.384 0.384] [0.329 0.329 0.329 0.329 0.329 ... 0.388 0.388 0.388 0.388 0.388] [0.345 0.345 0.345 0.345 0.345 ... 0.392 0.392 0.392 0.392 0.392] [0.357 0.357 0.357 0.357 0.357 ... 0.4 0.4 0.4 0.4 0.4 ] [0.361 0.361 0.361 0.361 0.361 ... 0.408 0.408 0.408 0.408 0.408] ... [0.004 0.004 0.004 0.004 0.004 ... 0.004 0.004 0.004 0.004 0.004] [0.004 0.004 0.004 0.004 0.004 ... 0.012 0. 0. 0. 0. ] [0.004 0.004 0.004 0.004 0.004 ... 0.012 0. 0. 0. 0. ] [0.004 0.004 0.004 0.004 0.004 ... 0.012 0. 0. 0. 0. ] [0.004 0.004 0.004 0.004 0.004 ... 0.012 0. 0. 0. 0. ]] ith indices = [ 0 0 0 0 0 ... 252 252 252 252 252] jth indices = [ 0 1 2 3 4 ... 171 172 173 174 175]
'''Verify the searched elements'''
mtrx[idx_ids, jdx_ids].min()
0.3019608
mtrx[idx_ids, jdx_ids]
array([0.329, 0.329, 0.329, 0.329, 0.329, ..., 0.831, 0.737, 0.667, 0.58 , 0.459], dtype=float32)
mtrx[idx_ids, jdx_ids].shape
(54533,)
mtrx[idx_ids, jdx_ids].dtype
dtype('float32')