# coding: utf-8
import numpy as np
import matplotlib.pyplot as plt
import sklearn.decomposition
import matplotlib.cm as cm
import pandas as pd

wl=[0.48,0.55,0.59,0.7,0.86,0.95]

X=np.loadtxt('smass2_c_input.txt')
pca=sklearn.decomposition.PCA(n_components=4)
pca.fit(X)
pc=pca.components_
mn=pca.mean_
X_pca=pca.transform(X)
z=X-mn
pc1score=np.dot(z,pc[0])

plt.figure(figsize=(8,8))
plt.plot(wl[:-1],pc.T[:,0],'o-')
plt.plot(wl[:-1],pc.T[:,1],'o-')
plt.plot(wl[:-1],-pc.T[:,2],'o-')
plt.legend(['PC1','PC2','PC3'],fontsize=16)
plt.xlabel('Wavelength(nm)',fontsize=16)
plt.ylabel('Coefficients', fontsize=16)
plt.title('Principal Components',fontsize=16)
plt.grid()
plt.show()

plt.figure(figsize=(8,2))
plt.plot(wl[:-1],mn,'ko-')
plt.xlabel('Wavelength(nm)')
plt.title('Mean spectrum')
plt.ylabel('Normalized ref.')
plt.show()

smass=pd.read_csv('SMASS_ONC_filter.csv')
smass2=smass[smass.smass=='smass2']
smass2_c=smass2[(smass2.tax_bus=='C')|(smass2.tax_bd=='C')|(smass2.tax_bus=='Cb')|(smass2.tax_bd=='Cb')|(smass2.tax_bus=='B')|(smass2.tax_bd=='B')|(smass2.tax_bus=='Ch')|(smass2.tax_bd=='Ch')
                |(smass2.tax_bus=='Cgh')|(smass2.tax_bd=='Cgh')|(smass2.tax_bus=='Cg')|(smass2.tax_bd=='Cg')]

def get_rotation_matrix(theta):
    rad=theta/180.*np.pi
    rot = np.array([[np.cos(rad), -np.sin(rad)],
                  [np.sin(rad), np.cos(rad)]])
    return rot


plt.figure(figsize=(10,10))
ax=plt.subplot()
ax.patch.set_alpha(0.)
plt.grid()

size=50
a=0.7
plt.scatter(X_pca[(smass2_c.tax_bus=='C')|(smass2_c.tax_bd=='C'),1],X_pca[(smass2_c.tax_bus=='C')|(smass2_c.tax_bd=='C'),2],c='r',label='C',s=size,alpha=a)
plt.scatter(X_pca[(smass2_c.tax_bus=='Ch')|(smass2_c.tax_bd=='Ch'),1],X_pca[(smass2_c.tax_bus=='Ch')|(smass2_c.tax_bd=='Ch'),2],c='gold',label='Ch',s=size,alpha=a)
plt.scatter(X_pca[(smass2_c.tax_bus=='Cg')|(smass2_c.tax_bd=='Cg'),1],X_pca[(smass2_c.tax_bus=='Cg')|(smass2_c.tax_bd=='Cg'),2],c='darkgreen',label='Cg',s=size,alpha=a)
plt.scatter(X_pca[(smass2_c.tax_bus=='Cgh')|(smass2_c.tax_bd=='Cgh'),1],X_pca[(smass2_c.tax_bus=='Cgh')|(smass2_c.tax_bd=='Cgh'),2],c='lime',label='Cgh',s=size,alpha=a)
plt.scatter(X_pca[(smass2_c.tax_bus=='B')|(smass2_c.tax_bd=='B'),1],X_pca[(smass2_c.tax_bus=='B')|(smass2_c.tax_bd=='B'),2],c='b',label='B',s=size,alpha=a)
plt.scatter(X_pca[(smass2_c.tax_bus=='Cb')|(smass2_c.tax_bd=='Cb'),1],X_pca[(smass2_c.tax_bus=='Cb')|(smass2_c.tax_bd=='Cb'),2],c='c',label='c',s=size,alpha=a)
plt.legend()
plt.title('PCA: b to x',fontsize=18)
plt.ylabel('C-complex PC3',fontsize=18)
plt.xlabel('C-complex PC2',fontsize=18)
x=[-1,1]
m = np.tan(35/180.*np.pi)
y=np.array(x)*(m)
xy=np.dot(get_rotation_matrix(90),[x,y])
plt.plot(x,y,ls='dashed')
plt.plot(xy[0],xy[1],label='histogram x')

plt.legend(loc=4, fontsize=14)
plt.ylim(0.06,-0.06)
plt.xlim(0.1,-0.1)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
#plt.axis('equal')
plt.show()


