%% LABORATORIO 1 - METODI STATISTICI PER LA BIOINGEGNERIA

%% Direttive per la gestione del workspace e delle figure

clearvars  % - cancella tutte le variabili dal workspace
clc        % - Clear Command Window - cancella la command window
close all  % - chiude tutte le figure aperte in precedenza 

%% INTRODUZIONE
%% Indicizzazione 
% Creo una matrice semplice per capire come vengono indicizzati i vettori e
% le matrici in Matlab
mat2D=[1 1 1 1; 1 1 2 1; 1 1 1 1];

% Trovo l'indice dell'elemento uguale a 2
ind = find(mat2D==2);

% Traduco l'indice in formato riga/colonna
[r,c]=ind2sub(size(mat2D),ind);
% N.B. il risultato di ind2sub dipende dall'output che scelgo
ind2sub(size(mat2D),ind)

% Viceversa, a partire dagli indici di riga e colonna posso recuperare
% l'indice lineare che individua l'elemento che sto cercando
sub2ind(size(mat2D),r,c)

%% Vettorizzazione
% A volte e' utile vettorizzare una matrice 2D per effettuare operazioni di
% ricerca o di mascheramento
array1D = mat2D(:);

% Se conosco la dimensione della matrice iniziale, posso inserire gli
% elementi vettorizzati e recupeare l'ordine corretto 

mat2D_restored = zeros(size(mat2D));
mat2D_restored(:) = array1D;

%% Operazioni su vettori e matrici

v1=[4 2 6];
v2=[3 15 7];
% Calcolo la matrice A = v1 trasposto per v2 
% N.B. verificare le dimensioni della matrice A
A = v1'*v2;

% Calcolo la matrice A = v1 per v2 trasposto
A1 = v1*v2';

% Calcolo la norma del vettore v1
normv1 = sqrt(sum(v1.^2));
normv1 - norm(v1)

% Estrazione sottomatrice 2x2 - B = ultime 2 colonne e prima e ultima
% riga di A

B = A([1 end],[end-1 end]);

% Si utilizzino le seguenti matrici A e B

A =[-8     2     7     3
    -8    -9   -10     5
     9    -6   -10     3
    10    -3    -7    -1];

B = [2     5     6     1
    -3    -5    -7     4
     6    -1     0     2
    12     -2    0     4];


% Calcolo del rango
rank(A)
rank(B)

% Calcolo degli autovalori
eig(A)
eig(B)

% Calcolo della matrice D = A per inversa di A
D = A*inv(A)

% N.B. inaccurato, per calcolare A*inv(A) -> A/A
% Cosa dovrebbe risultare?
D1= A/A

% Mostro la differenza
D-D1

% N.B. inaccurato, per calcolare inv(B)*A -> B\A
E= inv(B)*A
E1= B\A

% Mostro la differenza
E-E1

% Calcolo prodotto elemento per elemento
F = A.*B;

%%
% Trovare indici riga e colonna di ogni elemento < 0
ind=find(F<0);

% Mostrare per ogni coppia il posizionamento nella matrice (senza ciclo for)
[r,c]=ind2sub(size(F),ind);

% Mostrare per ogni coppia il posizionamento nella matrice (con ciclo for)
for k=1:length(ind)
    [r,c]=ind2sub(size(F),ind(k))
end

% Operatore max e min
[maxA,indA]=max(A(:))
[r,c]=ind2sub(size(A),indA)

%N.B. cosa succede se non vettorizzo?
[maxA,indA]=max(A);
[r,c]=ind2sub(size(A),indA)

% ripetere per min
[minB,indB]=min(B(:));

%% Generazione di vettori

% intervallo [0,50] con passo 1 usando l'operatore  :
j1 = 0:1:50;
% intervallo [0,50] con passo 0.5 usando l'operatore  :
j2 = 0:.5:50;
% intervallo [501,1000] con passo 1 usando la funzione linspace;
% N.B. passo tra due campioni y=(x2-x1)/(n-1)
passo=1;
x2=1000;
x1=501;
n=(x2-x1)/passo+1;
j3 = linspace(x1,x2,n);

% concatenazione per righe (orizzontale) di j1 e j3 facendo uso delle parentesi quadre o della funzione horzcat
j4 = [j1 j3];
j4 = horzcat(j1,j3);

%  concatenazione per colonne (verticale) di j1 e j3 facendo uso delle
%  parentesi quadre e della funzione vertcat
%  NOTA: non è possibile concatenare  in colonna vettori che non abbiano lo stesso numero di righe. 
% Convertiamo i vettori riga in vettori colonna e ne facciamo il trasposto,
% così da ottenere due vettori colonna, concatenandoli poi verticalmente

j5 = [j1';j3']; %equivalente a j5 = j4'
j5 = vertcat(j1',j3');
% intervallo [1, 100] con 10 punti e con spaziatura logaritmica usando la funzione logspace;
j6=logspace(0,2,10);

%% Esercizio 5

%% Caricamento dati 
load("data_IVGTT.mat")

%% Conversione dati - da matrici a tabelle 

% CERCO MANUALMENTE NELLE VARIABILI QUELLE CHE MI INTERESSANO

% Salvataggio delle sole variabili di interesse
% - estraggo le variabili di interesse Gss_IVGTT (INDICE 7) e Gb_IVGTT
% (INDICE 10)
subset_data = raw_data(:,[7 10]);

% - estraggo le etichette di interesse Gss_IVGTT e Gb_IVGTT
subset_labels = labels([7 10]);

%% Visualizzazione 
% - Scatterplot 
close all
figure(1)
scatter(subset_data(:,1),subset_data(:,2))
xlabel(subset_labels{1,1})
ylabel(subset_labels{1,2})

% - Boxplot - figure separate
for k=1:size(subset_data,2)
    figure(k+1)
    boxplot(subset_data(:,k),'Labels',subset_labels{k})
end

figure(4)
% - Boxplot - stessa figura
for k=1:size(subset_data,2)
    subplot(1,2,k)
    boxplot(subset_data(:,k),'Labels',subset_labels{k})
end

%% Uso del find 
close all
% voglio trovare tutti i soggetti i cui campioni di Glucosio (corrispondenti 
% allo stato stazionario dell'esperimento IVGTT) con valori superiori a 90 [mg/dl] 

indGssHigh = find(subset_data(:,1)>90);
figure(1)

% Mostro anche i valori di Glucosio basale degli stessi soggetti 
boxplot(subset_data(indGssHigh,:),'Labels',subset_labels)

% Alternativa - invece di utilizzare find, utilizzo un array booleano
% (binario) che mi permette di selezionare i dati con operatori
% unari/binari (NOT, AND, OR, XOR)

maskGssHigh = subset_data(:,1)>90;
figure(2)
boxplot(subset_data(maskGssHigh,:),'Labels',subset_labels)
title('Variables of Subjects with Gss > 90 ml/dl')
figure(3)
boxplot(subset_data(not(maskGssHigh),:),'Labels',subset_labels)
title('Variables of Subjects with Gss <= 90 ml/dl')