
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import cartopy
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from sympy import linsolve, symbols
import warnings
import os
from matplotlib.ticker import FuncFormatter
warnings.filterwarnings("ignore")


# Step 4: Solution to the mathematical model equations
A = np.array([1, 8, 0, -1, 0, 6, 1, -3, 3]).reshape(3,3, order='F')
b = np.array([0, 10, 15])
x = np.linalg.solve(A,b)
print(f"Solution to Ax = b: x = {x}")

Solution to Ax = b: x = [1.5        2.16666667 0.66666667]


# Step 5: Interpretation of the Modeling Results
# #Verify the solution: a*c = b 
print(f"A@x = {np.round(A@x)}")

A@x = [ 0. 10. 15.]


I = np.identity(3) #3-dim identity matrix
A = np.array([[0.3, 0.2, 0.2],[0.2, 0.2, 0.1],[0.3, 0.2, 0.3]])
D = np.array([3000, 500, 1500])

x1 = np.linalg.solve(I-A, D)
print(f"Solution to (I-A)x1 = D: x1 = {x1}")

Solution to (I-A)x1 = D: x1 = [6875.         3090.27777778 5972.22222222]


I = np.identity(3) #3-dim identity matrix
A = np.array([[0.4102, 0.0301, 0.0257], 
              [0.0624, 0.3783, 0.1050],
             [ 0.1236, 0.1588, 0.1919]])
D = np.array([39.24,60.02,130.65])

x1 = np.linalg.solve(I-A, D)
print(f"Solution to (I-A)x1 = D: x1 = {x1}")

Solution to (I-A)x1 = D: x1 = [ 82.40040843 138.85496907 201.56522654]


#Demonstrate SVD for a simple 2X3 matrix 
A1 = np.array([[1, 0, -1],[1, -1, 0]])
A1

array([[ 1,  0, -1],
       [ 1, -1,  0]])


U, DE, V = np.linalg.svd(A1, full_matrices=False)
U

array([[-0.70710678, -0.70710678],
       [-0.70710678,  0.70710678]])

V

array([[-8.16496581e-01,  4.08248290e-01,  4.08248290e-01],
       [ 1.19879883e-16, -7.07106781e-01,  7.07106781e-01]])

DE

array([1.73205081, 1.        ])


D = np.diag(DE )#forms the SVD diagonal matrix 
D

array([[1.73205081, 0.        ],
       [0.        , 1.        ]])


#Verification of SVD: A = UDV'
#np.round() rounds floats to the nearest integer
#This yields the original data matrix a1
np.round(U@D@V)

array([[ 1., -0., -1.],
       [ 1., -1.,  0.]])


#Graphically show the U column vectors, aka EOFs
#Figure 5.2
fig, ax = plt.subplots(1,2, figsize=(12,6))

ax[0].plot(np.arange(U.shape[1])+1, U[:,1],
           'bo-', label = "EOF 1", linewidth = 5)
ax[0].plot(np.arange(U.shape[1])+1, U[:,0],
           'ko-', label = "EOF 2", linewidth = 5)

ax[0].set_ylim(-1, 1)
ax[0].text(1.4, -0.3, "EOF 1", color = 'blue', size = 15,fontweight="bold")
ax[0].text(1.4, -0.65, "EOF 2", color = 'black', size = 15,fontweight="bold")

ax[0].set_xlabel("Spatial Position $x$", size = 15)
ax[0].set_ylabel("Mode Values [Dimensionless]", size = 15)
ax[0].set_title("Spatial Modes: EOFs",fontweight="bold")
ax[0].tick_params(labelsize=15)
ax[0].set_xticks([1, 1.25, 1.5, 1.75, 2])
ax[0].set_xticklabels([1, 1.25, 1.5, 1.75, 2])
# set_xticks allows you to choose where ticks show coordinate-wise
# along the x-axis, set_xticklabels sets the labels associated with
# those ticks, sequentially. y-axis counterparts also exist:
# set_yticks, set_yticklabels


ax[1].plot(np.arange(V.shape[1])+1, V[0,:],
           'go-', label = "PC1", linewidth = 5)
