brain_glmm.m

%%
clear all
close all
clc
%% Add path
path = genpath('/Users/iricchi/Documents/MIP_LabImmersion');
addpath(path)
pos = 'server';

%% set what to run
figs = 0;
subj = 50;
sihl = 0;

%% Parameters to run the code
Ns = 50; % 20, 1 - number of subjects
 
% spatial normalization 
sp_norm = 1;

% type of analysis: choose between 'all', 'persubj', 'perses'
% N.B. perses considers only 1 subj
an_t = 'all'; 
%% load dataset
atlas = 'AAL';
Data_ = load_dataset(pos,atlas);

% Normalize entirely the dataset (all subj) spatially
if (sp_norm ==1)
    Data = Data_'; % along the spatial dimension
    Data = zscore(Data);
    Data = Data';
else
    Data = zscore(Data_);
end

%% Normalize the data per subject
val = length(Data)/Ns;
%data = zeros(size(Data));
%for i = 1:Ns
%    datan = zscore(Data(1+(i-1)*val:i*val,:));
%    data(1+(i-1)*val:i*val,:) = datan;
%end
%


%% generate graphs
m = size(Data,1); %graph size
n = size(Data,2); %number of signals
k = 7; %number of clusters


%% train a glmm on data y
iterations = 50;
%[L1_its, L2_its, gamma_its, mu1, mu2, log_likelihood1] = glmm(y, iterations);

% pick one
%[Ls, gamma_hats, mus, log_likelihood] = glmm_multi_smooth(y', iterations,k);
norm_par = 1.2; % 2 per subj, 1.2 all
delta = 10; % 25 per subj, 10 all
%spread = 0.1;

if(strcmp(an_t,'persubj'))
   
    Ls_t = zeros(n,n,k,Ns);
    gamma_hats_t = zeros(val,k,Ns);
    mus_t = zeros(n,k,Ns);
    for i = 0:Ns-1
%         if (i == 2 || i == 10 || i ==20 || i == 35 || i ==39 || i ==46)
%             delta = 11;
%         else
%             delta = 20;
%         end
       data = Data(val*i+1:val*(i+1),:);
       try
           [Ls_t(:,:,:,i+1), gamma_hats_t(:,:,i+1), mus_t(:,:,i+1)] = glmm_final(data,iterations,k,[],[],norm_par,delta);
       catch ME
           disp(['Subject ', num2str(i+1), ' did not work!'])
           warning(ME.message)
       end
       
    end
    
    save('Ls_t_RS_opt', 'Ls_t')
    save('gamma_hats_t_opt', 'gamma_hats_t')
    save('mus_t_RS_opt', 'mus_t')
    
elseif(strcmp(an_t,'perses'))
    Ls_s1_s = zeros(n,n,k,2);  % 2 sessions
    gamma_hats_s1_s = zeros(val/2,k,2);
    mus_s1_s = zeros(n,k,2);
    
    for x = 1:2
        data = Data(val/2*(x-1)+1:val/2*x,:);
        [Ls_s1_s(:,:,:,x), gamma_hats_s1_s(:,:,x), mus_s1_s(:,:,x)] = glmm_final(data,iterations,k,[],[],norm_par,delta);
    end
    
    
    
elseif(strcmp(an_t,'all'))
    [Ls, gamma_hats, mus] = glmm_final(Data,iterations, k, [], [], norm_par,delta);
    
    disp(sum(gamma_hats,1));
    
    %disp(sum(gamma_its,2));
    
else
    disp('Parameter an_t should be "all", "persubj" or "perses"!')
end
% 

% for i = 1:k
%     [precision(i), recall(i), f(i), NMI_score(i), num_of_edges(i)] = graph_learning_perf_eval(squeeze(Lap(:,:,identify(i))), squeeze(Ls(:,:,i)));
% end
% 
% if (norm(mus(:,1)' - center1,'fro') < norm(mus(:,2)' - center1,'fro'))
%     [precision1, recall1, f1, NMI_score1, num_of_edges1] = graph_learning_perf_eval(Lap1, Ls(:,:,1));
%     [precision2, recall2, f2, NMI_score2, num_of_edges2] = graph_learning_perf_eval(Lap2, Ls(:,:,2));
% else
%     [precision1, recall1, f1, NMI_score1, num_of_edges1] = graph_learning_perf_eval(Lap2, Ls(:,:,1));
%     [precision2, recall2, f2, NMI_score2, num_of_edges2] = graph_learning_perf_eval(Lap1, Ls(:,:,2));
% end
% if params not good: [ optPrec, optRec] = precisionRecall(Lap1,L1_its(:,:,10))



%%
%figure(1)
%plot(log_likelihood, 'b')
%hold on
%plot(log_likelihood1, 'r')

%%
if(figs == 1)
    figure(2)
    %subplot(1,2,1)
    %param.nb_pts = 360;
    [G1]=gsp_sphere(n);
    L1 = Ls(:,:,1);
    W1 = diag(diag(L1)) - L1;
    W1(W1<1e-4) = 0;
    L1 = diag(sum(W1)) - W1;
    G1.L = L1;
    G1.W = W1;
    %G1.W(mu1<0.4, mu1<0.4) = 0; %plotting only for values that are not always 0 (because those will be strongly correlated, but meaningless
    %G1.W(W1<0.2) = 0;
    gsp_plot_graph(G1)

    %%
    %subplot(1,2,2)
    L2 = Ls(:,:,5);
    W2 = diag(diag(L2)) - L2;
    W2(W2<1e-4) = 0;
    L2 = diag(sum(W2)) - W2;
    %G.W(mu2<0.4, mu2<0.4) = 0; %plotting only for values that are not always 0 (because those will be strongly correlated, but meaningless
    %G2.W(W2<0.9) = 0;
    %gsp_plot_graph(G)
    figure()
    imagesc(W2)
    title(["W2 for norm = ", norm_par, "delta = ", delta])

    %%

    L2 = Ls(:,:,20);
    W2 = diag(diag(L2)) - L2;
    W2(W2<1e-4) = 0;
    L2 = diag(sum(W2)) - W2;
    figure()
    imagesc(W2)
    colorbar()
end
%% imshow the gammas

if (strcmp(an_t,'all')) 

    figure()
    imagesc(gamma_hats)
    savefig(['gammas_delta',num2str(delta),'norm', num2str(norm_par),'.fig'])
    save('gamma_hats', 'gamma_hats')
    save('Ls', 'Ls')
    save('mus','mus')
end