# -*- coding: utf-8 -*-
"""
Created on Tue Mar 12 11:54:09 2024

@author: aleja
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


datos=[]
datos = pd.read_csv('d00_te.csv')

valor_medio = datos.mean(axis=0)
variable_centralizada = datos - valor_medio

desviacion_estandar = datos.std(axis=0)

X = variable_centralizada / desviacion_estandar
Xn=X.to_numpy()
n=Xn.shape[0]
m=Xn.shape[1]

r=1/(n-1)*np.dot(np.transpose(Xn), Xn)

v,delta, vt=np.linalg.svd(r)

delta=np.diag(delta)

Q=np.dot(np.linalg.inv(np.sqrt(delta)),vt)

#z=np.dot(np.dot(np.linalg.inv(np.sqrt(delta)),vt),np.transpose(Xn))
z=np.dot(Q,np.transpose(Xn))

#print(z.shape)

zd= np.zeros_like(z)
for i in range(z.shape[0]):  # Para cada fila
    for j in range(1, z.shape[1]):  # Para cada elemento en la fila
        zd[i, j] = z[i, j] - z[i, j - 1]

zd[:, 0] = z[:, 0]#copias la primera columna, para que no sea 0

rd=1/(n-1)*np.dot(zd,np.transpose(zd))

p,omega,pt=np.linalg.svd(rd)

w=np.dot(p,Q)

s=np.dot(w,np.transpose(Xn))

a=31

wa =w[:31, :]
we=w[31:, :]


pa = p[:31, :]

sa=np.dot(wa,np.transpose(Xn))
sa=np.transpose(sa)
se=np.dot(we,np.transpose(Xn))
se=np.transpose(se)

T2a=[]
for i in range (n):
    
    multi=np.dot(sa[i,:],np.transpose(sa[i,:]))
    
    T2a=np.append(T2a,multi)
    

Taumbral=np.percentile(T2a,99)

T2e=[]
for i in range (n):
    
    multi=np.dot(se[i,:],np.transpose(se[i,:]))
    
    T2e=np.append(T2e,multi)


Teumbral=np.percentile(T2e,99)

ejex=np.arange(n)

plt.title('SFA T2a')
plt.xlabel('t')
plt.ylabel('Ta^2')
plt.scatter(ejex, T2a, c ="green")
plt.axhline(y=Taumbral, color="r", label="Percentil 99")
plt.show()


plt.title('SFA T2e')
plt.xlabel('t')
plt.ylabel('T2e')
plt.scatter(ejex, T2e, c ="blue")
plt.axhline(y=Teumbral, color="r", label="Percentil 99")
plt.show()  




sa=np.transpose(sa)
sad= np.zeros_like(sa)
for i in range(sa.shape[0]):  # Para cada fila
    for j in range(1, sa.shape[1]):  # Para cada elemento en la fila
        sad[i, j] = sa[i, j] - sa[i, j - 1]

sad[:, 0] = sa[:, 0]#copias la primera columna, para que no sea 0


#print(np.shape(sad))
se=np.transpose(se)
sed= np.zeros_like(se)
for i in range(se.shape[0]):  # Para cada fila
    for j in range(1, se.shape[1]):  # Para cada elemento en la fila
        sed[i, j] = se[i, j] - se[i, j - 1]

sed[:, 0] = se[:, 0]#copias la primera columna, para que no sea 0

#print(np.shape(se))

omega=np.diag(omega)
omega_a=(omega[:31,:31])
sad=np.transpose(sad)

S2a=[]
for i in range (n):
    
    multi1=np.dot(np.dot(sad[i,:],np.linalg.inv(omega_a)),np.transpose(sad[i,:]))
      
    S2a=np.append(S2a,multi1)
#print(T2)

S2aumbral=np.percentile(S2a,99)


omega_e=(omega[31:,31:])
sed=np.transpose(sed)

S2e=[]
for i in range (n):
    
    multi1=np.dot(np.dot(sed[i,:],np.linalg.inv(omega_e)),np.transpose(sed[i,:])) 
      
    S2e=np.append(S2e,multi1)
 
S2eumbral=np.percentile(S2e,99)

ejex=np.arange(n)

plt.title('SFA TS2a')
plt.xlabel('t')
plt.ylabel('Sa^2')
plt.scatter(ejex, S2a, c ="green")
plt.axhline(y=S2aumbral, color="r", label="Percentil 99")
plt.show()


plt.title('SFA S2e')
plt.xlabel('t')
plt.ylabel('S2e')
plt.scatter(ejex, S2e, c ="blue")
plt.axhline(y=S2eumbral, color="r", label="Percentil 99")
plt.show()  


np.savez("fich_datos_SFA",wa=wa,we=we,Taumbral=Taumbral,Teumbral=Teumbral,
         valor_medio=valor_medio,desviacion_estandar=desviacion_estandar, w=w,
         omega_a=omega_a,omega_e=omega_e,S2aumbral=S2aumbral,S2eumbral=S2eumbral,Q=Q,p=p,pa=pa)