#//// rotate data
rot=get_rotation_matrix(90-35)
dataC=[X_pca[(smass2_c.tax_bus=='C')|(smass2_c.tax_bd=='C'),1],X_pca[(smass2_c.tax_bus=='C')|(smass2_c.tax_bd=='C'),2]]
dataCh=[X_pca[(smass2_c.tax_bus=='Ch')|(smass2_c.tax_bd=='Ch'),1],X_pca[(smass2_c.tax_bus=='Ch')|(smass2_c.tax_bd=='Ch'),2]]
dataCg=[X_pca[(smass2_c.tax_bus=='Cg')|(smass2_c.tax_bd=='Cg'),1],X_pca[(smass2_c.tax_bus=='Cg')|(smass2_c.tax_bd=='Cg'),2]]
dataCgh=[X_pca[(smass2_c.tax_bus=='Cgh')|(smass2_c.tax_bd=='Cgh'),1],X_pca[(smass2_c.tax_bus=='Cgh')|(smass2_c.tax_bd=='Cgh'),2]]
dataB=[X_pca[(smass2_c.tax_bus=='B')|(smass2_c.tax_bd=='B'),1],X_pca[(smass2_c.tax_bus=='B')|(smass2_c.tax_bd=='B'),2]]
dataCb=[X_pca[(smass2_c.tax_bus=='Cb')|(smass2_c.tax_bd=='Cb'),1],X_pca[(smass2_c.tax_bus=='Cb')|(smass2_c.tax_bd=='Cb'),2]]
rotxC=np.dot(rot,dataC)
rotxCh=np.dot(rot,dataCh)
rotxCg=np.dot(rot,dataCg)
rotxCgh=np.dot(rot,dataCgh)
rotxB=np.dot(rot,dataB)
rotxCb=np.dot(rot,dataCb)
rotxXY=np.dot(rot,[x,y])
rotxxy=np.dot(rot,xy)

data=rotxCb[0]
for i in range(0,len(data)):
    x=data[i]
    
histC, bins0=np.histogram(rotxC[0], bins=50,range=(-0.08,0.08))
histCh, bins=np.histogram(rotxCh[0], bins=50,range=(-0.08,0.08))
histCg, bins=np.histogram(rotxCg[0], bins=50,range=(-0.08,0.08))
histCgh, bins=np.histogram(rotxCgh[0], bins=50,range=(-0.08,0.08))
histB, bins=np.histogram(rotxB[0], bins=50,range=(-0.08,0.08))
histCb, bins=np.histogram(rotxCb[0], bins=50,range=(-0.08,0.08))
plt.figure(figsize=(8,6))
ax=plt.subplot()
label=['C','Ch','Cg','Cgh','B','Cb']
color=['r','yellow','darkgreen','lime','b','lightblue']
hatchlist = ['/','/',None,None,None,None]
plt.hist([rotxC[0],rotxCh[0],rotxCg[0],rotxCgh[0],rotxB[0],rotxCb[0]],stacked=True,bins=50,label=label,color=color,ec='black',hatch='/')

ax.patch.set_alpha(0.)
plt.title('histogram')
plt.ylabel('',fontsize=12)
plt.xlabel('P axis',fontsize=12)
plt.xlim(0.06,-0.06)
plt.legend()
plt.show()

ryu=np.loadtxt('Ryugu_global_BoxA_l2er_b2x.txt')
ryu_pc=pca.transform(ryu)

plt.figure(figsize=(10,10))
plt.grid()
plt.hist2d(ryu_pc[:,1],ryu_pc[:,2],bins=[np.linspace(-0.1,0.1,100),np.linspace(-0.1,0.1,100)],cmap=cm.gray_r)
plt.scatter(X_pca[(smass2_c.tax_bus=='C')|(smass2_c.tax_bd=='C'),1],X_pca[(smass2_c.tax_bus=='C')|(smass2_c.tax_bd=='C'),2],c='r',label='C',s=size,alpha=a)
plt.scatter(X_pca[(smass2_c.tax_bus=='Ch')|(smass2_c.tax_bd=='Ch'),1],X_pca[(smass2_c.tax_bus=='Ch')|(smass2_c.tax_bd=='Ch'),2],c='gold',label='Ch',s=size,alpha=a)
plt.scatter(X_pca[(smass2_c.tax_bus=='Cg')|(smass2_c.tax_bd=='Cg'),1],X_pca[(smass2_c.tax_bus=='Cg')|(smass2_c.tax_bd=='Cg'),2],c='darkgreen',label='Cg',s=size,alpha=a)
plt.scatter(X_pca[(smass2_c.tax_bus=='Cgh')|(smass2_c.tax_bd=='Cgh'),1],X_pca[(smass2_c.tax_bus=='Cgh')|(smass2_c.tax_bd=='Cgh'),2],c='lime',label='Cgh',s=size,alpha=a)
plt.scatter(X_pca[(smass2_c.tax_bus=='B')|(smass2_c.tax_bd=='B'),1],X_pca[(smass2_c.tax_bus=='B')|(smass2_c.tax_bd=='B'),2],c='b',label='B',s=size,alpha=a)
plt.scatter(X_pca[(smass2_c.tax_bus=='Cb')|(smass2_c.tax_bd=='Cb'),1],X_pca[(smass2_c.tax_bus=='Cb')|(smass2_c.tax_bd=='Cb'),2],c='c',label='c',s=size,alpha=a)
#plt.scatter(X_pca[:,1],X_pca[:,2],marker='x',c='gray',alpha=0.3)
plt.legend()
plt.title('PCA: b to x',fontsize=18)
plt.ylabel('C-complex PC3',fontsize=18)
plt.xlabel('C-complex PC2',fontsize=18)

plt.legend(loc=4, fontsize=14)
plt.ylim(0.05,-0.05)
plt.xlim(0.1,-0.1)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)

#plt.savefig('SMASS2_C_PC2-3+met+PB+b_v3.pdf',bbox_inches="tight")
plt.show()