In [1]:
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import scale
from sklearn.decomposition import PCA as sklearnPCA
from matplotlib.mlab import PCA
import FindPeaks
import syntheticspectra
import PCAsynthetic
%matplotlib inline
In [3]:
spectralmatrix = syntheticspectra.spectralmatrix
In [14]:
averagespectrum = PCAsynthetic.get_hyper_peaks(spectralmatrix, threshold = 0.01)
plt.plot(averagespectrum)
Out[14]:
In [20]:
featurematrix = PCAsynthetic.makefeaturematrix(spectralmatrix, averagespectrum)
featurematrix[10:13,:]
Out[20]:
In [27]:
featurematrix_std = PCAsynthetic.stdfeature(featurematrix, axis = 0)
#along axis 0 = running vertically downwards, across rows; 1 = columns
mean = featurematrix_std.mean(axis=0)
variance = featurematrix_std.std(axis=0)
print(mean, variance)
In [61]:
#define number of principal components
sklearn_pca = sklearnPCA(n_components=9)
#matrix with each sample in terms of the PCs
SkPC = sklearn_pca.fit_transform(featurematrix_std)
#covariance matrix
Skcov = sklearn_pca.get_covariance()
#score matrix
#Skscore = sklearn_pca.score_samples(featurematrix_std)
#explained variance
Skvariance = sklearn_pca.explained_variance_
Skvarianceratio = sklearn_pca.explained_variance_ratio_
In [65]:
Skvarianceratio
Out[65]:
In [66]:
Skvariance
Out[66]:
U’SU = L
U = orthonormal matrix, characteristic vectors
S = covariance matrix
L = diagonal matrix, characteristic roots
Get characteristic roots: |𝑺−𝒍𝑰|=𝟎
Get characteristic vector: |𝑺−𝒍𝑰|𝒕𝒊=𝟎
The projection of sample n onto principal component i: z$_i$ = u$^{’}_{i}$[x$_{n}$-x$_{avg}$]
In [40]:
mean_vec = np.mean(featurematrix_std, axis=0)
#need to take transpose, since rowvar = true by default
cov_mat = np.cov(featurematrix_std.T)
#solve for characteristic roots and vectors
eig_vals, eig_vecs = np.linalg.eig(cov_mat)
In [67]:
#check that the loadings squared sum to 1:
Lsquared = sum(eig_vecs**2)
In [69]:
mlPCA = PCA(featurematrix_std)
#get projections of samples into PCA space
mltrans = mlPCA.Y
#reshape
mltransreshape = mltrans.reshape((256,256,9))
mlloadings = mlPCA.Wt
#mltrans[513,:] should be the same as mltransreshape[2,1,:]
mlloadings.shape
Out[69]:
In [72]:
#projection of first sample, on to the first PC
P11 = np.dot(eig_vecs[:,0], featurematrix_std[0,:]-mean_vec)
mlP11 = mlPCA.Y[0,0]
SkP11 = SkPC[0,0]
P12 = np.dot(eig_vecs[:,1], featurematrix_std[0,:]-mean_vec)
mlP12 = mlPCA.Y[0,1]
SkP12 = SkPC[0,1]
P152 = np.dot(eig_vecs[:,1], featurematrix_std[15,:]-mean_vec)
mlP152 = mlPCA.Y[15,1]
SkP152 = SkPC[15,1]
print(P11, mlP11, SkP11)
print(P12, mlP12, SkP12)
print(P152, mlP152, SkP152)
In [81]:
print(mlloadings[0,7])
print(eig_vecs[0,7])
In [111]:
#Reshape spectral matrix
IRmatrix=spectralmatrix.reshape(65536,559)
print(IRmatrix[1,:].shape)
In [112]:
#make sure we've reshaped correctly
plt.plot(reshapespect[555,:])
Out[112]:
In [115]:
IRmatrix=np.concatenate((IRmatrix[:,20:60], IRmatrix[:,230:270], IRmatrix[:,420:460], IRmatrix[:,100:140],IRmatrix[:,305:345], IRmatrix[:,470:510], IRmatrix[:,158:198], IRmatrix[:,354:394], IRmatrix[:,512:552] ), axis=1)
#IRmatrix=np.concatenate((IRmatrix[:,30:40], IRmatrix[:,240:260], IRmatrix[:,430:450], IRmatrix[:,90:130],IRmatrix[:,395:335], IRmatrix[:,460:500], IRmatrix[:,148:188], IRmatrix[:,364:384], IRmatrix[:,522:542] ), axis=1)
In [116]:
IRmatrix_std = PCAsynthetic.stdfeature(IRmatrix, axis = 0)
IRmean = IRmatrix_std.mean(axis=0)
IRvariance = IRmatrix_std.std(axis=0)
print(IRvariance)
In [117]:
IRmlPCA = PCA(IRmatrix_std)
#get projections of samples into PCA space
IRmltrans = IRmlPCA.Y
#reshape
IRmlloadings = IRmlPCA.Wt
IRmltrans.shape
Out[117]:
In [118]:
IRmltransreshape=IRmltrans.reshape(256,256,360)
In [126]:
score1image = IRmltransreshape[:,:,0]
score2image = IRmltransreshape[:,:,1]
score3image = IRmltransreshape[:,:,2]
score4image = IRmltransreshape[:,:,3]
score5image = IRmltransreshape[:,:,4]
score6image = IRmltransreshape[:,:,5]
score7image = IRmltransreshape[:,:,6]
score8image = IRmltransreshape[:,:,7]
score9image = IRmltransreshape[:,:,8]
In [4]:
plt.imshow(syntheticspectra.Cmatrix)
Out[4]:
In [120]:
plt.imshow(score1image)
Out[120]:
In [121]:
plt.imshow(score2image)
Out[121]:
In [122]:
plt.imshow(score3image)
Out[122]:
In [124]:
plt.imshow(score4image)
Out[124]:
In [125]:
plt.imshow(score5image)
Out[125]:
In [127]:
plt.imshow(score6image)
Out[127]:
In [128]:
plt.imshow(score7image)
Out[128]:
In [129]:
plt.imshow(score8image)
Out[129]:
In [130]:
plt.imshow(score9image)
Out[130]:
In [ ]: