#
# EPage 205
# Exercise 4.17 EPage 165
# Data from Appendix D, Table 10 EPage 722
#
library(car)

save_figs = FALSE

source('../../Data/data_loaders.R')
C = load_appendix_pitprop_data()

N = 180 # number of samples
m = 13 # number of measurements


#
# Part (a):
#

# Split the correlation matrix into its parts:
#
XTX = C[1:m, 1:m] # X^T X
XTy = C[1:m, m+1] # X^T y

# The VIFs:
#
VIFs = diag(solve(XTX))
print(VIFs)

# The eigenvalues:
#
ES = eigen(XTX)
print('Eigenvalues')
print(ES$values)


#
# Part (b):
#
ev_m = ES$vectors[, m]
ev_mm1 = ES$vectors[, m-1]

P = cbind( ev_m, ev_mm1 )
print( P )

# Print the invariants:
#
vn = colnames(C)[1:m]
sp = sprintf('%+.3f %s', ev_m, vn)
eq_m = paste(sp, sep='', collapse=' ')
print(sprintf('%s = 0', eq_m))

sp = sprintf('%+.3f %s', ev_mm1, vn)
eq_mm1 = paste(sp, sep='', collapse=' ')
print(sprintf('%s = 0', eq_mm1))

n_to_keep = 3

ev_order = order(abs(ev_m), decreasing=TRUE)
sp = sprintf('%+.3f %s', ev_m[ev_order][1:n_to_keep], vn[ev_order][1:n_to_keep])
eq_m = paste(sp, sep='', collapse=' ')
print(sprintf('%s = 0', eq_m))

ev_order = order(abs(ev_mm1), decreasing=TRUE)
sp = sprintf('%+.3f %s', ev_mm1[ev_order][1:n_to_keep], vn[ev_order][1:n_to_keep])
eq_m = paste(sp, sep='', collapse=' ')
print(sprintf('%s = 0', eq_m))


# Part (c):
#

# Consider possible linear combinations of the last two eigenvectors:
#
xs = seq( -0.5, +1.5, length.out=100 )
my_norm = c()
for( xi in 1:length(xs) ){
    my_norm = c( my_norm, mean(abs(xs[xi] * ev_m + (1-xs[xi]) * ev_mm1)) ) # I expect that at the correct linear combiation most coefficients are zero
}
plot( xs, my_norm, type='l' )
grid()

xi = which.min(my_norm) # find the index and x value of the smallest
print(sprintf('blend eigenvectors with %f', xs[xi]))

print('blended eigenvector looks like:')
print(round( xs[xi] * ev_m + (1-xs[xi]) * ev_mm1, 2 ))

print( c( vn[3], vn[4] ) ) # What are these variables: (MS = moist) = -(TG = testseg)