Chapter 1: Basic Matrix Operations

Create your working directory named LinAlg and go there

Run the following, replacing “alyss/OneDrive/Documents” with your directory, to set the working directory: setwd(“/Users/alyss/OneDrive/Documents/LinAlg”).

Run getwd() to verify that you are in the directory/folder.

To save a figure

Before plotting the figure, run “setEPS(),” then run “postscript(”figname.eps”, width = w, height = h),” replacing figname, w, and h with the desired parameters. The figure will not show up as an output, and will instead be a file on your computer.

Some toy examples of R

1 + 2
## [1] 3
1 - 3
## [1] -2
x = c(1, 2)
y = c(3, 0)
plot(x, y)

plot(x,y, xlim = c(-3, 3),
     ylim = c(-1, 4),
     col = 'blue',
     pch = 16, type = 'o',
     lty = 2)

plot(1, col = 'red' , cex = 2, pch = 3)

x = c(1, 2)
y = c(3, 2)
plot(x, y, xlim = c(-2, 2), ylim = c(-3,3),
     col = 'blue', type = 'o')

plot(sin, c(-9,9), lwd =3, col = 'purple',
     main = 'Sam\'s color plot')

A = matrix(c(1,1,1,-1), nrow=2)
A
##      [,1] [,2]
## [1,]    1    1
## [2,]    1   -1
A1 = matrix(c(1, 1, 1, -1), 
            byrow = TRUE, 
            ncol = 2, nrow = 2)
A1
##      [,1] [,2]
## [1,]    1    1
## [2,]    1   -1
A2 = matrix(1, nrow =5, ncol =4)
A2
##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    1    1    1    1
## [3,]    1    1    1    1
## [4,]    1    1    1    1
## [5,]    1    1    1    1
B = matrix(c(1,3,2,0), nrow=2)
B
##      [,1] [,2]
## [1,]    1    2
## [2,]    3    0

Matrix/Vector Operations

#Matrix subtraction
A = matrix(c(1,0,3,-1), nrow=2) 
A
##      [,1] [,2]
## [1,]    1    3
## [2,]    0   -1
#R arranges matrix entries according to column
B = matrix(c(1,3,2,4), nrow=2)
B
##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4
A - B
##      [,1] [,2]
## [1,]    0    1
## [2,]   -3   -5
#Dot product
#install.packages('geometry')
library(geometry)
a = c(1, 1)
b = c(1, 0)
# Calculating dot product using dot()
dot(a, b)
## [1] 1
#Another way to compute dot product is to 
a%*%b
##      [,1]
## [1,]    1
#calculate the angle between two vectors
#Euclidean norm or called amplitude of a vector 
am = norm(a, type="2")  
am
## [1] 1.414214
bm = norm(b, type="2") 
bm
## [1] 1
#Angle in deg between two vectors
angleab = acos(dot(a, b)/(am*bm))*180/pi
angleab
## [1] 45
#Cross product of two 3D vectors
#install.packages('pracma')
library(pracma)
## Warning: package 'pracma' was built under R version 4.4.3
## 
## Attaching package: 'pracma'
## The following objects are masked from 'package:geometry':
## 
##     cart2pol, cart2sph, dot, pol2cart, polyarea, sph2cart
x = 1:3
y = 4:6
cross(x, y) 
## [1] -3  6 -3
dot(x, y)
## [1] 32
#Scalar times a matrix
A = matrix(c(1,0,3,-1), nrow=2) 
A
##      [,1] [,2]
## [1,]    1    3
## [2,]    0   -1
3*A
##      [,1] [,2]
## [1,]    3    9
## [2,]    0   -3
a
## [1] 1 1
A
##      [,1] [,2]
## [1,]    1    3
## [2,]    0   -1
b*A
##      [,1] [,2]
## [1,]    1    3
## [2,]    0    0
b
## [1] 1 0
#Matrix multiplication by command %*%
A = matrix(c(1,0,3,-1), nrow=2) 
B = matrix(c(1,3,2,4), nrow=2)
A%*%B
##      [,1] [,2]
## [1,]   10   14
## [2,]   -3   -4
B%*%A
##      [,1] [,2]
## [1,]    1    1
## [2,]    3    5
#Matrix transpose
C = matrix(c(1,2,3,4), ncol =2, byrow=T)
C
##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4
t(C)
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4
#Generate a diagonal matrix
D = diag(c(2, 1, -3))
round(D, digits = 0)
##      [,1] [,2] [,3]
## [1,]    2    0    0
## [2,]    0    1    0
## [3,]    0    0   -3
#Generate a 3-dimensional identity matrix
I = diag(3)
I
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
#Generate a 2-by-3 zero matrix
M = matrix(0, 2, 3)
M
##      [,1] [,2] [,3]
## [1,]    0    0    0
## [2,]    0    0    0
#Compute the inverse of a matrix
A = matrix(c(1,1,1,-1), nrow=2)
invA = solve(A)
invA
##      [,1] [,2]
## [1,]  0.5  0.5
## [2,]  0.5 -0.5
#Verify the inverse
invA %*% A
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

