%% AAL3 Summarization of voxel-based results - Requires MATLAB version R19b and newer

clear; clc

% You need to change this to the path where you have the voxel-based framework folder
framework_to_run = 'C:\Users\Konstantinos Tsirlis\Dropbox\My Desktop\Mourao - Miranda\ABCD results\abcdAll_noA\framework\spls_fair20-5_fair_modeA_pub_noAnom_new_common_bmi';


% You need to change this to the path where you have the warp_atlas folder
cd('C:\Users\Konstantinos Tsirlis\Dropbox\My Desktop\Mourao - Miranda\Atlases_Mikael\warp_atlas\result')
aal3 = spm_vol('raal_AAL3v1_1mm.nii');
aal3 = spm_read_vols(aal3);
aal3 = reshape(aal3,1,[]);


% First you have to go inside the framework/model/res folder that you want.
cd([framework_to_run '\res']);

!del summ_table_aal3.xlsx



for kk = 1:length(dir('level*'))


    
cd(['level' num2str(kk)])



lala = dir('weightX_*');
lala = lala.name;
weights = spm_vol(lala);
weights = spm_read_vols(weights);
weights = reshape(weights,1,[]);
lala = dir('brain_*');
lala = lala.name;
loadings = spm_vol(lala);
loadings = spm_read_vols(loadings);
loadings = reshape(loadings,1,[]);




% areas of (1) non-zero weights, (2) positive weights, (3) negative weights,
% (4) net weight
for i = 1:max(aal3)
dummy = aal3==i;
aal3_weights{i,1} = weights(dummy);
aal3_weights{i,2} = aal3_weights{i,1}(:,aal3_weights{i,1}~=0);
aal3_weights{i,3} = length(aal3_weights{i,2});
aal3_weights{i,4} = aal3_weights{i,1}(:,aal3_weights{i,1}>0);
aal3_weights{i,5} = length(aal3_weights{i,4});
aal3_weights{i,6} = aal3_weights{i,1}(:,aal3_weights{i,1}<0);
aal3_weights{i,7} = length(aal3_weights{i,6});

aal3_loadings{i,1} = loadings(dummy);
aal3_loadings{i,2} = aal3_loadings{i,1}(:,aal3_loadings{i,1}~=0);
aal3_loadings{i,3} = length(aal3_loadings{i,2});
aal3_loadings{i,4} = aal3_loadings{i,1}(:,aal3_loadings{i,1}>0);
aal3_loadings{i,5} = length(aal3_loadings{i,4});
aal3_loadings{i,6} = aal3_loadings{i,1}(:,aal3_loadings{i,1}<0);
aal3_loadings{i,7} = length(aal3_loadings{i,6});

end


all_non_zero_weights = sum(cell2mat(aal3_weights(:,3)));
all_non_zero_loadings = sum(cell2mat(aal3_loadings(:,3)));


lala = [];
for i = 1:max(aal3)  
lala = [lala aal3_weights{i,2}];
end
net_total_weights = sum(abs(lala));
net_positive_weights = sum(lala(lala>0));
flag_pos_weights = 0;
if net_positive_weights==0
    flag_pos_weights = 1;
end
net_negative_weights = sum(lala(lala<0));
flag_neg_weights = 0;
if net_negative_weights==0
    flag_neg_weights = 1;
end

lala = [];
for i = 1:max(aal3)
lala = [lala aal3_loadings{i,2}];
end
net_total_loadings = sum(abs(lala));
net_positive_loadings = sum(lala(lala>0));
flag_pos_loadings = 0;
if net_positive_loadings==0
    flag_pos_loadings = 1;
end
net_negative_loadings = sum(lala(lala<0));
flag_neg_loadings = 0;
if net_negative_loadings==0
    flag_neg_loadings = 1;
end

% Change the path accordingly
load('C:\Users\Konstantinos Tsirlis\Dropbox\My Desktop\Mourao - Miranda\Atlases_Mikael\warp_atlas\katia_aal3.mat')



