%% Compute TF-Power Memory effect 
%load folder names of subjects
% subjects = [1:22 24:28 30:35 37 39:47 49];
% trialType = 'memory';
% atlasType = 'BA';
% rsawidth = 300;
% rsastep = 50; 
% maxT = 1500;

function subjAvgTFR(trialType,atlasType,atlasIdx,rsawidth,rsastep,maxT,nperm) %trialType can be "memory","category","complexity","memorability"
    
    %set base Dir
    tic
    fs = '\';
    dataDir = '..\Data';
    
    %set colormap Spectral
    cmap = flipud(cbrewer('div','RdBu',51)); cmap(cmap>1) = 1; cmap(cmap<0) = 0;
    
    %% ANALYSIS parameters
    all_subjects = [1:22 24:28 30:37 39:49]; %%%  exclude 23 and 29 for negative d' and 38 for different stimset and 2nd time implant
    load(['..\subject_list.mat']);

    %% RSA Sliding window settings
    %set directory names for storing (ers settings)
    rsa_tag = [num2str(rsawidth) 'width_' num2str(rsastep) 'step_' num2str(maxT) 'ms'];
    rsa_tag = [num2str(rsawidth) 'width_' num2str(rsastep) 'step_-0.5-' num2str(maxT) 'ms']; %tag for tf_data with baseline and extended to 2000ms

    %set directory names for storing (roi selection)
    if strcmp(atlasType,'BA')
        roi_tag = [atlasType char(replace(join(string(atlasIdx)),' ','_'))];
    elseif strcmp(atlasType,'HO')
        %load dir_labels and not numbers for HO naming
        %load(['..\HarvardOxford\labels_both.mat']);
        %fileTag_labels = labels(atlasIdx,3);
        %roi_tag = [atlasType '_' char(replace(join(fileTag_labels),' ','_'))];
    end

    %% find all subjects that contribute to current ROI analysis
    [subjects,numChans,channels] = FindSubjOfROI(all_subjects,atlasType,atlasIdx,'both','all');
    disp([num2str(length(subjects)) ' subjects contribute ' num2str(sum(numChans)) ' channels in ' roi_tag])
    %for all further analyses run from this script, only subjects who
    %contribute channels are included!

    %% compute freq desriptives for each subject
    for num = 1:length(subjects)
        subject = subjects(num);
        fileName = subNames{subject};
        disp(fileName)

        %load RSA prepped TF representations for encoding and
        %recognition
        load([dataDir fs fileName fs 'tf_data' fs 'tfRSA_' rsa_tag '_enc.mat'])
        load([dataDir fs fileName fs 'tf_data' fs 'tfRSA_' rsa_tag '_ret.mat'])
        encoding = encRSA;
        recognition = retRSA; clear encRSA retRSA;

        %get trial indices of trialType
        if strcmp(trialType, 'memory')
                %get remembered and forgotten trials
                enc_type1 = find(encoding.trialinfo(:,2)== 1); %rem
                enc_type2 = find(encoding.trialinfo(:,2)~= 1); %forg

                ret_type1 = find(recognition.trialinfo(:,1) == 1);
                ret_type2 = find(recognition.trialinfo(:,1) == 2);
        end

        %% Load electrode positions and group elecs of same roi
        %get indices of channels in each ROI 
        load([dataDir fs fileName fs 'coordinates' fs 'electrodes_bp_atlas_extended.mat'])

        if strcmp(atlasType,'BA') 
            chanIdx = find(ismember(electrodes.BA,atlasIdx));
            chanLabels = electrodes.BA(chanIdx);
        elseif strcmp(atlasType,'HO')
