addpath(genpath('.'))

CONTRA_OR_IPSI_GRIP   = {'contra', 'ipsi'};

ECOG_SITES = {'contraGrip_ECoG_sites', 'ipsiGrip_ECoG_sites'};

comparisons = {'PLV', 'POW', 'PLV_BEHAV', 'POW_BEHAV'};


for comp = 1:numel(comparisons)
    
    currComp = comparisons{comp};
    if strcmp(currComp, 'PLV_BEHAV') || strcmp(currComp, 'POW_BEHAV') 
        behavVars = {'RT', 'pkFrc', 'pkVel'}; 
    else
        behavVars = {''};    
    end
    
    for b = 1:numel(behavVars)
        behavVar = behavVars{b};
    

        for ipsiGrip = 0:1 % ipsiGrip: 0 = contralateral grips, 1 = ipsilateral grips

            % ECOG_sites: If 1 take the best ECOG contacts found during
            % contralateral gripping, if 2 take the best ECoG contacts during
            % ipsilateral gripping
            for site = 1:2 

                try
                    load(['/matfiles/', currComp, '_', behavVar, CONTRA_OR_IPSI_GRIP{ipsiGrip+1}, 'Grip_', ECOG_SITES{site}], 'data', 'data_perm', 'info')
                catch
                   warning(['File not present: ', '/matfiles/', currComp, '_', behavVar, CONTRA_OR_IPSI_GRIP{ipsiGrip+1}, 'Grip_', ECOG_SITES{site}]) 
                   continue
                end


                fig_h = figure;  
                numSub    = size(data,1);
                % Plot t-scores
                matPlot = squeeze(bsxfun(@rdivide, nanmean(data, 1), (nanstd(data, [], 1) / sqrt(numSub)) ));

                if strcmp(currComp, 'PLV')  
                    % Compare the PLVs against the previously generated
                    % permutation distribution
                    pwDiff             = data - mean(data_perm,4);
                    t_score_rel_toPerm = mean(pwDiff) ./ ( (std(pwDiff)) / sqrt(numSub) );
                    matPlot            = squeeze(t_score_rel_toPerm);

                    if size(data,2) > 1
                        zscores_orig = (squeeze(mean(data)) - mean(squeeze(mean(data_perm,1)),3)) ./ std(squeeze(mean(data_perm,1)),[],3);
                        zscores_perm = bsxfun(@rdivide, bsxfun(@minus, squeeze(mean(data_perm)), mean(squeeze(mean(data_perm)),3)), std(squeeze(mean(data_perm)),[],3));
                    else
                        zscores_orig = (squeeze(mean(data)) - mean(squeeze(mean(data_perm,1)),2)) ./ std(squeeze(mean(data_perm,1)),[],2);
                        zscores_perm = bsxfun(@rdivide, bsxfun(@minus, squeeze(mean(data_perm)), mean(squeeze(mean(data_perm)),2)), std(squeeze(mean(data_perm)),[],2));                    
                    end
                    p_orig = 2 * (1 - normcdf(abs(zscores_orig), 0, 1)); % get p-values from the zscore, abs and *2 to make it 2-tailed
                    p_perm =  2 * (1 - normcdf(abs(zscores_perm), 0, 1)); % get p-values from the zscore, abs and *2 to make it 2-tailed                
                    [~, clusPos_Z_Stat, ~, ~] = getSignifClusters(p_orig, zscores_orig, p_perm, zscores_perm, 0.05, 0.05);
                else
                    % Compare the baseline-normalized power or the differences
                    % in PLV or the differences in power to zero
                    NUM_PERMS = 2000;  % use a high enough number, at least 1000
                    [clusPos_Z_Stat, ~] = permstats(data, zeros(size(data)), NUM_PERMS);                                        
                end


                % To display the representative intervals, stretch the intervals 
                % between events such that they correspond to the average intervals
                % by using interpolation
                INTERPOL_FACTOR = 2;
                matPlot_interpol = interp2(matPlot, INTERPOL_FACTOR);
                len.oldX = size(matPlot, 2);
                len.newX = size(matPlot_interpol, 2);
                len.oldY = size(matPlot, 2);
                len.newY = size(matPlot_interpol, 2);
                centerEvts =  round(info.centerEvts * (len.newX/len.oldX));
                centerEvts(1) = 1;
                centerEvts(5) = len.newX;
                [resizedMat, centerEvts_resized] = stretchEvts(matPlot_interpol, centerEvts);

                imagesc(resizedMat); hold on                       
                cb = colorbar;      
                cb.Label.String = 't-scores';


                switch currComp, 
                    case 'PLV'
                        colormap(Const.greyRed_div_CM)                 
                        caxis([-5 5])    
                    case 'POW'
                        cm=cbrewer('div', 'RdBu', 80); 
                        colormap(flipud(cm))
                        caxis([-8 8])  
                        cb.Label.String = 'Relative change [%]';                          
                    case 'POW_BEHAV'
                        cm=cbrewer('div', 'RdBu', 80); 
                        colormap(flipud(cm))
                        caxis([-4 4])       
                    case 'PLV_BEHAV'
                        cm=cbrewer('div', 'PRGn', 100); 
                        colormap(cm)
                        caxis([-4 4])       
                end

                oldCenterEvts = centerEvts;
                if sum(clusPos_Z_Stat(:))>0 
                    clusPos_Z_Stat_interpol      = interp2(double(clusPos_Z_Stat), INTERPOL_FACTOR);
                    [clusPos_Z_Stat_interpol, ~] = stretchEvts(clusPos_Z_Stat_interpol, oldCenterEvts);
                        
                    % Interpolation yields values between 0-1
                    % but the contour function can only draw lines based on
                    % binary matrices containing only 0s and 1s, hence we
                    % need to cut it dichotomize the interpolation result
                    if strcmp(currComp, 'POW') % if we would have the threshold of 0.4 for the power plots,  
                    % the clusters would look like the positive and negative cluster (in gamma and beta respectively)
                    % were one cluster. To keep separate lines around the
                    % clusters, the threshold for drawing the contour was increased to 0.8
                        clusPos_Z_Stat_interpol = clusPos_Z_Stat_interpol > 0.8;
                    else
                        clusPos_Z_Stat_interpol = clusPos_Z_Stat_interpol > 0.4;  % for the other plots a slightly 
                        % smaller value results in smooth contours
                    end
                   contour(clusPos_Z_Stat_interpol, 1, 'linecolor','k', 'linewidth', 2)
                end
                centerEvts = centerEvts_resized;
                axis xy               

                % Plot the events as vertical dashed lines
                for e = 1:numel(centerEvts)
                    plot([1, 1] * centerEvts(e), [0, 200], '--', 'Color', 'k', 'LineWidth', 2)
                end

                set(gca, 'XTick', centerEvts)
                xLabeling = {'Le/Ri cue', 'GO', '     Grip onset', '& end'};

                set(gca, 'XTickLabel', xLabeling)
                yTicks = 1:2:size(info.freqsOfInt,1);
                yTicks_shifted = (yTicks-0.5)  * (len.newY/len.oldY);
                yTicks_shifted(1) = 1;
                set(gca, 'YTick', yTicks_shifted)
                newLabels    = mean(info.freqsOfInt(yTicks,:),2);
                newLabels(1) = info.freqsOfInt(1,1);
                set(gca, 'YTickLabel', newLabels)

                titling = strrep(['n=', num2str(numel(fileNames)), ' ', CONTRA_OR_IPSI_GRIP{ipsiGrip+1}, ... 
                            ' gripping, ', ECOG_SITES{site}, behavVar], '_', ' ');
                title(titling)

                set(gca, 'FontSize', 16) 

                yPos = 2;
                xPos = 3;
                load(['/matfiles/RiLe_Cue_to_GoSignal_n=', num2str(numel(fileNames)), '.mat'], 'RiLe_Cue_to_GoSignal')
                ms_width = (centerEvts(2)-centerEvts(1)) / (RiLe_Cue_to_GoSignal * 1000);
                plot([xPos, xPos+500*ms_width],[yPos, yPos], 'LineWidth', 3, 'Color', 'k')
                text(xPos, yPos +2.5, '500 ms', 'FontSize', 14)

                % Save the figure
                fig_name = sprintf(['plots/GroupStats_%s_%s_%s_%s'],  ...
                                     currComp, behavVar, CONTRA_OR_IPSI_GRIP{ipsiGrip+1}, ECOG_SITES{site});
                print(fig_h, '-dpng', [fig_name, '.png'], '-r300');

            end
        end
        close all
    end
end
