# -*- coding: utf-8 -*-
"""
Created on Thu Apr 18 18:05:05 2024

@author: aleja
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import skew
from scipy.stats import norm, kurtosis
from scipy.stats import f , chi2

datos=[]

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



#############por columnas#########################

conjunto_muestras = []
conjunto_medias=[]
conjunto_desviacion_estandar=[]
conjunto_correlacion=[]
conjunto_skewness=[]
conjunto_kurtosis=[]
kv=50

for i in range(1, n-(kv-1)):
    muestra = datos.iloc[i:i+kv, :] #el iloc es por los datos en panda
    muestra = np.array(muestra, dtype=np.float64)  # convert to numpy array of float64
    conjunto_muestras.append(muestra)  # append the numpy array

    media = np.mean(muestra,axis=0)  # calculate mean of the entire matrix
    conjunto_medias.append(media)  # append the mean to the list

    desviacion_estandar=np.std(muestra,axis=0)
    desviacion_estandar=desviacion_estandar**2
    conjunto_desviacion_estandar.append(desviacion_estandar)  
    
    
   # correlacion = np.corrcoef(muestra) 
    #tri_sup = np.triu(correlacion, k=1)
    #tri_elementos = tri_sup[tri_sup > 0.3]
    #conjunto_correlacion.append(tri_elementos)
    
    skewness = skew(muestra,axis=0)
    conjunto_skewness.append(skewness)
    
    kurtosiss = kurtosis(muestra, axis=0)
    conjunto_kurtosis.append(kurtosiss)
    
#para que tengan forma array of float64
conjunto_muestras = np.concatenate(conjunto_muestras, axis=0)#para comprobar
conjunto_medias=np.array(conjunto_medias,dtype=np.float64)
conjunto_desviacion_estandar=np.array(conjunto_desviacion_estandar,dtype=np.float64)
#conjunto_correlacion=np.array(conjunto_correlacion,dtype=np.float64)
conjunto_skewness=np.array(conjunto_skewness,dtype=np.float64)
conjunto_kurtosis=np.array(conjunto_kurtosis,dtype=np.float64)



# Crea una lista de las 4 matrices
matrices = [conjunto_medias, conjunto_desviacion_estandar,conjunto_skewness,conjunto_kurtosis]

# Para intercalar, tomaremos una columna de cada matriz a la vez y las concatenaremos
#datos_nueva = np.hstack([matrices[i][:, j].reshape(n-kv, 1) for j in range(52) for i in range(4)])
datos_nueva=np.concatenate((conjunto_medias, conjunto_desviacion_estandar,conjunto_skewness), axis=1)
datos_nueva=np.concatenate((datos_nueva,conjunto_kurtosis), axis=1)
#datos_nueva=np.concatenate((datos_nueva,conjunto_kurtosis,conjunto_correlacion), axis=1)
n,m=datos_nueva.shape
########ahora hacemos PCA de forma normal


#2 normalizar a media conocida y varianza 1
    
valor_medio = datos_nueva.mean(axis=0)
variable_centralizada = datos_nueva - valor_medio

desviacion_estandar = datos_nueva.std(axis=0)

Xn = variable_centralizada / desviacion_estandar

#3 cakcualar la covarianza de x

r=1/(n-1)*np.dot(np.transpose(Xn), Xn)    #52x52 es el tamaño que buscamos
#r=np.cov                                   #960x960



#4 calcular los valores propios, valores singulares de la matriz

v,s, vt=np.linalg.svd(r) #u mat vect sing izq, s valores singulares R ordenados enforma descend, v vec sing derec

#5 componentes principales

t=np.dot(Xn,v)

#print(np.shape(t)) 


variabilidad=90
variabilidadc=0
a=0

traza=sum(s)

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

#print('hace ',a,'iteraciones') 

p=v[:,:a]              #matriz de tamaño xa
sa=np.diag(s[:a])

ta=np.dot(Xn, p)
#print(np.shape(ta))    #960xa



T2=[]
for i in range (n):
    
    multi=np.dot(np.dot(np.dot(np.dot(Xn[i,:], p), np.linalg.inv(sa)), np.transpose(p)), np.transpose(Xn[i,:]))

    T2=np.append(T2,multi)
#print(T2)

Tumbral=np.percentile(T2,99.9)
#print(Tumbral)


#calcular la estadistica Q
Q=[]

for i in range(n):
    R=Xn[i,:]-np.dot(np.dot(Xn[i,:],p),np.transpose(p))
    Rt=np.transpose(R) #para que de la forma adecuada
    Qi=np.dot(Rt,R)
    Q=np.append(Q,Qi)
    
#print(np.shape(Q))
#print(Q)

Qumbral=np.percentile(Q,99.9)


mu = np.mean(Q)
varQ = np.std(Q)#duda esto es la varianza o la desviación tipica??
h = 2 * (mu ** 2) / varQ
g =varQ / (2*mu)
QA = g * chi2.ppf(0.99 , h)
#hacer graficas

ejex=np.arange(n)

plt.title('SPA 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()


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


###guardar datos en un fichero

np.savez("fich_datos_SPA",p=p,sa=sa,Qumbral=Qumbral,Tumbral=Tumbral,valor_medio=valor_medio,desviacion_estandar=desviacion_estandar,kv=kv)




    
    
