#!/usr/bin/env python
import scipy as sc
import matplotlib.pyplot as plt
import scipy.linalg as la
# 1. Comprimento da sépala em cm
# 2. Largura da sépala em cm
# 3. Comprimento da pétala em cm
# 4. Largura da pétala em cm
# 5. Classe: 1 - Iris Setosa
#            2 - Iris Versicolor
#            3 - Iris Virginica

iris=sc.asarray([[5.1,3.5,1.4,0.2,1],
[4.9,3.0,1.4,0.2,1],
[4.7,3.2,1.3,0.2,1],
[4.6,3.1,1.5,0.2,1],
[5.0,3.6,1.4,0.2,1],
[5.4,3.9,1.7,0.4,1],
[4.6,3.4,1.4,0.3,1],
[5.0,3.4,1.5,0.2,1],
[4.4,2.9,1.4,0.2,1],
[4.9,3.1,1.5,0.1,1],
[5.4,3.7,1.5,0.2,1],
[4.8,3.4,1.6,0.2,1],
[4.8,3.0,1.4,0.1,1],
[4.3,3.0,1.1,0.1,1],
[5.8,4.0,1.2,0.2,1],
[5.7,4.4,1.5,0.4,1],
[5.4,3.9,1.3,0.4,1],
[5.1,3.5,1.4,0.3,1],
[5.7,3.8,1.7,0.3,1],
[5.1,3.8,1.5,0.3,1],
[5.4,3.4,1.7,0.2,1],
[5.1,3.7,1.5,0.4,1],
[4.6,3.6,1.0,0.2,1],
[5.1,3.3,1.7,0.5,1],
[4.8,3.4,1.9,0.2,1],
[5.0,3.0,1.6,0.2,1],
[5.0,3.4,1.6,0.4,1],
[5.2,3.5,1.5,0.2,1],
[5.2,3.4,1.4,0.2,1],
[4.7,3.2,1.6,0.2,1],
[4.8,3.1,1.6,0.2,1],
[5.4,3.4,1.5,0.4,1],
[5.2,4.1,1.5,0.1,1],
[5.5,4.2,1.4,0.2,1],
[4.9,3.1,1.5,0.1,1],
[5.0,3.2,1.2,0.2,1],
[5.5,3.5,1.3,0.2,1],
[4.9,3.1,1.5,0.1,1],
[4.4,3.0,1.3,0.2,1],
[5.1,3.4,1.5,0.2,1],
[5.0,3.5,1.3,0.3,1],
[4.5,2.3,1.3,0.3,1],
[4.4,3.2,1.3,0.2,1],
[5.0,3.5,1.6,0.6,1],
[5.1,3.8,1.9,0.4,1],
[4.8,3.0,1.4,0.3,1],
[5.1,3.8,1.6,0.2,1],
[4.6,3.2,1.4,0.2,1],
[5.3,3.7,1.5,0.2,1],
[5.0,3.3,1.4,0.2,1],
[7.0,3.2,4.7,1.4,2],
[6.4,3.2,4.5,1.5,2],
[6.9,3.1,4.9,1.5,2],
[5.5,2.3,4.0,1.3,2],
[6.5,2.8,4.6,1.5,2],
[5.7,2.8,4.5,1.3,2],
[6.3,3.3,4.7,1.6,2],
[4.9,2.4,3.3,1.0,2],
[6.6,2.9,4.6,1.3,2],
[5.2,2.7,3.9,1.4,2],
[5.0,2.0,3.5,1.0,2],
[5.9,3.0,4.2,1.5,2],
[6.0,2.2,4.0,1.0,2],
[6.1,2.9,4.7,1.4,2],
[5.6,2.9,3.6,1.3,2],
[6.7,3.1,4.4,1.4,2],
[5.6,3.0,4.5,1.5,2],
[5.8,2.7,4.1,1.0,2],
[6.2,2.2,4.5,1.5,2],
[5.6,2.5,3.9,1.1,2],
[5.9,3.2,4.8,1.8,2],
[6.1,2.8,4.0,1.3,2],
[6.3,2.5,4.9,1.5,2],
[6.1,2.8,4.7,1.2,2],
[6.4,2.9,4.3,1.3,2],
[6.6,3.0,4.4,1.4,2],
[6.8,2.8,4.8,1.4,2],
[6.7,3.0,5.0,1.7,2],
[6.0,2.9,4.5,1.5,2],
[5.7,2.6,3.5,1.0,2],
[5.5,2.4,3.8,1.1,2],
[5.5,2.4,3.7,1.0,2],
[5.8,2.7,3.9,1.2,2],
[6.0,2.7,5.1,1.6,2],
[5.4,3.0,4.5,1.5,2],
[6.0,3.4,4.5,1.6,2],
[6.7,3.1,4.7,1.5,2],
[6.3,2.3,4.4,1.3,2],
[5.6,3.0,4.1,1.3,2],
[5.5,2.5,4.0,1.3,2],
[5.5,2.6,4.4,1.2,2],
[6.1,3.0,4.6,1.4,2],
[5.8,2.6,4.0,1.2,2],
[5.0,2.3,3.3,1.0,2],
[5.6,2.7,4.2,1.3,2],
[5.7,3.0,4.2,1.2,2],
[5.7,2.9,4.2,1.3,2],
[6.2,2.9,4.3,1.3,2],
[5.1,2.5,3.0,1.1,2],
[5.7,2.8,4.1,1.3,2],
[6.3,3.3,6.0,2.5,3],
[5.8,2.7,5.1,1.9,3],
[7.1,3.0,5.9,2.1,3],
[6.3,2.9,5.6,1.8,3],
[6.5,3.0,5.8,2.2,3],
[7.6,3.0,6.6,2.1,3],
[4.9,2.5,4.5,1.7,3],
[7.3,2.9,6.3,1.8,3],
[6.7,2.5,5.8,1.8,3],
[7.2,3.6,6.1,2.5,3],
[6.5,3.2,5.1,2.0,3],
[6.4,2.7,5.3,1.9,3],
[6.8,3.0,5.5,2.1,3],
[5.7,2.5,5.0,2.0,3],
[5.8,2.8,5.1,2.4,3],
[6.4,3.2,5.3,2.3,3],
[6.5,3.0,5.5,1.8,3],
[7.7,3.8,6.7,2.2,3],
[7.7,2.6,6.9,2.3,3],
[6.0,2.2,5.0,1.5,3],
[6.9,3.2,5.7,2.3,3],
[5.6,2.8,4.9,2.0,3],
[7.7,2.8,6.7,2.0,3],
[6.3,2.7,4.9,1.8,3],
[6.7,3.3,5.7,2.1,3],
[7.2,3.2,6.0,1.8,3],
[6.2,2.8,4.8,1.8,3],
[6.1,3.0,4.9,1.8,3],
[6.4,2.8,5.6,2.1,3],
[7.2,3.0,5.8,1.6,3],
[7.4,2.8,6.1,1.9,3],
[7.9,3.8,6.4,2.0,3],
[6.4,2.8,5.6,2.2,3],
[6.3,2.8,5.1,1.5,3],
[6.1,2.6,5.6,1.4,3],
[7.7,3.0,6.1,2.3,3],
[6.3,3.4,5.6,2.4,3],
[6.4,3.1,5.5,1.8,3],
[6.0,3.0,4.8,1.8,3],
[6.9,3.1,5.4,2.1,3],
[6.7,3.1,5.6,2.4,3],
[6.9,3.1,5.1,2.3,3],
[5.8,2.7,5.1,1.9,3],
[6.8,3.2,5.9,2.3,3],
[6.7,3.3,5.7,2.5,3],
[6.7,3.0,5.2,2.3,3],
[6.3,2.5,5.0,1.9,3],
[6.5,3.0,5.2,2.0,3],
[6.2,3.4,5.4,2.3,3],
[5.9,3.0,5.1,1.8,3]],dtype=float)