for i = 1:max(aal3)
summ_table_weights(i,1) = length(aal3_weights{i,1});
summ_table_weights(i,2) = 100*aal3_weights{i,3}/summ_table_weights(i,1); % proportion of the region
summ_table_weights(i,3) = 100*sum(abs(aal3_weights{i,2}))/net_total_weights;
if flag_pos_weights
    summ_table_weights(i,4) = 0;
else
    summ_table_weights(i,4) = 100*sum(aal3_weights{i,4})/net_positive_weights;
end
if flag_neg_weights
    summ_table_weights(i,5) = 0;
else
    summ_table_weights(i,5) = 100*sum(aal3_weights{i,6})/net_negative_weights;
end
end
summ_table_weights(isnan(summ_table_weights(:,2)),2)=0;


for i = 1:max(aal3)
summ_table_loadings(i,1) = length(aal3_loadings{i,1});
summ_table_loadings(i,2) = 100*aal3_loadings{i,3}/summ_table_loadings(i,1); % proportion of the region
summ_table_loadings(i,3) = 100*sum(abs(aal3_loadings{i,2}))/net_total_loadings;
if flag_pos_loadings
    summ_table_loadings(i,4) = 0;
else
    summ_table_loadings(i,4) = 100*sum(aal3_loadings{i,4})/net_positive_loadings;
end
if flag_neg_loadings
    summ_table_loadings(i,5) = 0;
else
    summ_table_loadings(i,5) = 100*sum(aal3_loadings{i,6})/net_negative_loadings;
end
end
summ_table_loadings(isnan(summ_table_loadings(:,2)),2)=0;





%%% THE 2 ADDED COLUMNS
summ_table_weights(:,6) = summ_table_weights(:,4)./summ_table_weights(:,1);
summ_table_weights(:,7) = summ_table_weights(:,5)./summ_table_weights(:,1);
summ_table_loadings(:,6) = summ_table_loadings(:,4)./summ_table_loadings(:,1);
summ_table_loadings(:,7) = summ_table_loadings(:,5)./summ_table_loadings(:,1);






aal3_areas = {katia.textdata{:,2}}';
picu = {'Region name' 'Number of voxels in region' '% of non-zeros weights in the regions' '% of total weights in the regions' '% of total positive weights in the regions' '% of total negative weights in the regions' '% of total positive weights in the regions/no. of voxels' '% of total negative weights in the regions/no. of voxels'};



writecell(picu,[framework_to_run '\res\summ_table_aal3.xlsx'],'Sheet',['Weights_level' num2str(kk)]);
writecell(aal3_areas,[framework_to_run '\res\summ_table_aal3.xlsx'],'Range','A2','Sheet',['Weights_level' num2str(kk)]);
writematrix(summ_table_weights,[framework_to_run '\res\summ_table_aal3.xlsx'],'Range','B2','Sheet',['Weights_level' num2str(kk)]);

writecell(picu,[framework_to_run '\res\summ_table_aal3.xlsx'],'Sheet',['Loadings_level' num2str(kk)]);
writecell(aal3_areas,[framework_to_run '\res\summ_table_aal3.xlsx'],'Range','A2','Sheet',['Loadings_level' num2str(kk)]);
writematrix(summ_table_loadings,[framework_to_run '\res\summ_table_aal3.xlsx'],'Range','B2','Sheet',['Loadings_level' num2str(kk)]);


% save(['aal3_areas_level' num2str(kk)],'aal3_areas')
clear aal3_weights summ_table_weights aal3_loadings summ_table_loadings




cd ..
end







%% Schaefer-Buck-Choi Summarization of voxel-based results - Requires MATLAB version R19b and newer - 17 Networks




clear; clc

load('C:\Users\Konstantinos Tsirlis\Dropbox\My Desktop\Mourao - Miranda\ABCD results\ABCD Github Repo\warp_atlas_github\schaefer_17_sum.mat')

