# -*- coding: utf-8 -*-
"""
Created on Mon Jun 17 16:03:03 2024

@author: aleja
"""


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

#subir un poco umbral Q


variables = joblib.load('fich_datos_DSPCA.pkl')

Qumbral=variables['Qumbral']
Tumbral=variables['Tumbral']
valor_medio=variables['valor_medio']
desviacion_estandar=variables['desviacion_estandar']
var_t=variables['var_t']
var_res=variables['var_res']
p=variables['p']

model_spca=variables['model_spca']


datos = pd.read_csv('d12_te.csv')

shape = datos.shape
n=len(datos.axes[0])
m=len(datos.axes[1])


variable_centralizada = datos - valor_medio

X = variable_centralizada / desviacion_estandar

Xn=X.to_numpy()
Xn_1 = np.roll(Xn, 1, axis=1)#desplazamientos de los elementos de la matriz
Xn_1[:, 0] = Xn[:, 0]

Xn_2 = np.roll(Xn_1, 1, axis=1)
Xn_2[:, 0] = Xn_1[:, 0]


Xn_3 = np.roll(Xn_2, 1, axis=1)
Xn_3[:, 0] = Xn_2[:, 0]


XD = np.concatenate((Xn,Xn_1, Xn_2, Xn_3), axis=1)


n_components=31

t_spca = model_spca.transform(XD)
P_spca=model_spca.components_.T
Xn_spca = model_spca.inverse_transform(t_spca)

T2=[]
for i in range (n):
    
    multi=(np.dot(np.dot(t_spca[i,:],np.linalg.pinv(var_t)), np.transpose(t_spca[i,:])))

    T2=np.append(T2,multi)


residuo=Xn_spca-XD

Q=[]
for i in range (n):
    
    multi1=(np.dot(np.dot(residuo[i,:],np.linalg.pinv(var_res)), np.transpose(residuo[i,:])))

    Q=np.append(Q,multi1)



ejex=np.arange(n)
plt.title('Estadistico T^2 de Hotelling')
plt.xlabel('t')
plt.ylabel('T^2')
plt.scatter(ejex, T2, c ="green")
plt.axhline(y=Tumbral, color="r", label="Percentil 99")
plt.show()

ejex=np.arange(n)
plt.title('Estadistico Q')
plt.xlabel('t')
plt.ylabel('Q')
plt.scatter(ejex, Q, c ="blue")
plt.axhline(y=Qumbral, color="r", label="Percentil 99")
plt.show()




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

  
con=0    
t_falloQ=0
for i in range (0,n):
    if Q[i]>Qumbral:
         con=con+1
    else:
         con=0
        
    if con==k:
        fallo =1
        t_falloQ=i-9
        break
    
#print('tiempo para Q',t_falloQ)
#falsas alarmas
cont=0
for i in range(159):
    if T2[i]>Tumbral:
        cont=cont+1

falsa_alarmaT=cont/160*100

print('falsas alarmas para T2',falsa_alarmaT)

cont2=0
for i in range(159):
    if Q[i]>Qumbral:
        cont2=cont2+1

falsa_alarmaQ=cont2/160*100

#print('falsas alarmas para Q',falsa_alarmaQ)

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

alarma_detectadaT=con/800*100
print(' alarmas detectadas para T2',alarma_detectadaT)

con2=0
for i in range (160,n):
    if Q[i]>Qumbral:
        con2=con2+1

alarma_detectadaQ=con2/800*100
#print(' alarmas detectadas para Q',alarma_detectadaQ)


print('tiempo para Q',t_falloQ)
print('falsas alarmas para Q',falsa_alarmaQ)
print(' alarmas detectadas para Q',alarma_detectadaQ)
#################### identificacion fallo

#R_fallo=XD[t_falloQ,:]
R_fallo=XD[t_falloQ,:]-np.dot(np.dot(XD[t_falloQ,:],p),np.transpose(p))

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

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