Examples: 9/9/2024 class

par(mfrow = c(1, 2))
H = matrix(rnorm(64), nrow = 8)
image(H)
G = solve(H)
image(G)

#dev.off()

filled.contour(H)
grid()

m = 9
n = 8
A = matrix(2, m, n)
A
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
##  [1,]    2    2    2    2    2    2    2    2
##  [2,]    2    2    2    2    2    2    2    2
##  [3,]    2    2    2    2    2    2    2    2
##  [4,]    2    2    2    2    2    2    2    2
##  [5,]    2    2    2    2    2    2    2    2
##  [6,]    2    2    2    2    2    2    2    2
##  [7,]    2    2    2    2    2    2    2    2
##  [8,]    2    2    2    2    2    2    2    2
##  [9,]    2    2    2    2    2    2    2    2
A = matrix(2, nrow = 9, ncol = 8)
A #print A
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
##  [1,]    2    2    2    2    2    2    2    2
##  [2,]    2    2    2    2    2    2    2    2
##  [3,]    2    2    2    2    2    2    2    2
##  [4,]    2    2    2    2    2    2    2    2
##  [5,]    2    2    2    2    2    2    2    2
##  [6,]    2    2    2    2    2    2    2    2
##  [7,]    2    2    2    2    2    2    2    2
##  [8,]    2    2    2    2    2    2    2    2
##  [9,]    2    2    2    2    2    2    2    2
A = matrix(1:72, nrow = 9, ncol = 8)
A #print A
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
##  [1,]    1   10   19   28   37   46   55   64
##  [2,]    2   11   20   29   38   47   56   65
##  [3,]    3   12   21   30   39   48   57   66
##  [4,]    4   13   22   31   40   49   58   67
##  [5,]    5   14   23   32   41   50   59   68
##  [6,]    6   15   24   33   42   51   60   69
##  [7,]    7   16   25   34   43   52   61   70
##  [8,]    8   17   26   35   44   53   62   71
##  [9,]    9   18   27   36   45   54   63   72
#visualize a matrix
image(A)

rnorm(40)
##  [1]  0.15488363  0.49858426  0.12387345 -0.60166366  0.48270959  1.27186149
##  [7] -1.32102699  0.25126418  1.22669472  0.47526113  1.23224267  1.50728211
## [13] -1.99499852 -0.26266593 -0.76228718  1.78290625  0.45740153  0.20637489
## [19]  1.24944627  1.31834462 -0.56074856 -1.47073660 -1.40757971 -0.07549309
## [25]  0.76095421  0.82699793  0.89565507 -0.28663244 -0.87078693  0.46445922
## [31]  0.90468661  0.49825656 -0.05694803 -0.48386256 -0.87560556 -0.24542624
## [37]  1.52198005 -0.29141272 -0.04053436 -1.28781869
hist(rnorm(4000), breaks = 99)

A = matrix(rnorm(800*1000), nrow = 800, ncol = 1000)
image(A)

Editing a matrix

#Delete rows or/and columns
A = matrix(1:9, nrow = 3)
A
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9
image(A)