% Change the path to the framework you have created.
framework_to_run = 'C:\Users\Konstantinos Tsirlis\Dropbox\My Desktop\Mourao - Miranda\ABCD results\abcdAll_noA\framework\spls_fair20-5_fair_modeA_pub_noAnom_new_common_bmi';



% Choose the atlas you want.
cd('C:\Users\Konstantinos Tsirlis\Dropbox\My Desktop\Mourao - Miranda\Atlases_Mikael\warp_atlas\result')
schaefer = spm_vol('Schaefer_Buck_Choi_17Networks_atlas.nii');
schaefer = spm_read_vols(schaefer);
schaefer = reshape(schaefer,1,[]);




schaefer_dummy = schaefer;
for i = 1:234
schaefer_dummy(schaefer_dummy==picu{i,1}) = picu{i,3};
end
tabulate(schaefer_dummy)
schaefer = schaefer_dummy;



% First you have to go inside the framework/model/res folder that you want.
cd([framework_to_run '\res']);

!del summ_table_schaefer_17networks.xlsx



schaefer_areas_ids = (1:17)';


for kk = 1:length(dir('level*'))


    
cd(['level' num2str(kk)])



lala = dir('weightX_*');
lala = lala.name;
weights = spm_vol(lala);
weights = spm_read_vols(weights);
weights = reshape(weights,1,[]);
lala = dir('brain_*');
lala = lala.name;
loadings = spm_vol(lala);
loadings = spm_read_vols(loadings);
loadings = reshape(loadings,1,[]);




% areas of (1) non-zero weights, (2) positive weights, (3) negative weights,
% (4) net weight
for i = 1:length(schaefer_areas_ids)
dummy = schaefer==schaefer_areas_ids(i);
schaefer_weights{i,1} = weights(dummy);
schaefer_weights{i,2} = schaefer_weights{i,1}(:,schaefer_weights{i,1}~=0);
schaefer_weights{i,3} = length(schaefer_weights{i,2});
schaefer_weights{i,4} = schaefer_weights{i,1}(:,schaefer_weights{i,1}>0);
schaefer_weights{i,5} = length(schaefer_weights{i,4});
schaefer_weights{i,6} = schaefer_weights{i,1}(:,schaefer_weights{i,1}<0);
schaefer_weights{i,7} = length(schaefer_weights{i,6});

schaefer_loadings{i,1} = loadings(dummy);
schaefer_loadings{i,2} = schaefer_loadings{i,1}(:,schaefer_loadings{i,1}~=0);
schaefer_loadings{i,3} = length(schaefer_loadings{i,2});
schaefer_loadings{i,4} = schaefer_loadings{i,1}(:,schaefer_loadings{i,1}>0);
schaefer_loadings{i,5} = length(schaefer_loadings{i,4});
schaefer_loadings{i,6} = schaefer_loadings{i,1}(:,schaefer_loadings{i,1}<0);
schaefer_loadings{i,7} = length(schaefer_loadings{i,6});

end


all_non_zero_weights = sum(cell2mat(schaefer_weights(:,3)));
all_non_zero_loadings = sum(cell2mat(schaefer_loadings(:,3)));
sum(weights~=0)


lala = [];
for i = 1:length(schaefer_areas_ids)  
lala = [lala schaefer_weights{i,2}];
end
net_total_weights = sum(abs(lala));
net_positive_weights = sum(lala(lala>0));
flag_pos_weights = 0;
if net_positive_weights==0
    flag_pos_weights = 1;
end
net_negative_weights = sum(lala(lala<0));
flag_neg_weights = 0;
if net_negative_weights==0
    flag_neg_weights = 1;
end

lala = [];
for i = 1:length(schaefer_areas_ids)
lala = [lala schaefer_loadings{i,2}];
end
net_total_loadings = sum(abs(lala));
net_positive_loadings = sum(lala(lala>0));
flag_pos_loadings = 0;
if net_positive_loadings==0
    flag_pos_loadings = 1;
end
net_negative_loadings = sum(lala(lala<0));
flag_neg_loadings = 0;
if net_negative_loadings==0
    flag_neg_loadings = 1;