%             HO labels are ROI names, so read the atlas labels and indices
%             atlasLabels = labels(atlasIdx,2); %get the HO labels of the current ROI selection
%             chanIdx = find(ismember(electrodes.HO,atlasLabels)); %find all electrodes for selection
%             chanLabels = electrodes.HO(chanIdx);
%             fileTag_labels = labels(atlasIdx,3);
        end

        %% compute subj avg across condition (trialType)
        %iterate over ROIs to calc subj avg over roi channels

        %avg over channels 
        cfg = []; 
        cfg.channel = chanIdx;
        cfg.avgoverchan = 'yes';
        encoding_roi = ft_selectdata(cfg,encoding);
        retrieval_roi = ft_selectdata(cfg,recognition);

        % now split by condition
        %calculate the subjAvg TFR for remembered and forgotten per ROI
        %ENCODING
        cfg = [];
        cfg.keeptrials = 'no'; %get average over trials!

        %type1 
        cfg.trials = enc_type1;
        freq_Etype1 = ft_freqdescriptives(cfg, encoding_roi);

        %type2
        cfg.trials = enc_type2;
        freq_Etype2   = ft_freqdescriptives(cfg, encoding_roi);          

        %ret type1
        cfg.trials = ret_type1;
        freq_Rtype1 = ft_freqdescriptives(cfg, retrieval_roi);

        %ret type2
        cfg.trials = ret_type2;
        freq_Rtype2 = ft_freqdescriptives(cfg, retrieval_roi);

        %also compute subject-specific encoding type1-type2 |
        %recognition type1-type2
        freq_sub_E = freq_Etype1;
        freq_sub_E.powspctrm = freq_Etype1.powspctrm - freq_Etype2.powspctrm;

        freq_sub_R = freq_Rtype1;
        freq_sub_R.powspctrm = freq_Rtype1.powspctrm - freq_Rtype2.powspctrm;

        %% save to subj directory
        saveDir = [dataDir fs fileName fs 'stats' fs 'PowerAnalysis' fs rsa_tag fs trialType fs roi_tag];
        status = mkdir(saveDir);

        save([saveDir fs 'freq_Etype1.mat'],'freq_Etype1') %encoding 1
        save([saveDir fs 'freq_Etype2.mat'],'freq_Etype2') %encoding 2
        save([saveDir fs 'freq_Rtype1.mat'],'freq_Rtype1') %retrieval 1
        save([saveDir fs 'freq_Rtype2.mat'],'freq_Rtype2') %retrieval 2
        save([saveDir fs 'freq_Etype12.mat'],'freq_sub_E') %encoding 1-2
        save([saveDir fs 'freq_Rtype12.mat'],'freq_sub_R') %retrieval 1-2
    end
    
    %% compute grad avg over subjects per roi, split by trialType, then clusterstats
        
    %create empty datafile and fill by iterating ov subjects
    subj_data = {};
    info.subjects = {};

    for num = 1:length(subjects)
        subject = subjects(num);
        %check if folder exists
        fileName = subNames{subject};
        info.subjects{end+1,1} = fileName;
        tfDir = [dataDir fs fileName fs 'stats' fs 'PowerAnalysis' fs rsa_tag fs trialType fs roi_tag];
        load([tfDir fs 'freq_Etype1.mat'])
        load([tfDir fs 'freq_Etype2.mat'])
        load([tfDir fs 'freq_Rtype1.mat'])
        load([tfDir fs 'freq_Rtype2.mat'])

        %differences within subject
        load([tfDir fs 'freq_Etype12.mat'])
        load([tfDir fs 'freq_Rtype12.mat'])

        subj_data(end+1,1:6) = {freq_Etype1; freq_Etype2; freq_Rtype1; freq_Rtype2; freq_sub_E; freq_sub_R};
    end
    
    %manully set "chan" label to area name
    for k = 1:size(subj_data,1)
        for m = 1:6
            subj_data{k,m}.label = {roi_tag};
        end
    end

    %% calculate grand average across subjects
    cfg = [];
    cfg.keepindividual = 'yes';
    freq_Etype1 = ft_freqgrandaverage(cfg, subj_data{:,1});
    freq_Etype2 = ft_freqgrandaverage(cfg, subj_data{:,2});
    freq_Rtype1 = ft_freqgrandaverage(cfg, subj_data{:,3});
    freq_Rtype2 = ft_freqgrandaverage(cfg, subj_data{:,4});

    freq_Etype12 = ft_freqgrandaverage(cfg, subj_data{:,5});
    freq_Rtype12 = ft_freqgrandaverage(cfg, subj_data{:,6});

    %avg over subjects
    grand_Etype1 = squeeze(nanmean(freq_Etype1.powspctrm,1));
    grand_Etype2 = squeeze(nanmean(freq_Etype2.powspctrm,1));
    grand_Rtype1 = squeeze(nanmean(freq_Rtype1.powspctrm,1));
    grand_Rtype2 = squeeze(nanmean(freq_Rtype2.powspctrm,1));
    grand_Etype12 = squeeze(nanmean(freq_Etype12.powspctrm,1));
    grand_Rtype12 = squeeze(nanmean(freq_Rtype12.powspctrm,1));

    %compute diff of Remembered - Forgotten
    encoding_grand = grand_Etype1 - grand_Etype2;
    retrieval_grand = grand_Rtype1 - grand_Rtype2;

    %compute SME - Recognition effect (Encoding Memory -
    %Recognition memory effect!)
    encoding_recognition = grand_Etype12 - grand_Rtype12;

    grands = {grand_Etype1; grand_Etype2; encoding_grand; grand_Rtype1; grand_Rtype2; retrieval_grand; encoding_recognition};
    %set titles and labels according to trialtype
    if strcmp(trialType,'memory')
        t1 = 'Rem'; t2 = 'Forg';
        titles = {'Encoding Remembered','Encoding Forgotten','Encoding Rem-Forg','Recognition Remembered','Recognition Forgotten','Recognition Rem-Forg'};
    elseif strcmp(trialType,'category')
        t1 = 'In'; t2 = 'Out';
        titles = {'Encoding Indoor','Encoding Outdoor','Encoding In-Out','Recognition Indoor','Recognition Outdoor','Recognition In-Out'};
    elseif strcmp(trialType,'complexity')
        t1 = 'High'; t2 = 'Low';
        titles = {'Encoding HighComp','Encoding LowComp','Encoding High-Low','Recognition HighComp','Recognition LowComp','Recognition High-Low'};
    elseif strcmp(trialType,'memorability')
        t1 = 'High'; t2 = 'Low';
        titles = {'Encoding HighMem','Encoding LowMem','Encoding High-Low','Recognition HighMem','Recognition LowMem','Recognition High-Low'};
    end
    %plot figure
    figure('Renderer', 'painters', 'Position', [20 20 1400 700]); 
    xlab = freq_Etype1.time*1000 + (rsawidth/2);
    for j = 1:6
        subplot(2,3,j); imagesc(grands{j}); set(gca,'Ydir','normal'); 
        xlabel('time in ms'); xticks(1:size(freq_Etype1.powspctrm,4)); xticklabels(xlab); xtickangle(90);
        ylabel('freq'); yticklabels([10 20 30 80 130]);
        colormap(cmap);  colorbar(); caxis([-0.5 0.5]);
        title(titles{j});
        %plot 0 line
        gcf; xline(0); hold on
    end
    sgtitle([replace(roi_tag,'_',' ') ' | ' num2str(length(info.subjects)) ' subs | ' num2str(sum(numChans)) ' electrodes'])
    
    %% save results to disk
    statsDir = ['..\STATS' fs 'PowerAnalysis_rsafeatures' fs rsa_tag fs trialType fs roi_tag];
    status = mkdir(statsDir);
    
    figName = [roi_tag '_grand_average'];
    saveas(gcf, [statsDir fs figName '.png'])
    %saveas(gcf, [statsDir fs figName '.svg'])

    %save grandaverages
    save([statsDir fs 'ga_freq_E_' t1 '.mat'],'freq_Etype1')
    save([statsDir fs 'ga_freq_E_' t2 '.mat'],'freq_Etype2')
    save([statsDir fs 'ga_freq_R_' t1 '.mat'],'freq_Rtype1')
    save([statsDir fs 'ga_freq_R_' t2 '.mat'],'freq_Rtype2')
    close(gcf);

    %plot and save Encoding-Recognition (Condition) Difference
    figure; imagesc(encoding_recognition); set(gca,'Ydir','normal'); 
    xlabel('time in ms'); xticks(1:size(freq_Etype1.powspctrm,4)); xticklabels(xlab); xtickangle(90);
    ylabel('freq'); yticks([10:10:50]); yticklabels([10 20 30 80 130]);
    colormap(cmap);  colorbar(); caxis([-1 1]);
    title(['Encoding-Recognition X ' trialType])

    saveas(gcf, [statsDir fs 'Encoding-Recognition.png'])
    %saveas(gcf, [statsDir fs 'Encoding-Recognition.svg'])
    close(gcf);

    %% permutation statistics on grand averages across subjects
    cfg = [];
    cfg.channel          = {roi_tag};
    cfg.latency          = [freq_Etype1.time(1) maxT/1000];
    cfg.frequency        = 'all';
    cfg.method           = 'montecarlo';
    cfg.statistic        = 'ft_statfun_depsamplesT';
    cfg.correctm         = 'cluster';
    cfg.clusteralpha     = 0.05;
    cfg.clusterstatistic = 'maxsum';
    %cfg.minnbchan        = 2;
    cfg.tail             = 0;
    cfg.clustertail      = 0;
    cfg.alpha            = 0.025;
    cfg.numrandomization = nperm;

    subj = size(subj_data,1);
    design = zeros(2,2*subj);
    for i = 1:subj
      design(1,i) = i;
    end
    for i = 1:subj
      design(1,subj+i) = i;
    end
    design(2,1:subj)        = 1;
    design(2,subj+1:2*subj) = 2;

    cfg.design   = design;
    cfg.uvar     = 1;
    cfg.ivar     = 2;

    [enc_stat] = ft_freqstatistics(cfg, freq_Etype1, freq_Etype2);
    [rec_stat] = ft_freqstatistics(cfg, freq_Rtype1, freq_Rtype2);

    [enc_rec_stat]= ft_freqstatistics(cfg,freq_Etype12, freq_Rtype12);

    %save stats
    save([statsDir fs 'enc_stats.mat'],'enc_stat')
    save([statsDir fs 'ret_stats.mat'],'rec_stat')
    save([statsDir fs 'enc_rec_stats.mat'], 'enc_rec_stat')
        
    %% plot results of freqstats (t-values!)
    % search for all TF bins with prob < 0.05 and mark them in stats plot
    
    figure('Renderer','painters','Position',[75 200 1200 800]);
    for phase = 1:2 %enc and rec
        switch phase
            case 1
                stat_results = enc_stat;
                title_tag = {'Encoding SME'};
            case 2
                stat_results = rec_stat;
                title_tag = {'Recognition Rem vs Forg'};
        end
    
        %plot stats 
        subplot(1,2,phase);
        imagesc(squeeze(stat_results.stat)); set(gca,'Ydir','normal'); c = colorbar; c.Label.String = 't-value'; colormap(cmap);
        xlabel('time in ms'); xticks(1:size(freq_Etype1.powspctrm,4)); xticklabels(xlab); xtickangle(90);
        ylabel('freq'); yticks([10:10:50]); yticklabels([10 20 30 80 130]); 
        hold on;

        temp = zeros(size(stat_results.prob,2),size(stat_results.prob,3)); %init empty mask
        temp(find(stat_results.prob < cfg.clusteralpha)) = 1;           %set all TF bins 1 which have prob < clusteralpha

        cluster_edges = bwperimtrace(temp);
        for i = 1:1:length(cluster_edges)
            plot(cluster_edges{i}(:,1),cluster_edges{i}(:,2),'Color','black','LineWidth',3)
            hold on
        end
        
        title(title_tag)
    end
    sgtitle([replace(roi_tag,'_',' ') ' | ' num2str(length(subjects)) ' subs | ' num2str(sum(numChans)) ' chans']);
    saveas(gcf,[statsDir fs 'stats.png'])

    %%% for each subj, load avg TF spectrum for remembered and forgotten
    %%% (or other contrasts) trials in encoding and recognition
    
    %%% compute for each timepoint and frequency whether age predicts the
    %%% magnitude of power | also do it over timepoints of sig ERS clusters
    %%% (post-hoc)
    
    %% now compute freqstatistics across age and performance!
    %load subvars again (recognition replaced above)
    load(['..\subject_list.mat']);

    %do it separately for encoding and recognition memory effects -> subj
    %data 5 and 6 (freq_Etype12, freq_Rtype12
    
    cfg = [];
    cfg.channel          = {roi_tag};
    cfg.latency          = [freq_Etype1.time(1) maxT/1000];
    cfg.frequency        = 'all';
    cfg.method           = 'montecarlo';
    cfg.statistic        = 'ft_statfun_correlationT'; %pearson!
    cfg.type             = 'spearman';
    cfg.correctm         = 'cluster';
    cfg.clusteralpha     = 0.05;
    cfg.clusterstatistic = 'maxsum';
    %cfg.minnbchan        = 2;
    cfg.tail             = 0;
    cfg.clustertail      = 0;
    cfg.alpha            = 0.025;
    cfg.numrandomization = nperm;

    %age
    cfg.design = age(subjects);
    cfg.ivar=1;
    
    [enc_age_stat] = ft_freqstatistics(cfg, freq_Etype12);
    [rec_age_stat] = ft_freqstatistics(cfg, freq_Rtype12);
    
    save([statsDir fs 'enc_age_stats.mat'],'enc_age_stat')
    save([statsDir fs 'ret_age_stats.mat'],'rec_age_stat')
    
    %performance
    cfg.design = recognition(subjects);
    cfg.ivar=1;
    
    [enc_perf_stat] = ft_freqstatistics(cfg, freq_Etype12);
    [rec_perf_stat] = ft_freqstatistics(cfg, freq_Rtype12);
    
    save([statsDir fs 'enc_perf_stats.mat'],'enc_perf_stat')
    save([statsDir fs 'ret_perf_stats.mat'],'rec_perf_stat')
        
    %plot results
    %% plot results of freqstats (t-values!)
    % search for all TF bins with prob < 0.05 and mark them in stats plot
    
    for param = 1:2 %age, performance
        figure('Renderer','painters','Position',[75 200 1200 800]);
        for phase = 1:2 %enc and rec
            
            if phase == 1 && param == 1 %enc,age
                    stat_results = enc_age_stat;
                    title_tag = {'Encoding SME * Age'};
                    param_label = 'Age';
            elseif phase == 2 && param == 1
                    stat_results = rec_age_stat;
                    title_tag = {'Recognition Rem vs Forg * Age'};
                    param_label = 'Age';
            elseif phase == 1 && param == 2
                    stat_results = enc_perf_stat;
                    param_label = 'Performance';
                    title_tag = {'Encoding SME * Performance'};
            elseif phase == 2 && param == 2
                    stat_results = rec_perf_stat;
                    title_tag = {'Recognition Rem vs Forg * Performance'};
                    param_label = 'Performance';
            end
            
            %plot stats 
            subplot(1,2,phase);
            imagesc(squeeze(stat_results.stat)); set(gca,'Ydir','normal'); c = colorbar; c.Label.String = 't-value'; colormap(cmap);
            xlabel('time in ms'); xticks(1:size(freq_Etype1.powspctrm,4)); xticklabels(xlab); xtickangle(90);
            ylabel('freq'); yticks([10:10:50]); yticklabels([10 20 30 80 130]); 
            caxis([-3 3]);
            hold on;

            temp = zeros(size(stat_results.prob,2),size(stat_results.prob,3)); %init empty mask
            temp(find(stat_results.prob < cfg.clusteralpha)) = 1;           %set all TF bins 1 which have prob < clusteralpha

            
            cluster_edges = bwperimtrace(temp);
            for i = 1:1:length(cluster_edges)
                plot(cluster_edges{i}(:,1),cluster_edges{i}(:,2),'Color','black','LineWidth',3)
                hold on
            end

            title(title_tag)
        end
        set(gca,'color','white')
        sgtitle([replace(roi_tag,'_',' ') ' | ' upper(trialType) ' X ' param_label ' |  ' num2str(length(subjects)) ' subs | ' num2str(sum(numChans)) ' chans']);
        saveas(gcf,[statsDir fs trialType 'X' param_label '.png'])
    end
end

% 
% %%% plot stats manually
%move to roi_tag dir
%e.g. [..\STATS\PowerAnalysis_rsafeatures\300width_50step_1500ms\memory\BA17_18_19]
% rsa_tag = '300width_50step_1500ms';
% rsa_tag = '300width_50step_-0.5-1500ms';
% roi_tag = 'BA39_40'; direction = 'pos';
% stats_tag = 'item-specific\diff_memcat_row_1000rep'
% cd(['..\STATS\PowerAnalysis_rsafeatures\' rsa_tag '\memory\' roi_tag])
% 
% load('enc_stats.mat')
% load('ret_stats.mat')
% cmap = flipud(cbrewer('div','RdBu',100)); cmap(cmap>1) = 1; cmap(cmap<0) = 0;
% 
% % find timepoints of clusters
% %[clusteridx,enctime,rectime]= extract_ersclusteridx_032122(rsa_tag,stats_tag,'memory',roi_tag,direction);
% 
% %plot vertical bars to indicate clustertime onset and offsets
% enconset = min(enctime); encoffset = max(enctime);
% reconset = min(rectime); recoffset = max(rectime);
% 
% %plot
% freqlabels = {'5','10','15','20','25','30','55','80','105','130'}; freqticks = [5:5:54];
% 
% %if time from 0-1500
% timelabels = {'250','500','750','1000','1250'}; timeticks = [3:5:25];
% %if time from -0.5 to 1500
% timelabels = {'-250','0','250','500','750','1000','1250'}; timeticks = [1:5:35];
% 
% figure('renderer','painters','position',[135 580 550 200]); set(gcf,'color','white');
% subplot(1,2,1);
% imagesc(squeeze(enc_stat.stat(:,:,5:end))); set(gca,'Ydir','normal')
% xlabel('encoding time'); ylabel('freq'); xticks(timeticks); xticklabels(timelabels); xtickangle(90); yticks(freqticks); yticklabels(freqlabels); %ytickangle(90);
% cb = colorbar(); caxis([-3 3]); cb.Title.String = 't-value'; title('encoding sme'); hold on; xline(6);
% %xline(enconset,'--','LineWidth',2);
% %xline(encoffset,'--','LineWidth',2);
% 
% %mark sig cluster
% cluster_edges = bwperimtrace(squeeze(enc_stat.mask(:,:,5:end)));
% for i = 1:1:length(cluster_edges)
%     plot(cluster_edges{i}(:,1),cluster_edges{i}(:,2),'Color','black','LineWidth',3)
% end
% 
% subplot(1,2,2);
% imagesc(squeeze(rec_stat.stat(:,:,5:end))); set(gca,'Ydir','normal')
% xlabel('recognition time'); ylabel('freq'); xticks(timeticks); xticklabels(timelabels); xtickangle(90); yticks(freqticks); yticklabels(freqlabels); %ytickangle(90);
% cb = colorbar(); caxis([-3 3]); cb.Title.String = 't-value'; title('recognition rem vs forg'); hold on; xline(6);
% %xline(reconset,'--','LineWidth',2);
% %xline(recoffset,'--','LineWidth',2);
% 
% %mark sig cluster
% cluster_edges = bwperimtrace(squeeze(squeeze(rec_stat.negclusterslabelmat(:,:,5:end))==1));
% for i = 1:1:length(cluster_edges)
%     plot(cluster_edges{i}(:,1),cluster_edges{i}(:,2),'Color','black','LineWidth',3)
% end
% 
% colormap(cmap)

%%%%%%%%%%%%%%%%%%% 
%% PLOT Occipital Gamma Age effect as TF map & post-hoc correlation with
%%% age
% cmap = flipud(cbrewer('div','RdBu',100)); cmap(cmap>1) = 1; cmap(cmap<0) = 0;
% 
% load(['..\STATS\PowerAnalysis_rsafeatures\300width_50step_-0.5-1500ms\memory\BA17_18_19\enc_age_stats.mat'])
% load(['..\STATS\PowerAnalysis_rsafeatures\300width_50step_-0.5-1500ms\memory\BA17_18_19\enc_perf_stats.mat'])
% 
% %vals = squeeze(enc_age_stat.stat);
% 
% for params = 1:2
%     if params == 1
%         data = enc_age_stat;
%         vals = squeeze(data.stat);
%         val_label = 'age';
%     elseif params == 2
%         data = enc_perf_stat;
%         vals = squeeze(data.stat);
%         val_label = 'performance';
%     end
% 
%     figure;
%     imagesc(vals); set(gca,'Ydir','normal'); colormap(cmap); caxis([-3 3]); hold on;
%     title(val_label)
%     %mark only if sig cluster
%     probs = squeeze(data.prob);
%     cidx = find(probs<0.05);
%     
%     if ~isempty(cidx)
%     
%         map = squeeze(data.posclusterslabelmat);
%         map(map ~= 1) = 0;
% 
%         cluster_edges = bwperimtrace(map);
%         for i = 1:1:length(cluster_edges)
%             plot(cluster_edges{i}(:,1),cluster_edges{i}(:,2),'Color','black','LineWidth',3)
%         end
%         hold on;
%     end
%     
%     set(gcf,'color','white');
%     xline(10);
%     yticks(10:10:50); yticklabels({'10','20','30','80','130'})
%     xticks(5:5:35); xticklabels({'-250','0','250','500','750','1000','1250'});
% 
%     
%     if ~isempty(cidx)
% 
%         cluster_idx = find(map);
% 
%         %%%% get average value from each subject
%         fs = '\';
%         dataDir = '..\DATA';
%         load(['..\subject_list.mat']);
% 
%         [subjects,numChans] = FindSubjOfROI_030822([1:22 24:28 30:37 39:49],'BA',[17 18 19]);
%         disp([num2str(length(subjects)) ' subjects contribute ' num2str(sum(numChans)) ' channels in BA17_18_19'])
% 
%         sub_age = age(subjects);
%         sub_avgPower = []
% 
%         for k = 1:length(subjects)
%             subject = subjects(k);
%             fileName = subNames{subject};
%             disp(fileName)
% 
%             %load mean diff in power between rem and forg for encoding
%             load([dataDir fs fileName fs 'stats' fs 'PowerAnalysis\300width_50step_-0.5-1500ms\memory\BA17_18_19\freq_Etype12.mat'])
%             %get avg across sig groupcluster
%             powspctrm = squeeze(freq_sub_E.powspctrm);
%             sub_avg = nanmean(powspctrm(cluster_idx));
%             sub_avgPower(k) = sub_avg;
%         end
% 
%         plot_correlation_raw(sub_age',sub_avgPower','age','Gamma Power (Rem-Forg)','spearman')
%     end
% end