A[-3,]
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
A[,-1]
##      [,1] [,2]
## [1,]    4    7
## [2,]    5    8
## [3,]    6    9
A[-3,-1]
##      [,1] [,2]
## [1,]    4    7
## [2,]    5    8
A[1:2, 2:3] #Sub-matrix
##      [,1] [,2]
## [1,]    4    7
## [2,]    5    8
A[1:2, ] #keep all the columns
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
#Insert a row or column to a matrix
A = matrix(1:4, nrow = 2)
br = 5:6
bc = 7:8
rbind(A, br)
##    [,1] [,2]
##       1    3
##       2    4
## br    5    6
rbind(br, A)
##    [,1] [,2]
## br    5    6
##       1    3
##       2    4
rbind(rbind(A[1,], br), A[2,])
##    [,1] [,2]
##       1    3
## br    5    6
##       2    4
Abc = cbind(A, bc)
Abc
##          bc
## [1,] 1 3  7
## [2,] 2 4  8
Abc = cbind(A, bc)

cbind(Abc, A)#stack two matrices
##          bc    
## [1,] 1 3  7 1 3
## [2,] 2 4  8 2 4

Row or column statistics

library(matrixStats)
A = matrix(1:6, nrow = 2)
A
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6
rowSums(A)
## [1]  9 12
rowMeans(A)
## [1] 3 4
library(matrixStats)
rowCumsums(A) #cumulative sum
##      [,1] [,2] [,3]
## [1,]    1    4    9
## [2,]    2    6   12
colMeans(A)
## [1] 1.5 3.5 5.5
#install.packages('matrixStats')
library(matrixStats) #SD needs the library
rowSds(A) 
## [1] 2 2
colSds(A)
## [1] 0.7071068 0.7071068 0.7071068
#Sweep a matrix by a vector using subtraction
A = matrix(1:6, nrow=2, byrow = TRUE)
A
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
u = 1:3
u
## [1] 1 2 3
Br = sweep(A, 2, u) #2 means sweep every row
Br
##      [,1] [,2] [,3]
## [1,]    0    0    0
## [2,]    3    3    3
v= 1:2
Bc = sweep(A, 1, v) #1 means sweep every column
Bc
##      [,1] [,2] [,3]
## [1,]    0    1    2
## [2,]    2    3    4
c = colMeans(A) #means of each column
sweep(A, 2, c) #anomaly data matrix
##      [,1] [,2] [,3]
## [1,] -1.5 -1.5 -1.5
## [2,]  1.5  1.5  1.5
A
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
sin(A)#function operation on each matrix element
##            [,1]       [,2]       [,3]
## [1,]  0.8414710  0.9092974  0.1411200
## [2,] -0.7568025 -0.9589243 -0.2794155
A^2 #not equal to A%*%A
##      [,1] [,2] [,3]
## [1,]    1    4    9
## [2,]   16   25   36
#Sweep a matrix by a vector using multiplication
A = matrix(1:6, nrow=2, byrow = TRUE)
w = 1:2
A
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
w
## [1] 1 2
w*A
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    8   10   12
A*w # yields the same result as w*A
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    8   10   12
w3 = 1:3
#sweep each row by transposing A
t(w3*t(A))
##      [,1] [,2] [,3]
## [1,]    1    4    9
## [2,]    4   10   18
w3*A #multiplication sweep by row-dimensions not matching
##      [,1] [,2] [,3]
## [1,]    1    6    6
## [2,]    8    5   18
A/w #sweeping by division
##      [,1] [,2] [,3]
## [1,]    1  2.0    3
## [2,]    2  2.5    3

Conversions between a Vector and a Matrix

