import pandas as pd
import numpy as np
from sklearn.decomposition import PCA
from sklearn import preprocessing
import matplotlib.pyplot as plt

dataset = pd.read_csv('MC_6Core.csv',index_col=0,decimal=',')  # costruisci dataset da csv

print(dataset.head())
print(dataset.shape)

scaled_data = preprocessing.scale(dataset) #normalizza trasposta del dataset

pca=PCA() #crea oggetto PCA
pca.fit(scaled_data) #fa la pca
pca_data = pca.transform(scaled_data) #prendi Coordinate PCA dai dati normalizzati

##PLOT

#SCREE PLOT
per_var = np.round(pca.explained_variance_ratio_* 100, decimals=1)
labels = ['PC' + str(x) for x in range(1, len(per_var)+1)]

plt.bar(x=range(1,len(per_var)+1), height=per_var, tick_label=labels)
plt.ylabel('Percentage of Explained Variance')
plt.xlabel('Principal Component')
plt.title('Scree Plot')
plt.show()

#the following code makes a fancy looking plot using PC1 and PC2
pca_df = pd.DataFrame(pca_data,index=dataset.index, columns=labels)

plt.figure(figsize=(12, 10))  # Aggiunto per ingrandire il grafico
plt.scatter(pca_df.PC1, pca_df.PC2)
plt.title('PCA MC 6 Core')
plt.xlabel('PC1 - {0}%'.format(per_var[0]))
plt.ylabel('PC2 - {0}%'.format(per_var[1])) 

# Prendi i componenti (autovettori) dalla PCA
components = pca.components_

# Definisci un fattore di scala per i vettori
scaling_factor = 2

# Disegna le frecce per ciascuna misura
for i, (comp_pc1, comp_pc2) in enumerate(zip(components[0], components[1])):
    # Vettore per la misura i-esima
    plt.arrow(0, 0, comp_pc1 * scaling_factor, comp_pc2 * scaling_factor,
              color='red', width=0.01, head_width=0.1, alpha=0.8)
    
    # Etichetta il vettore con il nome della misura
    plt.text(comp_pc1 * scaling_factor * 1.1, comp_pc2 * scaling_factor * 1.1, 
             dataset.columns[i], color='red', ha='center', va='center', fontsize=12)

# Aggiungi le etichette per i punti dati
for sample in pca_df.index:
    plt.annotate(sample, (pca_df.PC1.loc[sample], pca_df.PC2.loc[sample]),
                 textcoords="offset points", xytext=(0,10), ha='center')

# Aggiungi una griglia e linee per gli assi per una migliore visualizzazione
plt.grid(True)
plt.axhline(0, color='grey', linestyle='--')
plt.axvline(0, color='grey', linestyle='--')
plt.axis('equal') # Mantieni la stessa scala per gli assi per evitare distorsioni

plt.show()

## get the name of the top 10 measurements (disturbi) that contribute
## most to pc1.
## first, get the loading scores
loading_scores = pd.Series(pca.components_[0], index=dataset.columns)
## now sort the loading scores based on their magnitude
sorted_loading_scores = loading_scores.abs().sort_values(ascending=False)

# get the names of the top 10 genes
top_metrics = sorted_loading_scores.index.values

## print the gene names and their scores (and +/- sign)
print(loading_scores[top_metrics])