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

@author: aleja
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import SparsePCA
import joblib

datos=[]
#1 abrir archivo.csv
 

datos = pd.read_csv('d00_te.csv')
shape = datos.shape
n=len(datos.axes[0])
m=len(datos.axes[1])


#2 normalizar a media conocida y varianza 1
    
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()

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

model_spca=SparsePCA(n_components=31, alpha=0.1) #cambiar el alpha  a 0.1
model_spca.fit(XD)

t_spca = model_spca.transform(XD)
var_t=np.cov(t_spca.T)

P_spca=model_spca.components_.T

Xn_spca = model_spca.inverse_transform(t_spca)

#r=1/(n-1)*np.dot(np.transpose(Xn), Xn) 
r=1/(n-1)*np.dot(np.transpose(XD), XD)
v,s, vt=np.linalg.svd(r)
a=0
variabilidad=90
variabilidadc=0
traza=sum(s)

while (variabilidadc < variabilidad and a <len(s)):
    variabilidadc += (s[a]/traza)*100
    a += 1

p = v[:, :a]

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)


Tumbral=np.percentile(T2,99)

residuo=Xn_spca-XD

var_res=np.cov(residuo.T)

Q=[]
for i in range (n):  #algunas de las matrices son singulares por eso pinv
    
    multi1=(np.dot(np.dot(residuo[i,:],np.linalg.pinv(var_res)), np.transpose(residuo[i,:])))

    Q=np.append(Q,multi1)




Qumbral=np.percentile(Q,99)

ejex=np.arange(n)
plt.title('PCA T2')
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('PCA T2')
plt.xlabel('t')
plt.ylabel('Q')
plt.scatter(ejex, Q, c ="blue")
plt.axhline(y=Qumbral, color="r", label="Percentil 99")
plt.show()


 
 
joblib.dump({
    'Qumbral': Qumbral,
    'P_spca':P_spca,
    'p':p,
    'Tumbral': Tumbral,
    'valor_medio': valor_medio,
    'desviacion_estandar': desviacion_estandar,
    'var_t':var_t,
    'var_res':var_res,
    'model_spca': model_spca
}, 'fich_datos_DSPCA.pkl')