#Conversions between a Vector and a Matrix
v = c(60, 58, 67, 70, 55, 53)
M = matrix(v, nrow = 2) #from vector to matrix
M
##      [,1] [,2] [,3]
## [1,]   60   67   55
## [2,]   58   70   53
c(M) #from matrix to vector by column
## [1] 60 58 67 70 55 53
c(t(M)) #from matrix to vector by row
## [1] 60 67 55 58 70 53
#Reduce the dimension of an nD array
x <- array(1:(2*3*4), dim=c(2,3,4))
dim(x)
## [1] 2 3 4
x #a stack of four 2-by-3 matrices
## , , 1
## 
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6
## 
## , , 2
## 
##      [,1] [,2] [,3]
## [1,]    7    9   11
## [2,]    8   10   12
## 
## , , 3
## 
##      [,1] [,2] [,3]
## [1,]   13   15   17
## [2,]   14   16   18
## 
## , , 4
## 
##      [,1] [,2] [,3]
## [1,]   19   21   23
## [2,]   20   22   24
#install.packages('R.utils')
library(R.utils)
## Warning: package 'R.utils' was built under R version 4.4.2
## Loading required package: R.oo
## Loading required package: R.methodsS3
## R.methodsS3 v1.8.2 (2022-06-13 22:00:14 UTC) successfully loaded. See ?R.methodsS3 for help.
## R.oo v1.27.0 (2024-11-01 18:00:02 UTC) successfully loaded. See ?R.oo for help.
## 
## Attaching package: 'R.oo'
## The following object is masked from 'package:R.methodsS3':
## 
##     throw
## The following objects are masked from 'package:methods':
## 
##     getClasses, getMethods
## The following objects are masked from 'package:base':
## 
##     attach, detach, load, save
## R.utils v2.12.3 (2023-11-18 01:00:02 UTC) successfully loaded. See ?R.utils for help.
## 
## Attaching package: 'R.utils'
## The following object is masked from 'package:utils':
## 
##     timestamp
## The following objects are masked from 'package:base':
## 
##     cat, commandArgs, getOption, isOpen, nullfile, parse, use, warnings
#flat all the other dim except the 3rd one
#flat the 1st and 2nd dim
y <- wrap(x, map=list(3, NA)) 
dim(y)
## [1] 4 6
y
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    2    3    4    5    6
## [2,]    7    8    9   10   11   12
## [3,]   13   14   15   16   17   18
## [4,]   19   20   21   22   23   24
#back to the original 3D array
array(t(y), dim = c(2,3,4))
## , , 1
## 
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6
## 
## , , 2
## 
##      [,1] [,2] [,3]
## [1,]    7    9   11
## [2,]    8   10   12
## 
## , , 3
## 
##      [,1] [,2] [,3]
## [1,]   13   15   17
## [2,]   14   16   18
## 
## , , 4
## 
##      [,1] [,2] [,3]
## [1,]   19   21   23
## [2,]   20   22   24

Solve linear equations

A = matrix(c(1,1,1,-1), nrow = 2)
b = c(20,4)
solve(A, b) #This is the result x1=12, and x2=8.
## [1] 12  8
A = matrix(c(1,2,-1,4),nrow =2)
A
##      [,1] [,2]
## [1,]    1   -1
## [2,]    2    4
res = eigen(A)
res$values
## [1] 3 2
res$vectors
##            [,1]       [,2]
## [1,]  0.4472136 -0.7071068
## [2,] -0.8944272  0.7071068
B = matrix(rnorm(6), nrow = 2)
B
##             [,1]       [,2]       [,3]
## [1,] -0.93849351 -1.0888761  1.6988421
## [2,]  0.05096902  0.3549143 -0.9520125
svd(B)
## $d
## [1] 2.4145878 0.3964105
## 
## $u
##            [,1]      [,2]
## [1,] -0.9193904 0.3933463
## [2,]  0.3933463 0.9193904
## 
## $v
##            [,1]       [,2]
## [1,]  0.3656485 -0.8130274
## [2,]  0.4724228 -0.2573105
## [3,] -0.8019463 -0.5222814

Spatial covariance matrix and SVD

