%Figure 3D: boxplots of change in node strenght relative to baseline

%get data

cd 'D:\Labor\Analyses\Toolbox\fMRI'
addpath('FSLNets')           
addpath('L1precision')      
addpath('pwling')            
addpath('../Tools/BCT')

folder_idx  = pwd; % get folder ID to store data
files_idx   = dir([folder_idx,'/*mat']); % list data2load
baseFileName = files_idx.name;
fullFileName = fullfile(folder_idx, baseFileName);
dirInfos = dir(fullFileName);
clear ('folder_idx','files_idx','baseFileName','fullFileName');

cd   'D:\Labor\Projects\Stroke_tDCS_rsfmri\B_rsfmri_stroke_tDCS_mice\data';
targetFolder=pwd;

clear 'data';
load data.mat;
load labels.mat;

group_name = {'Control','PT','PT_tDCS'};

%% get individual node values

tick  = 1;
for i = 1:size(data,2)
    group = data(i).group;


    for k = 1:size(data(i).subj,2)
        sbj = data(i).subj(k).ID;

        for t = 1:size(data(i).subj(k).tp,2)
            mat = data(i).subj(k).tp(t).mNet;
                if ~isempty(mat)

            for node = 1:length(labels.nr14)
                
                node_strength = sum(mat(labels.nr14(node),labels.nr14),'omitnan');
                master{tick,node+3} = node_strength;
               
            end

                master{tick,1} = group;
                master{tick,2} = sbj;
                master{tick,3} = num2str(t);
                tick = tick+1;

                end
        end
    end
end

header = {'group','subject','timepoint'};

for kx = 1:length(labels.ROI14_short)
    header(kx+3)=labels.ROI14_short(kx);
end

master_tab = cell2table(master,'VariableNames',header);
master_tab.group = categorical(master_tab.group);
master_tab.subject = categorical(master_tab.subject);
master_tab.timepoint = categorical(master_tab.timepoint);
%%
tack = 1;
for ix = 1:length(group_name)
    subdata= master_tab(master_tab.group == group_name{ix},:);
    IDlist = unique(subdata.subject);

    % get group mean baseline values 
    for ox = 1:length(labels.nr14)
        baseValues(ox) = mean(cell2mat(table2cell(subdata(subdata.timepoint=='1',ox+3))));
    end

    % loop over single subjects
    for mx = 1:length(IDlist)
        %determine if subject data at timepoint is available
        if sum(subdata.subject==IDlist(mx) & subdata.timepoint=='3')~=0
            IDdata = cell2mat(table2cell(subdata(subdata.subject==IDlist(mx) & subdata.timepoint=='3',4:end)));
            
            % determine if subject has individual baseline
            if sum(subdata.subject==IDlist(mx) & subdata.timepoint=='1')~= 0
                baseData = cell2mat(table2cell(subdata(subdata.subject==IDlist(mx) & subdata.timepoint=='1',4:end)));
                deltaValues = IDdata - baseData;
            else
                deltaValues = IDdata - baseValues;
            end

            Delta{tack,1} = group_name{ix};
            Delta{tack,2} = string(IDlist(mx));
            for lx = 1:length(labels.ROI14_short)
                Delta{tack,lx+3} = deltaValues(lx);
            end
            tack = tack+1;
       end
    end
end

%% get control shift

    for ox = 1:length(labels.nr14)
        
        baseValues(1,ox) = mean(cell2mat(table2cell(...
            master_tab(master_tab.group == 'Control' & ...
            master_tab.timepoint =='1',ox+3)...
            )));

        baseValues(2,ox) = mean(cell2mat(table2cell(...
            master_tab(master_tab.group == 'Control' & ...
            master_tab.timepoint =='3',ox+3)...
            )));
    end
Delta_base = baseValues(2,:)-baseValues(1,:);

% Delta(:,4:end) = cell2mat(Delta(:,4:end)) - Delta_base;
%%

% get concanated list of data
for kx = 1:length(labels.ROI14_short)
    subseg = Delta(:,[1:3 kx+3]);
    DeltaList([1+(kx-1)*size(subseg,1):(kx*size(subseg,1))],:)= subseg;
    DeltaList([1+(kx-1)*size(subseg,1):(kx*size(subseg,1))],3)= labels.ROI14_short(kx);
end

% transform cell to Table
DeltaHead = {'group','subject','ROI','value'};
Delta_Tab = cell2table(DeltaList,'VariableNames',DeltaHead);
Delta_Tab.ROI = categorical(Delta_Tab.ROI);

bpColors (1,:) = hex2rgb('777777');
bpColors (2,:) = hex2rgb('D95319');
bpColors (3,:) = hex2rgb('0072BD');

%%

PT_Delta = Delta;
GroupID = categorical(PT_Delta(:,1));

rawdata = cell2mat(PT_Delta(GroupID=='PT_tDCS',3:end));

for r = 1:size(rawdata,2)
    pd = fitdist(rawdata(:,r),'Normal');
    ci = paramci(pd);
    statPT(4,r) = pd.mean;
    statPT(5:6,r) = ci(:,1);
end

statPT ([1 4],:) = statPT([1 4],:)- Delta_base;

tabNodePT = array2table(statPT,'VariableNames',labels.ROI14_short)


%%  Permutaion Based Statistics (PBS) of node strength

eval = squeeze(mean(cell2mat(PT_Delta(GroupID=='PT',3:end)))) - ...
    squeeze(mean(cell2mat(PT_Delta(GroupID=='PT_tDCS',3:end)))) ;

% number of permutations
n_permutes = 1000;
conc_data= cell2mat(PT_Delta(:,3:end));