end




% load('C:\Users\Konstantinos Tsirlis\Dropbox\My Desktop\Mourao - Miranda\Atlases_Mikael\warp_atlas\katia_schaefer.mat')




for i = 1:length(schaefer_areas_ids)
summ_table_weights(i,1) = length(schaefer_weights{i,1});
summ_table_weights(i,2) = 100*schaefer_weights{i,3}/summ_table_weights(i,1); % proportion of the region
summ_table_weights(i,3) = 100*sum(abs(schaefer_weights{i,2}))/net_total_weights;
if flag_pos_weights
    summ_table_weights(i,4) = 0;
else
    summ_table_weights(i,4) = 100*sum(schaefer_weights{i,4})/net_positive_weights;
end
if flag_neg_weights
    summ_table_weights(i,5) = 0;
else
    summ_table_weights(i,5) = 100*sum(schaefer_weights{i,6})/net_negative_weights;
end
end
summ_table_weights(isnan(summ_table_weights(:,2)),2)=0;


for i = 1:length(schaefer_areas_ids)
summ_table_loadings(i,1) = length(schaefer_loadings{i,1});
summ_table_loadings(i,2) = 100*schaefer_loadings{i,3}/summ_table_loadings(i,1); % proportion of the region
summ_table_loadings(i,3) = 100*sum(abs(schaefer_loadings{i,2}))/net_total_loadings;
if flag_pos_loadings
    summ_table_loadings(i,4) = 0;
else
    summ_table_loadings(i,4) = 100*sum(schaefer_loadings{i,4})/net_positive_loadings;
end
if flag_neg_loadings
    summ_table_loadings(i,5) = 0;
else
    summ_table_loadings(i,5) = 100*sum(schaefer_loadings{i,6})/net_negative_loadings;
end
end
summ_table_loadings(isnan(summ_table_loadings(:,2)),2)=0;




%%% THE 2 ADDED COLUMNS
summ_table_weights(:,6) = summ_table_weights(:,4)./summ_table_weights(:,1);
summ_table_weights(:,7) = summ_table_weights(:,5)./summ_table_weights(:,1);
summ_table_loadings(:,6) = summ_table_loadings(:,4)./summ_table_loadings(:,1);
summ_table_loadings(:,7) = summ_table_loadings(:,5)./summ_table_loadings(:,1);



schaefer_areas_names = schaefer17network_names;
picu48 = {'Region name' 'Number of voxels in region' '% of non-zeros weights in the regions' '% of total weights in the regions' '% of total positive weights in the regions' '% of total negative weights in the regions' '% of total positive weights in the regions/no. of voxels' '% of total negative weights in the regions/no. of voxels'};



writecell(picu48,[framework_to_run '\res\summ_table_schaefer_17networks.xlsx'],'Sheet',['Weights_level' num2str(kk)]);
writecell(schaefer_areas_names,[framework_to_run '\res\summ_table_schaefer_17networks.xlsx'],'Range','A2','Sheet',['Weights_level' num2str(kk)]);
writematrix(summ_table_weights,[framework_to_run '\res\summ_table_schaefer_17networks.xlsx'],'Range','B2','Sheet',['Weights_level' num2str(kk)]);

writecell(picu48,[framework_to_run '\res\summ_table_schaefer_17networks.xlsx'],'Sheet',['Loadings_level' num2str(kk)]);
writecell(schaefer_areas_names,[framework_to_run '\res\summ_table_schaefer_17networks.xlsx'],'Range','A2','Sheet',['Loadings_level' num2str(kk)]);
writematrix(summ_table_loadings,[framework_to_run '\res\summ_table_schaefer_17networks.xlsx'],'Range','B2','Sheet',['Loadings_level' num2str(kk)]);


% save(['schaefer_areas_level' num2str(kk)],'schaefer_areas')
clear schaefer_weights summ_table_weights schaefer_loadings summ_table_loadings




cd ..
end