dat = matrix(c(0,-1,1,2,3,4), nrow=3)
dat
##      [,1] [,2]
## [1,]    0    2
## [2,]   -1    3
## [3,]    1    4
colMeans(dat)
## [1] 0 3
A = sweep(dat, 2, colMeans(dat))
A
##      [,1] [,2]
## [1,]    0   -1
## [2,]   -1    0
## [3,]    1    1
covm=(1/(dim(A)[2]))*A%*%t(A)
covm #is the covariance matrix.
##      [,1] [,2] [,3]
## [1,]  0.5  0.0 -0.5
## [2,]  0.0  0.5 -0.5
## [3,] -0.5 -0.5  1.0
u = c(1, 1, 0)
v = covm %*% u
v #u and v are in different directions
##      [,1]
## [1,]  0.5
## [2,]  0.5
## [3,] -1.0
#Eigenvectors of a covariance matrix 
ew = eigen(covm)
ew
## eigen() decomposition
## $values
## [1] 1.500000e+00 5.000000e-01 1.332268e-15
## 
## $vectors
##            [,1]          [,2]      [,3]
## [1,] -0.4082483 -7.071068e-01 0.5773503
## [2,] -0.4082483  7.071068e-01 0.5773503
## [3,]  0.8164966  8.881784e-16 0.5773503
#Verify the eigenvectors and eigenvalues
covm%*%ew$vectors[,1]/ew$values[1]
##            [,1]
## [1,] -0.4082483
## [2,] -0.4082483
## [3,]  0.8164966
w = ew$vectors[,1] #is an eigenvector

R code for SVD

#Develop a 2-by-3 space-time data matrix for SVD
A=matrix(c(1,-1,2,0,3,1),nrow=2)
A
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]   -1    0    1
#Perform SVD calculation
msvd=svd(A)
msvd
## $d
## [1] 3.784779 1.294390
## 
## $u
##            [,1]       [,2]
## [1,] -0.9870875 -0.1601822
## [2,] -0.1601822  0.9870875
## 
## $v
##            [,1]       [,2]
## [1,] -0.2184817 -0.8863403
## [2,] -0.5216090 -0.2475023
## [3,] -0.8247362  0.3913356
msvd$d
## [1] 3.784779 1.294390
msvd$u
##            [,1]       [,2]
## [1,] -0.9870875 -0.1601822
## [2,] -0.1601822  0.9870875
msvd$v
##            [,1]       [,2]
## [1,] -0.2184817 -0.8863403
## [2,] -0.5216090 -0.2475023
## [3,] -0.8247362  0.3913356
#One can verify that A=UDV', where V' is transpose of V.
verim=msvd$u%*%diag(msvd$d)%*%t(msvd$v)
verim
##      [,1]         [,2] [,3]
## [1,]    1 2.000000e+00    3
## [2,]   -1 5.551115e-17    1
round(verim) #This is the original data matrix A
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]   -1    0    1
covm = (1/(dim(A)[2]))*A%*%t(A)
eigcov = eigen(covm)
eigcov$values
## [1] 4.7748518 0.5584816
eigcov$vectors
##            [,1]       [,2]
## [1,] -0.9870875  0.1601822
## [2,] -0.1601822 -0.9870875
((msvd$d)^2)/(dim(A)[2])
## [1] 4.7748518 0.5584816
eigcov$values
## [1] 4.7748518 0.5584816

Regression examples

x1=c(1,2,3) #Given the coordinates of the 3 points
x2=c(2,1,3)
y=c(-1,2,1)
df=data.frame(x1,x2,y) #Put data into the data.frame format
fit <- lm(y ~ x1 + x2, data=df)
fit #Show the regression results
## 
## Call:
## lm(formula = y ~ x1 + x2, data = df)
## 
## Coefficients:
## (Intercept)           x1           x2  
##  -5.128e-16    1.667e+00   -1.333e+00
1.667*x1-1.333*x2  #Verify that 3 points determining a plane
## [1] -0.999  2.001  1.002
#Multilinear regression
u=c(1,2,3,1)
v=c(2,4,3,-1)
w=c(1,-2,3,4)
mydata=data.frame(u,v,w)
myfit <- lm(w ~ u + v, data=mydata)
summary(myfit) #Show the result
## 
## Call:
## lm(formula = w ~ u + v, data = mydata)
## 
## Residuals:
##    1    2    3    4 
##  1.0 -1.0  0.5 -0.5 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   1.0000     1.8708   0.535    0.687
## u             2.0000     1.2472   1.604    0.355
## v            -1.5000     0.5528  -2.714    0.225
## 
## Residual standard error: 1.581 on 1 degrees of freedom
## Multiple R-squared:  0.881,  Adjusted R-squared:  0.6429 
## F-statistic:   3.7 on 2 and 1 DF,  p-value: 0.345
#Multilinear regression example for more data
dat=matrix(rnorm(40),nrow=10, dimnames=list(c(letters[1:10]), c(LETTERS[23:26])))
fdat=data.frame(dat)
fit=lm(Z~ W + X + Y, data=fdat)
summary(fit)
## 
## Call:
## lm(formula = Z ~ W + X + Y, data = fdat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4197 -0.3044 -0.1890  0.3377  0.6640 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.089407   0.176831   0.506   0.6312  
## W            0.081320   0.253697   0.321   0.7594  
## X           -0.419109   0.167364  -2.504   0.0463 *
## Y            0.002747   0.185969   0.015   0.9887  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5191 on 6 degrees of freedom
## Multiple R-squared:  0.5732, Adjusted R-squared:  0.3598 
## F-statistic: 2.686 on 3 and 6 DF,  p-value: 0.14