nx=2;%*size(conc_data,2); % Anzahl an tests; Bonferroni Korrektur?
zval_5 = abs(norminv(0.05/nx)); % for p=0.05
zval_1 = abs(norminv(0.01/nx)); % for p=0.01
zval_01 = abs(norminv(0.001/nx)); % for p=0.01
zval=abs(norminv(0.025));

for i=50000:-1:1
    p=i/100000;
    pval(i)=abs(norminv(p));
end

for permi=1:n_permutes
    randorder   = randperm(size(conc_data,1));
    temp        = conc_data(randorder,:);
    perm_map(permi,:) = squeeze(mean(temp(1:floor(size(conc_data,1)/2),:))) - ...
    squeeze(mean(temp(ceil(size(conc_data,1)/2):end,:))); 
end

mean_h0     = squeeze(mean(perm_map,1));
std_h0      = std(perm_map,1);
zstat = abs(eval-mean_h0)./(std_h0);

idx=zstat>zval_5;
imx=zstat>zval;
mat=eval.*idx;

for i=1:size(zstat,1)
    for k=1:size(zstat,2)
        [c idx] = min(abs(pval-zstat(i,k)));
        pstat(i,k)=idx/100000;
    end  
end

FDR=5;
fr=FDR/100;
[hidx, crit_p, adj_ci_cvrg, adj_p]=fdr_bh(pstat,fr,'dep');
hidx

%%

[B,I] = sort(eval,'descend');
ROI_order = labels.ROI14_short(I);

Delta_Tab.ROI = categorical(Delta_Tab.ROI,ROI_order);

figure
b = boxchart(Delta_Tab.ROI,Delta_Tab.value,'GroupByColor',Delta_Tab.group,'BoxWidth',0.75);
for bx=1:3
b(bx).BoxFaceColor = bpColors(bx,:);
b(bx).MarkerColor = bpColors(bx,:);
end

ylabel('Delta Node Strength (Day 14)');
legend({'Control','PT sham','PT tDCS'})

ax = gca;
ax.FontName = 'Times';
ax.FontSize = 22;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Figure 3E-F: Group Mean Boxplots of sensorimotor network
%%
master_tab = DataPTrsfmri;

group_name = {'control','PT sham','PT tDCS'};


%% get control shift

    for ox = 1:length(group_name)
        
        baseValues(1,:) = squeeze(mean(cell2mat(table2cell(...
            master_tab(master_tab.Group == group_name(ox) & ...
            master_tab.Time =='baseline',3:end)...
            )),1,'omitnan'));

        baseValues(2,:) = squeeze(mean(cell2mat(table2cell(...
            master_tab(master_tab.Group == group_name(ox) & ...
            master_tab.Time =='Day 28',3:end)...
            )),1,'omitnan'));

        Delta_base(ox,:) = baseValues(2,:)-baseValues(1,:);
    end


% Delta(:,4:end) = cell2mat(Delta(:,4:end)) - Delta_base;

%%
rawdata = cell2mat(table2cell(master_tab(master_tab.Group=='control'...
    & master_tab.Time=="Day 28",3:end)));

for r = 1:size(rawdata,2)
    pd = fitdist(rawdata(:,r),'Normal');
    ci = paramci(pd);
    stat(1,r) = pd.mean;
    stat([2 3],r) = ci(:,1);
end

% statPT ([1 4],:) = statPT([1 4],:)- Delta_base;
% 
% tabNodePT = array2table(statPT,'VariableNames',labels.ROI14_short)


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Figure 3A-D: Group Mean Boxplots of sensorimotor network

%get data

cd 'D:\Labor\Analyses\Toolbox\fMRI'
addpath('FSLNets')           
addpath('L1precision')      
addpath('pwling')            
addpath('../Tools/BCT')

folder_idx  = pwd; % get folder ID to store data
files_idx   = dir([folder_idx,'/*mat']); % list data2load
baseFileName = files_idx.name;
fullFileName = fullfile(folder_idx, baseFileName);
dirInfos = dir(fullFileName);
clear ('folder_idx','files_idx','baseFileName','fullFileName');

cd   'D:\Labor\Projects\Stroke_tDCS_rsfmri\B_rsfmri_stroke_tDCS_mice\data';
targetFolder=pwd;

clear 'data';
load data.mat;
load labels.mat;

group_name = {'Control','PT','PT_tDCS'};

%%
for ri = 1:size(data,2)
   for li = 1:size(data(ri).subj,2)
        for ti = 1:size(data(ri).subj(li).tp,2)
            
                mnx=data(ri).subj(li).tp(ti).mNet;
                if ~isempty(mnx)
                mNet(:,:,(ri-1)*4+ti,li)=mnx;
                end
            
        end
%         rnx=data(ri).subj(li).mNet;
%         ROI(:,:,li,ri)=rnx;
   end
end

mNet(mNet==0)=nan;
clear ('mnx','rnx','li','ri');
mNet_mean = squeeze(mean(mNet,4,'omitnan'));

for i = 1:3
    for k=1:4
DeltaBaseline (:,:,k+(i-1)*4) = mNet_mean(:,:,k+(i-1)*4)-mNet_mean(:,:,1+(i-1)*4);
    end
end

ROI_Delta = DeltaBaseline(labels.nr14,labels.nr14,:);
%%
for i = 1:size(ROI_Delta,3)
        plotmat_values_rsfMRI(ROI_Delta(:,:,i),targetFolder,dirInfos,[-.75 .75],0,labels.ROI14_short);
        hold on
        axis square;
        plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
        plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
        title(num2str(i))
% title(strcat(name))%,'; mean node strength = ', num2str(norm_ROI_delta(ix*4+k)))); 
        ax=gca;
        ax.FontSize=20;
        ax.FontName='Times';
end