ax[1].plot(np.arange(V.shape[1])+1, V[1,:],
           'ro-', label = "PC2", linewidth = 5)

ax[1].set_ylim(-1, 1)
ax[1].text(1.9, 0, "PC 1", color = 'green', size = 15,fontweight="bold")
ax[1].text(1.9, -0.9, "PC 2", color = 'red', size = 15,fontweight="bold")
ax[1].tick_params(labelsize=15)
ax[1].set_xlabel("Temporal Position $t$", size = 15)
ax[1].set_xticks([1, 1.5, 2, 2.5, 3])
ax[1].set_xticklabels([1, 1.5, 2, 2.5, 3])
ax[1].set_title("Temporal Modes: PCs",fontweight="bold")
ax[1].axes.get_yaxis().set_visible(False)
# axes.get_yaxis().set_visible lets us change the y-axis' visibility.
# We want it gone here, because we're using the same y ticks from
# ax[0] and there wouldn't be room between the plots anyways. An
# x-axis counterpart also exists: axes.get_xaxis().set_visible

fig.tight_layout(pad=-1.5)
#plt.savefig("EOFs&PCs.jpg")


#Generate random data on a 10-by-15 grid with 20 time points 
#SVD for 2D spatial dimension and 1D time
np.random.seed(100)
dat = np.random.standard_normal(size=(10*15, 20))
# np.random.standard_normal creates an array of random values in the
# standard_normal distribution, size specifies the array's dimensions.

x = np.arange(10) + 1
y = np.arange(15) + 1
# Adding 1 to each of these numpy ranges increments each element
# by 1. Remember vectorization!

U, D, V = np.linalg.svd(dat, full_matrices=False)

np.shape(U)

(150, 20)


np.shape(V)

(20, 20)


len(D)

20


umat = U[:,0].reshape(15, 10, order='F')
# .reshape lets you change the numpy array's dimensions to the
# dimensions specified. Order is the way in which elements are
# taken from the original and added to the newly shaped array.
# 'F' is for Fortran-like (column major).
np.shape(umat)

(15, 10)


clev = np.linspace(umat.min(), umat.max(), 25)
# The number of levels we want for the contourf below.


#Figure 5.3
fig = plt.figure(figsize=(12,12))
# This simply creates the figure, the other plot stuff is
# below.

ax1 = plt.subplot2grid((2,1), (0,0))
# This implicitly creates a 2 by 1 array of plots in fig, and
# ax1 refers to the first. This way of declaring the plot has
# to be used instead of ways already articulated, because the
# colorbar on the contourf plot is incompatible with the old
# way. This method, however, is compatible.

def one_decimal_format(x, pos):
    """Format tick labels using one decimal place."""
    return f'{x:.1f}'

contf = ax1.contourf(x, y, umat, clev, cmap=plt.cm.get_cmap('bwr').reversed())
colbar = plt.colorbar(contf, drawedges=True, aspect=8)
colbar.formatter = FuncFormatter(one_decimal_format)
colbar.update_ticks()
#colbar.set_ticks([-.2, -.1, 0, .1, .2])
colbar.ax.tick_params(labelsize=15)
# This sets the ticks that show up on the colorbar manually,
# much like setting the ticks on the axes manually.

ax1.set_xlabel("Spatial Position: $x$", size = 15)
ax1.set_ylabel("Spatial Position: $y$", size = 15)
ax1.set_title("Spatial Pattern: EOF 1", size=20,fontweight="bold")
ax1.tick_params(labelsize=15)
ax1.grid()

ax2 = plt.subplot2grid((2,1), (1,0))
# This creates the second plot in our 2 by 1 array of plots.
# ax2 is the second.

ax2.plot(np.arange(V.shape[0]) + 1, V[:,0], 'ro-')
ax2.set_xlabel("Time", size = 15)
ax2.set_ylabel("PC1 Values [Dimensionless]", size = 15)
ax2.tick_params(labelsize=15)
ax2.set_title("Temporal Pattern: Time series", size=20,fontweight="bold")
ax2.grid()

fig.tight_layout(pad=1)

#fig.savefig("Spatial&TemporalPatterns.jpg", bbox_inches='tight')#