Image analysis with SVD

#Ref: imager: an R package for image processing
#https://dahtah.github.io/imager/imager.html 

setwd("/Users/alyss/OneDrive/Documents/LinAlg")
#install.packages('imager') #run this if not installed yet
library(imager) 
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:R.utils':
## 
##     extract
## The following object is masked from 'package:R.oo':
## 
##     equals
## The following objects are masked from 'package:pracma':
## 
##     and, mod, or
## 
## Attaching package: 'imager'
## The following object is masked from 'package:magrittr':
## 
##     add
## The following objects are masked from 'package:stats':
## 
##     convolve, spectrum
## The following object is masked from 'package:graphics':
## 
##     frame
## The following object is masked from 'package:base':
## 
##     save.image
dat <- load.image('data/SamPhoto.png') #355 KB file size
dim(dat)
## [1] 430 460   1   3
#[1] 430 460   1   3
#430 rows and 460 columns, 430*460 = 197,800 pixels
#1 photo frame, 3 RGB colors 
#If a video, 1 will become 150 frames or more

dat[1:3,1:4,1,1]
##           [,1]      [,2]      [,3]      [,4]
## [1,] 0.5176471 0.5058824 0.5215686 0.5490196
## [2,] 0.5215686 0.5137255 0.5215686 0.5294118
## [3,] 0.5215686 0.4901961 0.5058824 0.5215686
image(dat[, , 1, 3])

#Show part of the data

#plot the color figure
plot(dat, 
     xlim = c(0, 430), ylim = c(460, 0),
     main = 'A Color Photo of Sam') #switch 0 and 460 around to ensure the picture comes out right side up

#Make the photo black-and-white
#dat = dat[,,1,1:3]
graydat = grayscale(dat)
dim(graydat)
## [1] 430 460   1   1
#[1] 430 460   1   1
#430 rows and 460 columns,  1 photo frame, 1 grayscale [0, 1]
#plot the gray b/w photo
graydat[1:3, 1:5, 1, 1]
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.6075294 0.6126667 0.6215686 0.6293725 0.6298431
## [2,] 0.6072549 0.6112549 0.6186667 0.6178431 0.6271373
## [3,] 0.6053725 0.6098431 0.6107843 0.6144706 0.6187059
image(graydat[,,1,1])

graydat1 = graydat[,,1,1] - 
  0.1*matrix(rchisq(430*460, df =5), nrow = 430)
image(graydat1)

plot(graydat, 
     xlim = c(0, 430), ylim = c(460, 0),
     main = 'B/W Gray Sam')

#Plot the color and b/w photos together
par(mfrow = c(1, 2))
plot(dat, 
     xlim = c(0, 430), ylim = c(460, 0),
     main = 'Color Sam')
plot(graydat, 
     xlim = c(0, 430), ylim = c(460, 0),
     main = 'B/W Sam')

#dev.off()

