gen_timec_figs.m

%%
% Script for generating timeseries from gamma values (probability of
% occurence of a state) in the format presented in the paper

%%

clear all
close all
clc

%% adding path

path = genpath('/Users/iricchi/Documents/MIP_LabImmersion');
addpath(path)
clc
%% load data

p = '/Users/iricchi/Documents/MIP_LabImmersion/DataFromServer/';

to_load = {'WM_opt_gamma_hats', 'WM_k4_gamma_hats', 'S_opt_gamma_hats', ...
           'S_k3_gamma_hats', 'RM_opt_gamma_hats', 'E_opt_gamma_hats', ...
           'M_opt_gamma_hats', 'M_k5_gamma_hats', 'L_opt_gamma_hats'};
titles = {'WorkingMemory', 'WM k=4', 'Social', 'S k=3', ...
           'RelationalMemory', 'Emotion', 'Motor', 'M k=5', 'Language'};
       
load HCP_Paradigms_updated.mat

Ns = 50; % num subjects

for i = 1:length(to_load)
    ftl = [p,to_load{i}];
    gamma_hats = load(ftl).gamma_hats;
    ks = size(gamma_hats,2);
    val = length(gamma_hats)/Ns;
    
    % compute average gammas from all 50 subjects
    ave_gammas = zeros(val,ks);

    ind = zeros(Ns,val);
    ind(1,1) = 1;
    for j=2:Ns
        ind(j,1) = ind(j-1,1)+val;
    end

    for j = 2:val
       ind(:,j) = ind(:,j-1)+1 ;
    end

    for j = 1:val
        ave_gammas(j,:) = mean(gamma_hats(ind(:,j),:));
    end
    
    % consider the task paradigm 
    if strcmp(titles{i},'Motor') || strcmp(titles{i},'M k=5')
        T = task_paradigm.Motor;
    elseif strcmp(titles{i},'RelationalMemory')
        T = task_paradigm.Relational;
    elseif strcmp(titles{i},'WorkingMemory') || strcmp(titles{i},'WM k=4')
        T = task_paradigm.WM;
    elseif strcmp(titles{i},'Language')
        T = task_paradigm.Language;
    elseif strcmp(titles{i},'Social') || strcmp(titles{i},'S k=3')
        T = task_paradigm.Social;
    elseif strcmp(titles{i},'Emotion')
        T = task_paradigm.Emotion;
    end

    nk = length(unique(T));
 
    T_mat = zeros(length(T), nk);  
                               
    for ii = 1:nk
       T_mat(:,ii) = T == ii-1;
    end
   
    colors = jet(nk);
    y = [0 0 ks ks]; 
    
    % take the indices of tasks to plot squares
    figure()
    for jj = 1:nk
        vec = T_mat(:,jj);
        diffs = diff(vec);
        if vec(1) == 1
            diffs = [1;diffs];
        elseif vec(1) == 0
            diffs = [0;diffs];
        end
        
        lim1 = find(diffs == 1);
        lim2 = find(diffs == -1)-1;
        
        if length(lim1) - length(lim2) == 1
            lim2 = [lim2; length(vec)];
        end
        Y = [];
        X = [];
        for l = 1:length(lim1)
            Y = [Y,y'];
            x = [lim1(l) lim2(l) lim2(l) lim1(l)];
            X = [X, x'];
        end
        
        patch(X,Y,colors(jj,:),'FaceAlpha',0.35, 'EdgeColor', 'none')         
        hold on
    end  
    
    offset = 0;
    for k = 1:ks
        pl = plot(ave_gammas(:,ks+1-k)+offset,'black');
        pl.LineWidth = 2;
        hold on
        offset = offset + 1 ; 
    end
    
    
    yticks(linspace(1-0.4,ks-0.4,ks))
    %inverse plotting
    
    yticks(linspace(1-0.5,ks-0.5,ks))
    yticklabels(linspace(ks,1,ks))
    %yticklabels(linspace(1,ks,ks))
    xlim([0 length(T)])
    xlabel('Timecourse')
    ylabel('K meta-states')
    title(['gamma hats: ', titles{i}])
    set(gca,'FontSize', 15)
    
    
end


%% resting state
ks = 6;
load('RS_k6_gamma_hats')

gammas = gamma_hats(1:500, :);
tmp1 = gammas(:,2);
tmp2 = gammas(:,5);
tmp3 = gammas(:,4);
tmp4 = gammas(:,3);
tmp5 = gammas(:,6);
tmp6 = gammas(:,1);

gammas(:,1) = tmp1;
gammas(:,2) = tmp2;
gammas(:,3) = tmp3;
gammas(:,4) = tmp4;
gammas(:,5) = tmp5;
gammas(:,6) = tmp6;
offset = 0;
for k = 1:ks
    % inverse plotting
    pl = plot(gammas(:,ks+1-k)+offset,'black');
    pl.LineWidth = 1;
    hold on
    offset = offset + 1.1 ; 
end

yticks(linspace(1-0.4,ks-0.4,ks))
%inverse plotting
yticklabels(linspace(ks,1,ks))
xlim([0 length(gammas)])
ylim([0,offset])
xlabel('Timecourse')
ylabel('K meta-states')
title('gamma hats: RS ')
set(gca,'FontSize', 15)