# -*- coding: utf-8 -*-
"""
Created on Thu Mar 21 16:17:28 2024

@author: aleja
"""

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

variables = np.load('fich_datos_SFA.npz')

wa=variables['wa']
we=variables['we']
Taumbral=variables['Taumbral']
Teumbral=variables['Teumbral']
valor_medio=variables['valor_medio']
desviacion_estandar=variables['desviacion_estandar']
omega_a=variables['omega_a']
omega_e=variables['omega_e']
S2aumbral=variables['S2aumbral']
S2eumbral=variables['S2eumbral']
Q=variables['Q']
p=variables['p']
pa=variables['pa']
datos = pd.read_csv('d21_te.csv')



#2 normalizar a media conocida y varianza 1
variable_centralizada = datos - valor_medio

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



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)




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

ejex=np.arange(n)

plt.title('Estadistico T^2a de Hotelling')
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('Estadistico T^2e de Hotelling')
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))

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)


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)
 

ejex=np.arange(n)

plt.title('Estadistico S2a')
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('Estadistico S2e')
plt.xlabel('t')
plt.ylabel('S2e')
plt.scatter(ejex, S2e, c ="blue")
plt.axhline(y=S2eumbral, color="r", label="Percentil 99")
plt.show()  



#detectar fallos
t_falloT2a=0
con=0
k=10
fallo=0
for i in range (0,n):
    if T2a[i]>Taumbral:     #para Ta
        con=con+1
    else:
        con=0
        
    if con==k:
        fallo =1
        t_falloT2a=i-9
        break
    
    
#print('tiempo fallo para T2a',t_falloT2a)  


t_falloT2e=0
con=0
k=10
fallo=0
for i in range (0,n):
    if T2e[i]>Teumbral:     #para Te
        con=con+1
    else:
        con=0
        
    if con==k:
        fallo =1
        t_falloT2e=i-9
        break
    
    
#print('tiempo fallo para T2e',t_falloT2e)  


t_falloS2a=0
con=0
k=10
fallo=0
for i in range (0,n):
    if S2a[i]>S2aumbral:     #para Sa
        con=con+1
    else:
        con=0
        
    if con==k:
        fallo =1
        t_falloS2a=i-9
        break
    
    
#print('tiempo fallo para S2a',t_falloS2a)  


t_falloS2e=0
con=0
k=10
fallo=0
for i in range (0,n):
    if S2e[i]>S2eumbral:     #para Se
        con=con+1
    else:
        con=0
        
    if con==k:
        fallo =1
        t_falloS2e=i-9
        break
    
    
#print('tiempo fallo para S2e',t_falloS2e)  



#falsas alarmas
cont=0
for i in range(159):
    if T2a[i]>Taumbral:
        cont=cont+1

falsa_alarmaTa=cont/160*100

#print('falsas alarmas para T2a',falsa_alarmaTa)

cont=0
for i in range(159):
    if T2e[i]>Teumbral:
        cont=cont+1

falsa_alarmaTe=cont/160*100

#print('falsas alarmas para T2e',falsa_alarmaTe)

#ahora para las S2
cont=0
for i in range(159):
    if S2a[i]>S2aumbral:
        cont=cont+1

falsa_alarmaS2a=cont/160*100

#print('falsas alarmas para S2a',falsa_alarmaS2a)

cont=0
for i in range(159):
    if T2e[i]>Teumbral:
        cont=cont+1

falsa_alarmaS2e=cont/160*100

#print('falsas alarmas para S2e',falsa_alarmaS2e)




#alarmas detectadas
con=0
for i in range (160,n):
    if T2a[i]>Taumbral:
        con=con+1

alarma_detectadaTa=con/800*100
#print(' alarmas detectadas para T2a',alarma_detectadaTa)

con=0
for i in range (160,n):
    if T2e[i]>Teumbral:
        con=con+1

alarma_detectadaTe=con/800*100
#print(' alarmas detectadas para T2e',alarma_detectadaTe)

con=0
for i in range (160,n):
    if S2a[i]>S2aumbral:
        con=con+1

alarma_detectadaS2a=con/800*100
#print(' alarmas detectadas para S2a',alarma_detectadaS2a)

con=0
for i in range (160,n):
    if S2e[i]>S2eumbral:
        con=con+1

alarma_detectadaS2e=con/800*100
#print(' alarmas detectadas para S2e',alarma_detectadaS2e)


################3
print('tiempo fallo para T2a',t_falloT2a)  
print('falsas alarmas para T2a',falsa_alarmaTa)
print(' alarmas detectadas para T2a',alarma_detectadaTa)
print()
print('tiempo fallo para T2e',t_falloT2e)  
print('falsas alarmas para T2e',falsa_alarmaTe)
print(' alarmas detectadas para T2e',alarma_detectadaTe)
print()
print('tiempo fallo para S2a',t_falloS2a)  
print('falsas alarmas para S2a',falsa_alarmaS2a)
print(' alarmas detectadas para S2a',alarma_detectadaS2a)
print()
print('tiempo fallo para S2e',t_falloS2e)  
print('falsas alarmas para S2e',falsa_alarmaS2e)
print(' alarmas detectadas para S2e',alarma_detectadaS2e)

#diagrama de barras

#R_fallo=Xn[t_falloS2a,:]
Xn_hat=np.dot(np.dot(np.linalg.inv(Q), np.transpose(pa)),sa);
Xn_hat=np.transpose(Xn_hat)
R_fallo= Xn[t_falloT2e,:]-Xn_hat[t_falloT2e,:];

 

R_fallo2=np.square(R_fallo)
plt.bar(range(len(R_fallo2)), R_fallo2)
plt.title('SFA Diagrama de barras de R_fallo2 d12')
plt.show()


max_indice = np.argmax(R_fallo2)
print(f'posición {max_indice}')