#SVD analysis of the grayscale data
svdDat = svd(graydat)
SVDd = svdDat$d
percentD = 100*(SVDd^2)/sum(SVDd^2)
cumpercentD = cumsum(percentD)
modeK = 1:length(SVDd)
#dev.off()
plot(modeK[1:20], percentD[1:20], 
     type = 'o', col = 'blue',
     xlab = 'Mode number', pch = 16,
     ylab = 'Percentage of mode variance',
     main = 'Scree Plot of SVD B/W Photo Data')
K = 20
lam = (svdDat$d)^2
lamK=lam[1:K]
lamK
##  [1] 91599.58514  2169.70792   691.45565   454.54794   331.83214   218.25663
##  [7]   188.95894   161.15557   127.79007   124.10945    94.67419    78.18529
## [13]    61.25030    51.53846    44.45430    40.95168    38.09727    32.50086
## [19]    28.94096    27.76117
par(mar=c(4,4,2,4), mgp=c(2.2,0.7,0))
plot(1:K, 100*lamK/sum(lam), ylim=c(0,100), type="o", 
     ylab="Percentage of Variance [%]",
     xlab="EOF Mode Number", 
     cex.lab=1.2, cex.axis = 1.1, lwd=2, 
     main="Scree Plot of the First 20 Eigenvalues")
legend(3,30, col=c("black"),lty=1, lwd=2.0,
       legend=c("Percentange Variance"),bty="n",
       text.font=2,cex=1.0, text.col="black")
par(new=TRUE)
plot(1:K,cumsum(100*lamK/sum(lam)),
     ylim = c(90,100), type="o",
     col="blue",lwd=2, axes=FALSE,
     xlab="",ylab="")
legend(3,94.5, col=c("blue"),lty=1,lwd=2.0,
       legend=c("Cumulative Percentage Variance"),bty="n",
       text.font=2,cex=1.0, text.col="blue")
axis(4, col="blue", col.axis="blue", mgp=c(3,0.7,0))
mtext("Cumulative Variance [%]",col="blue", 
      cex=1.2, side=4,line=2)

#Reconstructing b/w photo from SVD modes
U = svdDat$u
V = svdDat$v
D = diag(svdDat$d)
dim(D)
## [1] 430 430
#100% reconstruction using all the modes
recon = U%*%D%*%t(V)
dim(recon)
## [1] 430 460
#[1] 430 460
image(recon)

#dev.off()
par(mfrow = c(2, 2))

kR = 430 #Recon from all 413 modes
R430 = as.cimg(U[,1:kR]%*%D[1:kR, 1:kR]%*%t(V[, 1:kR]))
plot(R430, main = "All 430 modes")
kB = 20 #The first 20 modes
R20 = as.cimg(U[,1:kB]%*%D[1:kB, 1:kB]%*%t(V[, 1:kB]))
plot(R20, main = "The first 20 modes")
k = 3 #Recon from the first 3 modes
R3 = as.cimg(U[,1:k]%*%D[1:k, 1:k]%*%t(V[, 1:k]))
plot(R3, main = "The first 3 modes")
k1 = 21; k2 = 100 #Recon from 21st to 100th modes
Rk1_k2 = as.cimg(U[,k1:k2]%*%D[k1:k2, k1:k2]%*%t(V[, k1:k2]))
plot(Rk1_k2, main = "21st to 100th modes")

#dev.off()

#Three monotone photos and their color photo
dim(dat) #4D array of a color photo
## [1] 430 460   1   3
#[1] 430 460   1   3
#dev.off()
par(mfrow = c(2, 2))
R = as.cimg(apply(t(dat[,, 1, 1]), 1, rev))
plot(R, main = 'Monotone color R')
G = as.cimg(apply(t(dat[,, 1, 2]), 1, rev))
plot(G, main = 'Monotone color G')
B = as.cimg(apply(t(dat[,, 1, 3]), 1, rev))
plot(B, , main = 'Monotone color B')
#Bind the three monotone photos into one 
trippy = imappend(list(R,G,B), "c") 
dim(trippy) #color figure data
## [1] 430 460   1   3
#[1] 430 460   1   3
plot(trippy, main ='Blend RGB Colors')

#dev.off()