os.chdir('/Users/momtaza/Desktop/RMathModel/data')
months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
      'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
pta = pd.read_csv("PSTANDtahiti.txt", sep="\s+", index_col=0, names=months)
# Here we use pd.read_csv to read a text file, using column 0 as the
# label for each row, and using the months list as labels for the columns
ptamonv = np.array(pta).reshape(780, 1)


# Different txt file
pda = pd.read_csv("PSTANDdarwin.txt", sep="\s+", index_col=0, names=months)
pdamonv = np.array(pda).reshape(780, 1)


ptada = np.append(ptamonv,pdamonv, axis=1).T
# np.append combines lists/arrays into one by adding the elements
# of the second to the end of the first. Specifying the axis, lets
# you determine how they'll be added to each other. In this case,
# appending one array of shape (1, 780) to another array of the
# same shape along axis 0 yields an array of shape (2, 780). Note
# that appending it like that requires that each arrays dimension
# 2 must be equal.

print(f"Shape of ptada: {ptada.shape}")

u2, s2, v2 = np.linalg.svd(ptada, full_matrices=False)

print(f"Shape of Left-Singular u: {u2.shape},\n"
      f"Shape of Singular Values array s: {s2.shape},\n"
      f"Shape of Right-Singular v: {v2.shape}")

Shape of ptada: (2, 780)
Shape of Left-Singular u: (2, 2),
Shape of Singular Values array s: (2,),
Shape of Right-Singular v: (2, 780)


xtime = np.arange(1951, 2016 - 1/12, 1/12)
# xtime represents all months from January 1951 to December 2015

wsoi1 = s2[0]*v2[0]
wsoi2 = s2[1]*v2[1]
# Vectorized product of singular values with corresponding PCs.

cwsoi1 = np.cumsum(wsoi1)
cwsoi2 = np.cumsum(wsoi2)
# np.cumsum does the cumulative sum of each element in the array,
# and returns an array of the same length where element i is the
# sum of all elements up to i in the original array.

cnegsoi = np.cumsum(pdamonv-ptamonv)


#Figure 5.4
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(15, 12))
# We're not limited to n by 1 dimensional arrays of plots.
# This figure has 4 plots in a 2 by 2 arrangement.

ax[0,0].plot(xtime, ptamonv, 'r-')
ax[0,0].set_xlabel("Year", size=10)
ax[0,0].set_ylabel("Pressure", size=10)
ax[0,0].tick_params(labelsize=10)
ax[0,0].set_title("Standarized Tahiti SLP Anomalies",fontweight="bold",size=15)
ax[0,0].set_ylim(wsoi1.min(), wsoi1.max())
ax[0,0].grid()

ax[0,1].plot(xtime, pdamonv, 'b-')
ax[0,1].set_xlabel("Year", size=10)
ax[0,1].set_ylabel("Pressure", size=10)
ax[0,1].tick_params(labelsize=10)
ax[0,1].set_title("Standarized Darwin SLP Anomalies",fontweight="bold",size=15)
ax[0,1].set_ylim(wsoi2.min(), wsoi2.max())
ax[0,1].grid()

ax[1,0].plot(xtime, ptamonv - pdamonv, 'k-')
ax[1,0].set_xlabel("Year", size=10)
ax[1,0].set_ylabel("Cumulative Weighted SOI 2", size=10)
ax[1,0].tick_params(labelsize=10)
ax[1,0].grid()

ax[1,1].plot(xtime, cwsoi1, 'r-', label="Cumulative WSOI 1")
ax[1,1].plot(xtime, cnegsoi, 'k-', label="Cumulative Negative SOI")
ax[1,1].tick_params(labelsize=10)
ax[1,1].set_xlabel("Year", size=10)
ax[1,1].set_ylabel("Cumulative SOI", size=10)
ax[1,1].legend(loc='best')
ax[1,1].grid()

#plt.savefig("DarwinTahitiSOI.jpg")


np.shape(ptada)

(2, 780)


svdptd = np.linalg.svd(ptada, full_matrices = False)
svdptd_u = svdptd[0]
svdptd_d = np.diag([svdptd[1][0], svdptd[1][1]])
svdptd_v = svdptd[2]