X=iris[:,0:4]
[r,c]=sc.shape(X) #calcula numero de r linhas e c colunas de X
Mx=sc.mean(X); #uma linha de c medias, para cada coluna de X
#Xm=X-sc.dot(sc.ones([r,1]),Mx) #ones(r,1)*Mx, para termos r linhas de Mx
Xm=X;
print('Matriz de covariância\n');
C=sc.cov(sc.transpose(Xm)) #Xm.T é o transposta de Xm, o python vê matematicamente
print(C)
[D,V]=la.eigh(C) #,eigvals_only=True) #autovalores e autovetores da matriz de covariância de Xm
i=sc.argsort(D)[::-1]  # Obtém índices para ordenação decrescente dos autovalores
l=D[i]
print('\nAutovalores ordenados\n')
print(l)
print('\nAutovetores\n')
print(V[:,i])
print('\nPorcentagem explicativa acumulada das variáveis\n')
pe=sc.cumsum(l)/sc.sum(l) #porcentagem explicativa acumulada de cada eixo
print(pe.T)
O=V[:,i[0:2]]            #toma os dois eixos principais de acordo com a ordem dos autovalores descentes

#==================================================================================%
#Veja que matematicamente, Y=O'X. Na visão de BD Y está transposto. Logo,          %
#para vermos como queremos (visão BD) temos que transpor Y. Assim, Y'=X'O''= X'O.  %
#Como temos X em notação de BD, o nosso Y em notação de BD será Y=XO               %
#==================================================================================%

#Y=sc.dot(Xm,O) Xm rotacionado pelos eixos principais.
Y=Xm@O
#plota 50 tipos de cada especie com suas cores r,g,b, tamanho 2

plt.plot(Y[0:50,0],Y[0:50,1],'or',Y[50:100,0],Y[50:100,1],'og',Y[100:150,0],Y[100:150,1],'ob')

plt.title("Rotação dos eixos principais de X. Ou seja, gerar um Y=O'X") #título do gráfico
plt.xlabel('Eixo 1') #nome para eixo x
plt.ylabel('Eixo 2') #nome para eixo y
plt.legend(['Setosa','Versicolor', 'Virginica'],loc='best') #coloca legendas no gráfico
#plt.legend('location','best') #na melhor posição
plt.show()