<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 1. Calculate an average spectrum to ID peaks
Step1:  2. Make a feature matrix, n x p, where n = number of samples, p = number of features
Step2:  3. Standardize
Step3:  4. Sklearn PCA
Step4:  5. Matrix decomposition (see A User's Guide to Principal Components by Jackson, 1991.)
U’SU  = L 
U = orthonormal matrix, characteristic vectors
S = covariance matrix
L = diagonal matrix, characteristic roots  
Get characteristic roots
Step5:  6. Matlab's PCA
Step6:  Check that all three give similar results
Step7:  PCA with the entire IR spectrum.
Step8:  Standardize matrix
<ASSISTANT_TASK:>
Python Code:
averagespectrum = PCAsynthetic.get_hyper_peaks(spectralmatrix, threshold = 0.01)
plt.plot(averagespectrum)
Explanation: 1. Calculate an average spectrum to ID peaks
End of explanation
featurematrix = PCAsynthetic.makefeaturematrix(spectralmatrix, averagespectrum)
featurematrix[10:13,:]
Explanation: 2. Make a feature matrix, n x p, where n = number of samples, p = number of features
End of explanation
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)
Explanation: 3. Standardize: zero mean, unit variance. (and check!)
End of explanation
#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_
Skvarianceratio
Skvariance
Explanation: 4. Sklearn PCA
End of explanation
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)
#check that the loadings squared sum to 1: 
Lsquared = sum(eig_vecs**2)
Explanation: 5. Matrix decomposition (see A User's Guide to Principal Components by Jackson, 1991.)
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}$]
End of explanation
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
Explanation: 6. Matlab's PCA
End of explanation
#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)
print(mlloadings[0,7])
print(eig_vecs[0,7])
Explanation: Check that all three give similar results
End of explanation
#Reshape spectral matrix 
IRmatrix=spectralmatrix.reshape(65536,559)
print(IRmatrix[1,:].shape)
#make sure we've reshaped correctly 
plt.plot(reshapespect[555,:])
Explanation: PCA with the entire IR spectrum.
End of explanation
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)
IRmatrix_std = PCAsynthetic.stdfeature(IRmatrix, axis = 0)
IRmean = IRmatrix_std.mean(axis=0)
IRvariance = IRmatrix_std.std(axis=0)
print(IRvariance)
IRmlPCA = PCA(IRmatrix_std)
#get projections of samples into PCA space
IRmltrans = IRmlPCA.Y
#reshape
IRmlloadings = IRmlPCA.Wt
IRmltrans.shape
IRmltransreshape=IRmltrans.reshape(256,256,360)
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]
plt.imshow(syntheticspectra.Cmatrix)
plt.imshow(score1image)
plt.imshow(score2image)
plt.imshow(score3image)
plt.imshow(score4image)
plt.imshow(score5image)
plt.imshow(score6image)
plt.imshow(score7image)
plt.imshow(score8image)
plt.imshow(score9image)
Explanation: Standardize matrix
End of explanation
<END_TASK>