# Verify that recontd = ptada
# print(svdptd_u, svdptd_d, svdptd_v)
recontd = svdptd_u @ svdptd_d @ svdptd_v
np.allclose(ptada, recontd) # True if matrices are (element-wise) equal

True


#Figure 5.5
fig, ax = plt.subplots(nrows=2,figsize=(18,9))

ax[0] = plt.subplot(2, 1, 1,
                    projection=cartopy.crs.PlateCarree(central_longitude=180))
# projection lets you change the nature of the plot. In this
# case, the plot is turned into a cartopy projection of the
# world. cartopy is a python package that allows for
# cartographic plotting.

ax[0].set_extent([0,359, -70, 70])
# set_extent lets you set the 'window' of the projection.
# The first 2 numbers are the min and max longitude, and
# the second 2 numbers are the min and max latitude.

ax[0].coastlines()
# coastlines draws the coasts of the landmasses within the window.

ax[0].plot(131-180, -12, 'bo', ms=26)
ax[0].text(110-180, -24, f"Darwin {-u2[1,0]:3.2f}", color='b', fontsize = 13)
ax[0].plot(211-180, -18, 'ro', ms=20)
ax[0].text(190-180, -30, f"Tahiti {-u2[0,0]:3.2f}", color='r', fontsize = 13)
ax[0].text(-105, -55, "El Niño Southern Oscillation Mode 1",
           color='purple', size=20,fontweight="bold")
# You can add as much floating text to a plot as you want.
# When setting the x and y ticks for a projection, you must include
# the crs type as well.
ax[0].set_xlabel("Longitude")
ax[0].set_ylabel("Latitude")
gl1 = ax[0].gridlines(crs=cartopy.crs.PlateCarree(), draw_labels=True,
    linewidth=1, color='gray', alpha=0.1, linestyle='-')

gl1.xlabels_top = False
gl1.ylabels_right = False
gl1.xlocator = mticker.FixedLocator([-225, -180, -135, -90, -45, 0, 45, 90, 135])
gl1.xformatter = LONGITUDE_FORMATTER
gl1.yformatter = LATITUDE_FORMATTER

ax[1] = plt.subplot(2, 1, 2,
                    projection=cartopy.crs.PlateCarree(central_longitude=180))
ax[1].set_extent([0,359, -70, 70])
ax[1].coastlines()
ax[1].plot(131-180, -12, 'ro', ms=20)
ax[1].text(110-180, -24, f"Darwin {u2[1,1]:3.2f}", color='r', fontsize = 13)
ax[1].plot(211-180,-18, 'ro', ms=26)
ax[1].text(190-180, -30, f"Tahiti {u2[1,0]:3.2f}", color='r', fontsize = 13)
ax[1].text(-105, -55, "El Niño Southern Oscillation Mode 2", 
           color='purple', size=20,fontweight="bold")

gl2 = ax[1].gridlines(crs=cartopy.crs.PlateCarree(), draw_labels=True,
                  linewidth=1, color='gray', alpha=0.1, linestyle='-')
gl2.xlabels_top = False
gl2.ylabels_right = False
gl2.xlocator = mticker.FixedLocator([-225, -180, -135, -90, -45, 0, 45, 90, 135])
gl2.xformatter = LONGITUDE_FORMATTER
gl2.yformatter = LATITUDE_FORMATTER
fig.tight_layout(pad=-1)

#plt.savefig("OscillationModesDarwinTahiti.jpg")

Introduction to Modern Mathematical Modeling with R:

A User's Manual to Train Mathematical Consultants

A Cambridge University Press book

by SSP Shen

Chapter 5: Mathematical Modeling By Linear Algebra¶

Kirchhoff’s Laws and Solution of an Electric Circuit¶

Leontief Production Model: A Balance of Output and Input¶

The Fundamental Idea of SVD: Space-Time-Energy Separation¶

SVD demo for generated data: A = UDV’¶

Example 5.2¶

SVD Analysis for El Niño Southern Oscillation Data¶

Figure 5.5¶