NHP-pRF/Analysis_code/PostFit/matlab/ck_StatsAndPlots.m

% ck_StatsAndPlots
% Takes the pre-processed fitting result tables 
% Calculates stats and creates figures

SaveFigs=false; % autosave generated figs
CloseFigs=false; % autoclose generated figs


%% ========================================================================
%  INITIATE ---------------------------------------------------------------
%  ========================================================================

%% Paths ==================================================================
fprintf('Setting up paths\n')
BaseFld = pwd;

cd ../../../
SHARED_ROOT_FLD = pwd;
cd(BaseFld)

ResFld = fullfile(SHARED_ROOT_FLD,'FitResults','MULTIMODAL','cv1');

% add toolboxes (colorbrewer,matplotlib,boundedline)
addpath(genpath(fullfile(SHARED_ROOT_FLD,'Toolboxes')));

def_cmap = 'Spectral';
ResType = 'mean';
TT = ['Tables_' ResType];

% figure saving folder
pngfld = fullfile(pwd,'fig_png'); [~,~] = mkdir(pngfld);
svgfld = fullfile(pwd,'fig_svg'); [~,~] = mkdir(svgfld);

%% Load ===================================================================
fprintf(['Loading results table. '...
    'Please be patient, this will take a while..\n']);
tic; load(fullfile(ResFld,TT)); t_load=toc;
fprintf(['Loading took ' num2str(t_load) ' s\n']);

%% Correct for result type ================================================
fprintf('Generalizing variable names...\n');
fprintf('MRI...'); 
eval(['tMRI = tMRI_' ResType '; clear tMRI_' ResType ';']);
fprintf('MUA...'); 
eval(['tMUA = tMUA_' ResType '; clear tMUA_' ResType ';']);
fprintf('LFP...'); 
eval(['tLFP = tLFP_' ResType '; clear tLFP_' ResType ';']);
fprintf('>> DONE\n');

%% Split MRI table by model ===============================================
fprintf('Splitting MRI table by model\n')
m = unique(tMRI.Model);
for mi = 1:length(m)
    M = tMRI(strcmp(tMRI.Model,m{mi}),:);
    T(mi).mod = M;
    T(mi).name = m{mi};
    modidx.(m{mi}) = mi;
end

%% Initiate some information ==============================================
fprintf('Generating ROI, model, and subject labels\n')
% proces ROI info
rois = {...
    'V1',   [34];...            % occipital / visual
    'V2',   [131];...           % occipital / visual
    'V3',   [60,93];...     	% occipital / visual
    'V3A'	[123];...			% occipital / visual
    'V4',   [20,39,75];...      % mid-visual
    'V6',   [73];...            % mid-visual
    'V6A',  [141,56];...
    'MT',   [95];...            % mid-visual
    'MST',  [99];...            % mid-visual
    'TEO',  [125];...           % temporal
    'TAa',  [152];...           % temporal
    'Tpt',  [97];...            % temporal (temporal/parietal)
    'TPO',  [159];...           % temporal
    'FST',  [53];...            % temporal
    'A1',	[80];...			% temporal
    'ML'	[25];...			% temporal	
    'AL',	[46];...			% temporal
    'PULV', [197,198];...       % subcortical
    'LGN',  [200,201];...       % subcortical
    'STR',  [175];...           % subcortical
    'MIP',  [89];...  
    'PIP',  [90];...
    'LIP',  [31,130];...        % parietal
    'VIP',  [30];...            % parietal
    '5',    [134];...         	% parietal
    '7',    [91,121];...        % parietal
    'SI',   [50,137];...        % Prim Somatosensory
    'SII',  [63];...            % Sec Somatosensory
    'F2',   [153];...           % premotor
    'F4',   [146];...           % premotor
    'F5',   [129];...           % premotor
    'F7',   [126];...           % premotor
    '8',    [32,51,57,148];...  % frontal FEF  8A: 51, 148; 8B: 32, 57 % combine these
    'CINp', [6,17,27];...       % posterior cingulate
    'CINa', [45,98];...         % anterior cingulate
    'OFC',  [55,107];...        % frontal
    'INS',  [18,87,128];...     % frontal
    'DLPFC',[47,76,127];...     % frontal
    'VMPFC',[5];...         	% frontal
};

roi=rois(:,2);
roilabels=rois(:,1);

MRI_MODELS={...
    'linhrf_cv1_mhrf','linhrf_cv1_dhrf';...
    'linhrf_cv1_mhrf_neggain','linhrf_cv1_dhrf_neggain';...
    'csshrf_cv1_mhrf','csshrf_cv1_dhrf';...
    'doghrf_cv1_mhrf','doghrf_cv1_dhrf';...
    };
ephys_MOD={'linear_ephys_cv1','linear_ephys_cv1_neggain',...
    'css_ephys_cv1','dog_ephys_cv1'};
MMS={...
    'LIN_m','LIN_d';...
    'LIN-N_m','LIN-N_d';...
    'CSS_m','CSS_d';...
    'DOG_m','DOG_d';...
    };

cSUBS = unique(tMRI.Monkey);
SUBS = cellstr(cSUBS);



%% ========================================================================
%  FMRI PRF ANALYSIS ------------------------------------------------------
%  ========================================================================

%% Fig 2C: # significant voxels per ROI (split by monkey and model) =======
RTHRES = 5; % R2 inclusion threshold

ff = figure;
set(ff,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
set(ff,'Position',[10 10 1600 1000]);

paneln=1; 
for s = 1:length(SUBS)
    fprintf(['SUB: ' SUBS{s} '\n'])
    for rowmod=1:4
            subplot(4,2,paneln); hold on;
            nvox=[];
            for r=1:length(roi)
                nvox = [nvox; sum(...
                    T(modidx.(MRI_MODELS{rowmod,1})).mod.R2 > RTHRES & ...
                    strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey, cSUBS(s)) & ...
                    ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                        ck_GetROIidx(roilabels(r),rois) ) ) ...
                     sum(strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey, cSUBS(s)) & ...
                    ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                        ck_GetROIidx(roilabels(r),rois) ) ) ];
            end
            for xval=1:size(nvox,1)
                bar(xval,nvox(xval,1));
            end
            title(['SUB ' SUBS{s} ', ' MMS{rowmod,1} ' nVox with R2 > ' ...
                num2str(RTHRES)],'interpreter','none');
            xlabel('ROI'); ylabel('nVox','interpreter','none');
            set(gca, 'Box','off','yscale','log');
            set(gca,'xtick',1:length(nvox),'xticklabels',roilabels,'TickDir','out');
            xtickangle(45)
            paneln=paneln+1;
    end
end
if SaveFigs
    saveas(ff,fullfile(pngfld, 'MRI_nVoxSign_ROI.png'));
    saveas(ff,fullfile(svgfld, 'MRI_nVoxSign_ROI.svg'));
end
if CloseFigs; close(ff); end

%% SFig 1: prop significant voxels per ROI (split by monkey and model) ====
ff2 = figure;
set(ff2,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
set(ff2,'Position',[10 10 1600 1000]);

paneln=1; 
for s = 1: length(SUBS)
    fprintf(['SUB: ' SUBS{s} '\n'])
    for rowmod=1:4
            subplot(4,2,paneln); hold on;
            nvox=[];
            for r=1:length(roi)
                nvox = [nvox; sum(...
                    T(modidx.(MRI_MODELS{rowmod,1})).mod.R2 > RTHRES & ...
                    strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey, cSUBS(s)) & ...
                    ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                        ck_GetROIidx(roilabels(r),rois) ) ) ...
                     sum(strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey, cSUBS(s)) & ...
                    ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                        ck_GetROIidx(roilabels(r),rois) ) ) ];
            end
            for xval=1:size(nvox,1)
                bar(xval,nvox(xval,1)./nvox(xval,2));
            end
            title(['SUB ' SUBS{s} ', ' MMS{rowmod,1} ' proportion Vox with R2 > ' ...
                num2str(RTHRES)],'interpreter','none');
            xlabel('ROI'); ylabel('Proportion Vox','interpreter','none');
            set(gca, 'Box','off','ylim',[0 .6]);
            set(gca,'xtick',1:length(nvox),'xticklabels',roilabels,'TickDir','out');
            xtickangle(45)
            paneln=paneln+1;
    end
end

if SaveFigs 
    saveas(ff2,fullfile(pngfld, 'MRI_propVoxSign_ROI.png')); 
    saveas(ff2,fullfile(svgfld, 'MRI_propVoxSign_ROI.svg')); 
end
if CloseFigs; close(ff2); end

%% MRI scatter plots & differences R2 per ROI =============================
RTHRES = 0; % include all voxels
for sidx = 0:2 % both monkeys and individuals   
    % scatter plots ---
    % Different ROIS in different colors
    f=figure;
    set(f,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
    set(f,'Position',[10 10 1300 600]);
    
    % Tell us what's happening
    if sidx
        fprintf(['Monkey ' SUBS{sidx} '\n']);
        marker = [' ' SUBS{sidx}];
    else
        fprintf(['Both monkeys\n']);
        marker = ' both';
    end
    
    % plot scatter per area
    figure(f); paneln=1;
    for rowmod=1:4
        for colmod=rowmod+1:4
            subplot(2,3,paneln); hold on;
            PerROI(sidx+1,paneln).Models = {MRI_MODELS{rowmod,1},MRI_MODELS{colmod,1}};
            PerROI(sidx+1,paneln).mR2 = []; PerROI(sidx+1,paneln).seR2 = [];
            for r=1:length(roi)
                s_R2 = T(modidx.(MRI_MODELS{rowmod,1})).mod.R2 >= RTHRES  | ...
                    T(modidx.(MRI_MODELS{colmod,1})).mod.R2 >= RTHRES;
                if sidx
                    s_Monkey = strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey,cSUBS(sidx));
                    SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                        ck_GetROIidx(roilabels(r),rois)) & s_Monkey;
                else
                    SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                        ck_GetROIidx(roilabels(r),rois));
                end
                scatter(...
                    T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS),...
                    T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS),...
                    100,'Marker','.');
                PerROI(sidx+1,paneln).roi(r).sub = marker(2:end);
                PerROI(sidx+1,paneln).roi(r).name = roilabels(r);
                PerROI(sidx+1,paneln).roi(r).R2 = [...
                    T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS),...
                    T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)];
                PerROI(sidx+1,paneln).mR2 = [PerROI(sidx+1,paneln).mR2; ...
                    nanmean(T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS)),...
                    nanmean(T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS))];
                PerROI(sidx+1,paneln).seR2 = [PerROI(sidx+1,paneln).seR2; ...
                    nanstd(T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS))./...
                    sqrt(length(T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS))),...
                    nanstd(T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS))./...
                    sqrt(length(T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)))];
                % stats per ROI
                [PerROI(sidx+1,paneln).roi(r).p,...
                    PerROI(sidx+1,paneln).roi(r).h,...
                    PerROI(sidx+1,paneln).roi(r).stats] = signrank(...
                    PerROI(sidx+1,paneln).roi(r).R2(:,1),...
                    PerROI(sidx+1,paneln).roi(r).R2(:,2));
                PerROI(sidx+1,paneln).roi(r).dir = ...
                    diff(PerROI(sidx+1,paneln).mR2(r,:))./...
                    abs(diff(PerROI(sidx+1,paneln).mR2(r,:)));
            end
            plot([-2 100],[-2 100],'k','Linewidth',1);
            title([MMS{rowmod,1} ' vs ' MMS{colmod,1}],'interpreter','none');
            xlabel(MMS{rowmod,1},'interpreter','none');
            ylabel(MMS{colmod,1},'interpreter','none');
            set(gca, 'Box','off', 'xlim', [-2 100], 'ylim',[-2 100]);
            paneln=paneln+1;
        end
    end
    sgtitle(['SUBJECT' marker])
    if SaveFigs; saveas(f,fullfile(pngfld, ...
            ['MRI_ModelComparison_ROI_R2' marker '.png'])); end
    if SaveFigs; saveas(f,fullfile(svgfld, ...
            ['MRI_ModelComparison_ROI_R2' marker '.svg'])); end
    if CloseFigs; close(f); end
end

% Report stats per ROI ----
fprintf('Compare CSS against P-LIN per ROI -----\n');
nBetter = sum([PerROI(1,2).roi(:).p] < 0.05 & [PerROI(1,2).roi(:).dir] > 0);
nWorse = sum([PerROI(1,2).roi(:).p] < 0.05 & [PerROI(1,2).roi(:).dir] < 0);
nSimilar = sum([PerROI(1,2).roi(:).p] > 0.05);
nROI = length(PerROI(1,2).roi);
fprintf('- Both together -\n');
fprintf([num2str(nBetter) '/' num2str(nROI) ' CSS better\n']);
fprintf([num2str(nWorse) '/' num2str(nROI) ' CSS worse\n']);
fprintf([num2str(nSimilar) '/' num2str(nROI) ' no difference\n']);

nBetter = sum([PerROI(2,2).roi(:).p] < 0.05 & [PerROI(2,2).roi(:).dir] > 0);
nWorse = sum([PerROI(2,2).roi(:).p] < 0.05 & [PerROI(2,2).roi(:).dir] < 0);
nSimilar = sum([PerROI(2,2).roi(:).p] > 0.05);
nROI = length(PerROI(2,2).roi);
fprintf(['- ' SUBS{1} ' -\n']);
fprintf([num2str(nBetter) '/' num2str(nROI) ' CSS better\n']);
fprintf([num2str(nWorse) '/' num2str(nROI) ' CSS worse\n']);
fprintf([num2str(nSimilar) '/' num2str(nROI) ' no difference\n']);

nBetter = sum([PerROI(3,2).roi(:).p] < 0.05 & [PerROI(3,2).roi(:).dir] > 0);
nWorse = sum([PerROI(3,2).roi(:).p] < 0.05 & [PerROI(3,2).roi(:).dir] < 0);
nSimilar = sum([PerROI(3,2).roi(:).p] > 0.05);
nROI = length(PerROI(3,2).roi);
fprintf(['- ' SUBS{2} ' -\n']);
fprintf([num2str(nBetter) '/' num2str(nROI) ' CSS better\n']);
fprintf([num2str(nWorse) '/' num2str(nROI) ' CSS worse\n']);
fprintf([num2str(nSimilar) '/' num2str(nROI) ' no difference\n']);

%% Fig 4A / SFig 3C: MRI scatter plots & differences R2 all voxels ========
RTHRES = 0; % include all voxels
for sidx = 0:2 % both monkeys and individuals   

    % All ROIS together in black
    fsc=figure;
    set(fsc,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
    set(fsc,'Position',[10 10 1300 600]);
    
    % Tell us what's happening
    if sidx
        fprintf(['Monkey ' SUBS{sidx} '\n']);
        marker = [' ' SUBS{sidx}];
    else
        fprintf(['Both monkeys\n']);
        marker = ' both';
    end
    
    % plot scatter over all voxels
    figure(fsc); paneln=1;
    for rowmod=1:4
        for colmod=rowmod+1:4
            subplot(2,3,paneln); hold on;
            XY=[];
            for r=1:length(roi) 
                s_R2 = T(modidx.(MRI_MODELS{rowmod,1})).mod.R2 >= RTHRES  | ...
                    T(modidx.(MRI_MODELS{colmod,1})).mod.R2 >= RTHRES;
                if sidx
                    s_Monkey = strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey,cSUBS(sidx));
                    SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                        ck_GetROIidx(roilabels(r),rois)) & s_Monkey;
                else
                    SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                        ck_GetROIidx(roilabels(r),rois));
                end
                XY=[XY; [T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS) ...
                    T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)] ];
            end
            binscatter(XY(:,1),XY(:,2),100); colorbar; 
            set(gca,'ColorScale','log');
            colormap(inferno)
            plot([-2 100],[-2 100],'k','Linewidth',1);
            title([MMS{rowmod,1} ' vs ' MMS{colmod,1}],'interpreter','none');
            xlabel(MMS{rowmod,1},'interpreter','none');
            ylabel(MMS{colmod,1},'interpreter','none');
            set(gca, 'Box','off', 'xlim', [-2 100], 'ylim',[-2 100]);
            caxis([1 1e4])  
            paneln=paneln+1;
        end
    end
    sgtitle(['SUBJECT' marker])
    if SaveFigs
        saveas(fsc,fullfile(pngfld,['MRI_ModelComparison_R2' marker '.png']));
        saveas(fsc,fullfile(svgfld,['MRI_ModelComparison_R2' marker '.svg'])); 
    end
    if CloseFigs; close(fsc); end
end

% stats over both subs
MR2 = [T(2).mod.R2 T(4).mod.R2 T(7).mod.R2 T(8).mod.R2];
sel=logical(sum(MR2>0,2));
[p,tbl,stats] = kruskalwallis(MR2(sel,:),{'css','dog','p-lin','u-lin'});
fprintf(['Kruskal Wallis ----- p = ' num2str(p) '\n']);
[c,m,h,gnames] = multcompare(stats);
fprintf('Tukey HSD -----\n')
for i=1:size(c,1)
    fprintf([gnames{c(i,1)} ' vs ' gnames{c(i,2)} ...
        ', p = ' num2str(c(i,6))  '\n'])
end

%% SFig 3A-B Difference in R2 across models per ROI =======================
for sidx = 0:2 % both monkeys and individuals
    % Tell us what's happening
    if sidx
        fprintf(['Monkey ' SUBS{sidx} '\n']);
        marker = [' ' SUBS{sidx}];
    else
        fprintf(['Both monkeys\n']);
        marker = ' both';
    end
    
    idx=1;
    for rowmod=1:4
        for colmod=rowmod+1:4
            diffmat{sidx+1,idx}=[];
            for r=1:length(roi)               
                if sidx
                    s_Monkey = strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey,cSUBS(sidx));
                    SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                        ck_GetROIidx(roilabels(r),rois)) & s_Monkey;
                else
                    SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                        ck_GetROIidx(roilabels(r),rois));
                end
                [n,x] = hist(...
                    T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)-...
                    T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS), 100);
                f = n./sum(SSS);
                
                m = mean(T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)-...
                    T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS));
                sd = std(T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)-...
                    T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS));
                se = sd ./ sqrt(sum(SSS));
                diffmat{sidx+1,idx} = [diffmat{sidx+1,idx}; m sd se];
            end
            idx=idx+1;
        end
    end
end
    
f2=figure;
set(f2,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
set(f2,'Position',[100 100 1800 1000]);

cc=1;
for rowmod=1:4
    for colmod=rowmod+1:4
        subplot(3,2,cc); hold on;
        for xval=1:length(diffmat{1,cc})
            bar(xval,diffmat{1,cc}(xval,1));
        end
        for xval=1:length(diffmat{1,cc})
            errorbar(xval,diffmat{1,cc}(xval,1),diffmat{cc}(xval,3),...
                'k-','Linestyle','none');
        end
        set(gca,'xtick',1:length(diffmat{cc}),...
            'xticklabels',roilabels,'ylim',[-2 3],'TickDir','out');
        xlabel('ROI'); ylabel('Diff R2');
        title([MMS{colmod,1} ' - ' MMS{rowmod,1}],...
            'interpreter','none');
        xtickangle(45)
        cc=cc+1;
    end
end
if SaveFigs
    saveas(f2,fullfile(pngfld, 'MRI_ModelComparison_ROI_R2.png'));
    saveas(f2,fullfile(svgfld, 'MRI_ModelComparison_ROI_R2.svg')); 
end
if CloseFigs; close(f2); end

%% only CSS and DoG vs P-LIN per monkey ===================================
f3=figure;
set(f3,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
set(f3,'Position',[100 100 1200 1000]);

subplot(2,2,1); hold on;
for xval=1:length(diffmat{2,2})
    bar(xval,diffmat{2,2}(xval,1));
end
for xval=1:length(diffmat{2,2})
    errorbar(xval,diffmat{2,2}(xval,1),diffmat{2,2}(xval,3),...
        'k-','Linestyle','none');
end
set(gca,'xtick',1:length(diffmat{2,2}),...
    'xticklabels',roilabels,'ylim',[-0.1 2.5],'TickDir','out');
xlabel('ROI'); ylabel('Diff R2');
title('M1: CSS - P-LIN','interpreter','none');
xtickangle(45)

subplot(2,2,2); hold on;
for xval=1:length(diffmat{2,3})
    bar(xval,diffmat{2,3}(xval,1));
end
for xval=1:length(diffmat{2,3})
    errorbar(xval,diffmat{2,3}(xval,1),diffmat{2,3}(xval,3),...
        'k-','Linestyle','none');
end
set(gca,'xtick',1:length(diffmat{2,3}),...
    'xticklabels',roilabels,'ylim',[-0.1 2.5],'TickDir','out');
xlabel('ROI'); ylabel('Diff R2');
title('M1: DoG - P-LIN','interpreter','none');
xtickangle(45)

subplot(2,2,3); hold on;
for xval=1:length(diffmat{3,2})
    bar(xval,diffmat{3,2}(xval,1));
end
for xval=1:length(diffmat{3,2})
    errorbar(xval,diffmat{3,2}(xval,1),diffmat{3,2}(xval,3),...
        'k-','Linestyle','none');
end
set(gca,'xtick',1:length(diffmat{3,2}),...
    'xticklabels',roilabels,'ylim',[-0.1 2.5],'TickDir','out');
xlabel('ROI'); ylabel('Diff R2');
title('M2: CSS - P-LIN','interpreter','none');
xtickangle(45)

subplot(2,2,4); hold on;
for xval=1:length(diffmat{3,3})
    bar(xval,diffmat{3,3}(xval,1));
end
for xval=1:length(diffmat{3,3})
    errorbar(xval,diffmat{3,3}(xval,1),diffmat{3,3}(xval,3),...
        'k-','Linestyle','none');
end
set(gca,'xtick',1:length(diffmat{3,3}),...
    'xticklabels',roilabels,'ylim',[-0.1 2.5],'TickDir','out');
xlabel('ROI'); ylabel('Diff R2');
title('M2: DoG - P-LIN','interpreter','none');
xtickangle(45)

if SaveFigs
    saveas(f3,fullfile(pngfld, 'MRI_SelModelComparison_ROI_R2.png')); 
    saveas(f3,fullfile(svgfld, 'MRI_SelModelComparison_ROI_R2.svg')); 
end
if CloseFigs; close(f3); end

%% Fig 4B NegGain Fits: Characterize ======================================
R2th = 5; % minimum R2
R2enh = 5; % R2 improvement

DoG = tMRI(...
    strcmp(tMRI.Model,'doghrf_cv1_mhrf'),:);
lin_n = tMRI(...
    strcmp(tMRI.Model,'linhrf_cv1_mhrf_neggain'),:);
lin = tMRI(...
    strcmp(tMRI.Model,'linhrf_cv1_mhrf'),:);

f_neg2 = figure;
set(f_neg2,'Position',[10 10 1200 1000]);

vox_sel = lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh & ...
    (lin_n.ecc<25 & lin.ecc<25 & DoG.ecc<25) & ....
    (lin_n.rfs<20 & lin.rfs<20 & DoG.rfs<20);

fprintf('-------------\n');

subplot(2,3,1); hold on;
histogram(lin_n.gain(vox_sel),-100:0.1:100,...
    'Normalization','probability','FaceColor','k','FaceAlpha',0.75);
xlabel('gain LIN-POSNEG');ylabel('nvoxels');
set(gca,'xlim',[-2.6 2.6]);
MM=median(lin_n.gain(vox_sel));
yy=get(gca,'ylim');
plot([MM MM], [0 1],'k','Linewidth',2)
set(gca,'ylim',[0 0.3],'TickDir','out');
title('Gain SELECTED voxels')
fprintf(['MEDIAN GAIN SEL: ' num2str(MM) '\n'])

% Wilcoxon 1-tailed < 1
[p,h,stats] = signrank(lin_n.gain(vox_sel),0,'tail','left');
fprintf(['Gain SEL < 0: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

% all vox
subplot(2,3,4); hold on;
histogram(lin_n.gain,-100:0.1:100,...
    'Normalization','probability','FaceColor','k','FaceAlpha',0.75);
xlabel('gain LIN-POSNEG');ylabel('nvoxels');
set(gca,'xlim',[-2.6 2.6]);
MM=median(lin_n.gain);
yy=get(gca,'ylim');
plot([MM MM], [0 1],'k','Linewidth',2)
set(gca,'ylim',[0 0.3],'TickDir','out');
title('Gain ALL voxels')
fprintf(['MEDIAN GAIN ALL: ' num2str(MM) '\n'])

% Wilcoxon 1-tailed < 1
[p,h,stats] = signrank(lin_n.gain,0,'tail','right');
fprintf(['Gain ALL < 0: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);
fprintf('-------------\n');

subplot(2,3,2); hold on;
bb = [lin_n.ecc(vox_sel) lin.ecc(vox_sel)];
plot([1 2],mean(bb),'o','Linewidth',2)
errorbar([1 2],mean(bb),std(bb),'k','Linewidth',2)
plot([1 2],mean(bb),'o','MarkerSize',10,'MarkerFaceColor','k','MarkerEdgeColor','k')
set(gca,'xtick',1:2,'xticklabels',{'LIN-N','LIN'},...
    'ylim',[0 15],'xlim',[0.8 2.2],'TickDir','out')
ylabel('Eccentricity');
title('Ecc')

subplot(2,3,3); hold on;
histogram(lin.ecc(vox_sel)-lin_n.ecc(vox_sel),-10:0.5:10,...
    'FaceColor','k','FaceAlpha',0.5);
xlabel('Ecc. Diff (POS-POSNEG)');ylabel('nvoxels');
MM=median(bb(:,2)-bb(:,1));
yy=get(gca,'ylim');
plot([MM MM], [0 yy(2)+40],'k','Linewidth',5)
set(gca,'ylim',[0 yy(2)+100],'TickDir','out');
title('Ecc Diff')
fprintf(['MEDIAN ECC DIFF: ' num2str(MM) '\n'])

bbecc=bb;

subplot(2,3,5); hold on;
bb = [lin_n.rfs(vox_sel) lin.rfs(vox_sel)];
plot([1 2],mean(bb),'o','Linewidth',2)
errorbar([1 2],mean(bb),std(bb),'k','Linewidth',2)
plot([1 2],mean(bb),'o','MarkerSize',10,'MarkerFaceColor','k','MarkerEdgeColor','k')
set(gca,'xtick',1:2,'xticklabels',{'LIN-N','LIN'},...
    'ylim',[0 5.1],'xlim',[0.8 2.2],'TickDir','out')
ylabel('Size');
title('Size')

subplot(2,3,6); hold on;
histogram(lin.rfs(vox_sel)-lin_n.rfs(vox_sel),-10:0.25:10,...
    'FaceColor','k','FaceAlpha',0.5);
xlabel('Size Diff (POS-POSNEG)');ylabel('nvoxels');
MM=median(bb(:,2)-bb(:,1));
yy=get(gca,'ylim');
plot([MM MM], [0 yy(2)+40],'k','Linewidth',5)
set(gca,'ylim',[0 yy(2)+110]);
set(gca,'xlim',[-5.1 5.1],'TickDir','out');
title('Size Diff')
fprintf(['MEDIAN SZ DIFF: ' num2str(MM) '\n'])

bbsz=bb;
sgtitle('pRFs POS LINEAR vs POSNEG LINEAR model')

if SaveFigs
    saveas(f_neg2,fullfile(pngfld, 'MRI_NEG-PRF2.png'));
    saveas(f_neg2,fullfile(svgfld, 'MRI_NEG-PRF2.svg')); 
end
if CloseFigs; close(f_neg2); end

% Wilcoxon 1-tailed < 0
[p,h,stats] = signrank(bbecc(:,1),bbecc(:,2));
fprintf(['Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

%% DoG Fits: Characterize =================================================
f_neg3 = figure;
set(f_neg3,'Position',[10 10 1300 1100]);
vox_sel = DoG.R2>R2th & DoG.R2>lin.R2+R2enh;

subplot(1,3,1); hold on;
histogram(DoG.normamp(vox_sel),-10:0.1:10,'FaceColor','k','FaceAlpha',0.5);
xlabel('INH nAMP');ylabel('nvoxels');
set(gca,'xlim',[-0.2 3.1]);
MM=median(DoG.normamp(vox_sel));
yy=get(gca,'ylim');
plot([MM MM], [0 yy(2)+40],'k','Linewidth',5)
set(gca,'ylim',[0 yy(2)+30],'TickDir','out');
title('NORMAMP')
fprintf(['MEDIAN NAMP: ' num2str(MM) '\n'])


% Wilcoxon 1-tailed > 1
[p,h,stats] = signrank(DoG.normamp(vox_sel),1,'tail','right');
fprintf(['nAmp > 1: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);
fprintf(['Median nAMP: ' num2str(median(DoG.normamp(vox_sel))) ', IQR: '...
    num2str(iqr(DoG.normamp(vox_sel))) '\n'])

subplot(1,3,2); hold on;
bb = [DoG.ecc(vox_sel) lin.ecc(vox_sel)];
bbecc=bb;
plot([1 2],mean(bb),'o','Linewidth',2)
errorbar([1 2],mean(bb),std(bb),'k','Linewidth',2)
plot([1 2],mean(bb),'o','MarkerSize',10,'MarkerFaceColor','k','MarkerEdgeColor','k')
set(gca,'xtick',1:2,'xticklabels',{'DoG','LIN'},...
    'ylim',[0 20],'xlim',[0.8 2.2],'TickDir','out')
ylabel('Eccentricity');
title('Ecc Diff')

subplot(1,3,3); hold on;
histogram(lin.ecc(vox_sel)-DoG.ecc(vox_sel),-10:0.5:10,...
    'FaceColor','k','FaceAlpha',0.5);
xlabel('Ecc. Diff (POS-DoG)');ylabel('nvoxels');
MM=median(bb(:,2)-bb(:,1));
yy=get(gca,'ylim');
plot([MM MM], [0 yy(2)+40],'k','Linewidth',5)
set(gca,'ylim',[0 yy(2)+110],'TickDir','out');
title('Ecc Diff')
fprintf(['MEDIAN ECC DIFF: ' num2str(MM) '\n'])

sgtitle('pRFs POS LINEAR vs DoG model')

if SaveFigs
    saveas(f_neg3,fullfile(pngfld, 'MRI_NEG-PRF3.png'));
    saveas(f_neg3,fullfile(svgfld, 'MRI_NEG-PRF3.svg'));
end
if CloseFigs; close(f_neg3); end

% Wilcoxon 1-tailed < 1
[p,h,stats] = signrank(bbecc(:,1),bbecc(:,2));
fprintf(['Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

%% Value of exponential parameter for CSS across ROIs =====================
RTHRES = 5;
fexp = figure;
set(fexp,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
set(fexp,'Position',[10 10 1600 300]);
hold on;

for s = 1:length(SUBS)
    fprintf(['SUB: ' SUBS{s} '\n'])
    exptv(s).roi={};
    for r=1:length(roi)
        exptvals = T(modidx.csshrf_cv1_mhrf).mod.expt(...
            T(modidx.csshrf_cv1_mhrf).mod.R2 > RTHRES & ...
            strcmp(T(modidx.csshrf_cv1_mhrf).mod.Monkey,cSUBS(s)) & ...
            ismember(T(modidx.csshrf_cv1_mhrf).mod.ROI,...
            ck_GetROIidx(roilabels(r),rois) ) );
        exptv(s).roi{r}=exptvals;
    end
end
all_expt=[];
for xval=1:length(roi)
    if size([exptv(1).roi{xval};exptv(2).roi{xval}],1) > 4
        bar(xval,mean([exptv(1).roi{xval};exptv(2).roi{xval}]));
        all_expt = [all_expt; exptv(1).roi{xval}; exptv(2).roi{xval}];
    else
        bar(xval,0);
        all_expt = [all_expt; exptv(1).roi{xval}; exptv(2).roi{xval}];
    end
end
for xval=1:length(roi)
    if size([exptv(1).roi{xval};exptv(2).roi{xval}],1) > 4
        errorbar(xval,...
            mean([exptv(1).roi{xval};exptv(2).roi{xval}]),...
            std([exptv(1).roi{xval};exptv(2).roi{xval}]),...
            'k','Linestyle','none');
        %     errorbar(xval,...
        %         mean([exptv(1).roi{xval};exptv(2).roi{xval}]),...
        %         std([exptv(1).roi{xval};exptv(2).roi{xval}])./...
        %         sqrt(length([exptv(1).roi{xval};exptv(2).roi{xval}])),...
        %         'k','Linestyle','none');
    end
end
title(['Avg EXPT parameter, R2 > ' ...
    num2str(RTHRES)],'interpreter','none');
ylabel('EXPT PM','interpreter','none');
set(gca,'xtick',1:length(roi),'xticklabels',roilabels,'TickDir','out');
xtickangle(45)
if SaveFigs
    saveas(fexp,fullfile(pngfld, 'MRI_CSS_EXPT-PM_ROI.png')); 
    saveas(fexp,fullfile(svgfld, 'MRI_CSS_EXPT-PM_ROI.svg')); 
end
if CloseFigs; close(fexp); end

% Wilcoxon 1-tailed < 1
[p,h,stats] = signrank(all_expt,1,'tail','left');
fprintf(['ALL VOXELS - Wilcoxon EXPT < 1 : z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

fprintf('\n------\nSTATS\n------\n')
for xval=1:length(roi)
    if size([exptv(1).roi{xval};exptv(2).roi{xval}],1) > 4
        [p,h,stats] = signrank([exptv(1).roi{xval};exptv(2).roi{xval}],1,'tail','left');
        if isfield(stats,'zval')  
            fprintf([roilabels{xval} ' VOXELS - Wilcoxon EXPT < 1 : z = ' ...
                num2str(stats.zval) ', p = ' num2str(p) '\n']);
        end
    else
        fprintf([roilabels{xval} ': Not enough values for statistics\n'])
    end
end

%% SFig 13: MRI scatter plot HRF & differences ============================
RTHRES = 0;

f3a=figure;
set(f3a,'Position',[100 100 1500 400]);
set(f3a,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
s_R2 = T(modidx.linhrf_cv1_mhrf).mod.R2>RTHRES;

f3b=figure;
set(f3b,'Position',[100 100 1300 1200]);
set(f3b,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
s_R2 = T(modidx.linhrf_cv1_mhrf).mod.R2>RTHRES;

for mm = 1:size(MRI_MODELS,1)
    
    figure(f3a)
    subplot(1,4,mm); hold on;
    XY=[];
    for r=1:length(roi)
        SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                ck_GetROIidx(roilabels(r),rois) );
%         scatter(T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS),...
%             T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS),100,'Marker','.');
        XY=[XY; T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS) ...
            T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS)];
    end
    binscatter(XY(:,1),XY(:,2),100); colorbar; 
    set(gca,'ColorScale','log');
    colormap(inferno)
    caxis([1 1e4]) 
    plot([-2 100],[-2 100],'k','Linewidth',1);   
    title(['mHRF vs cHRF (' MRI_MODELS{mm,1} ')'],'interpreter','none'); 
    xlabel('Monkey HRF R2'); ylabel('Canonical HRF R2');
    set(gca, 'Box','off', 'xlim', [-2 100], 'ylim',[-2 100]);

    diffmat2{1}=[]; dH=[];
    for r=1:length(roi)
        SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
                ck_GetROIidx(roilabels(r),rois) );
        [n,x] = hist(...
            T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS)-...
            T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS),100);
        f = n./sum(SSS);
    
        m = mean(T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS)-...
            T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS));
        sd = std(T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS)-...
            T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS));
        nvox = sum(SSS);
        se = sd ./ sqrt(nvox);
        diffmat2{1} = [diffmat2{1}; m sd se nvox];
        dH=[dH;...
            T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS) - ...
            T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS) ];
    end
    
    figure(f3b)
    subplot(4,1,mm); hold on
    for xval=1:length(diffmat2{1})
        bar(xval,diffmat2{1}(xval,1));
    end
    for xval=1:length(diffmat2{1})
        errorbar(xval,diffmat2{1}(xval,1),diffmat2{1}(xval,3),...
        'k-','Linestyle','none')
    end
    % mean (weighted mean taking nvox into account)
    mAll = sum((diffmat2{1}(:,1).*diffmat2{1}(:,4)))./...
        sum(diffmat2{1}(:,4)); 
    mAll = mean(dH); 
    stdAll = std(dH);
    
    
    % Wilcoxon 1-tailed < 1
    [p,h,stats] = signrank(dH,0);
    fprintf(['ALL VOXELS - Wilcoxon dHRF (' MRI_MODELS{mm,1} ') < 1 : z = ' ...
        num2str(stats.zval) ', p = ' num2str(p) '\n']);

    text(0.5, -1.2,[num2str(mAll) ' +/- ' num2str(stdAll)])
    set(gca,'xticklabels',[],'ylim',[-1.6 1.6],'TickDir','out');
    %xlabel('ROI'); 
    ylabel('Diff R2');
    title(['mHRF - cHRF (' MRI_MODELS{mm,1} ')'],'interpreter','none'); 
    %legend(roilabels,'interpreter','none','Location','NorthEast');
    set(gca,'xtick',1:length(diffmat2{1}),'xticklabels',roilabels);
    xtickangle(45)    
end
if SaveFigs
    saveas(f3a,fullfile(pngfld, 'HRF_Comparison.png')); 
    saveas(f3a,fullfile(svgfld, 'HRF_Comparison.svg')); 
    saveas(f3b,fullfile(pngfld, 'HRF_Comparison_binned.png')); 
    saveas(f3b,fullfile(svgfld, 'HRF_Comparison_binned.svg')); 
end
if CloseFigs; close(f3a); close(f3b); close all; end

%% MRI rf size depending on HRF ===========================================
R2th=5;
s_R2 = T(modidx.csshrf_cv1_mhrf).mod.R2 > R2th;

for mm = 1:size(MRI_MODELS,1)
    xy_sz{mm}=[]; xy_ecc{mm}=[];
    
    for r=1:length(roi)
        SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
            ck_GetROIidx(roilabels(r),rois) );
        xy_sz{mm} = [xy_sz{mm};...
            T(modidx.(MRI_MODELS{mm,1})).mod.rfs(SSS) ...
            T(modidx.(MRI_MODELS{mm,2})).mod.rfs(SSS)];
        xy_ecc{mm} = [xy_ecc{mm};...
            T(modidx.(MRI_MODELS{mm,1})).mod.ecc(SSS) ...
            T(modidx.(MRI_MODELS{mm,2})).mod.ecc(SSS)];   
    end
    xy_sz{mm}( isinf(xy_sz{mm})  ) = nan;
    xy_ecc{mm}( isinf(xy_ecc{mm})  ) = nan;

    % Wilcoxon
    fprintf(['--' MMS{mm,1} '--\n'])
    [p,h,stats] = signrank(xy_sz{mm}(:,1),xy_sz{mm}(:,2),'method','approximate');
    fprintf(['ALL VOXELS R2>TH - Wilcoxon HRF sz: z = ' ...
        num2str(stats.zval) ', p = ' num2str(p) '\n']);
    fprintf(['Mean:' num2str(nanmean(diff(xy_sz{mm},1,2))) ', std ' ...
        num2str(nanstd(diff(xy_sz{mm},1,2))) '\n'])
    
    [p,h,stats] = signrank(xy_ecc{mm}(:,1),xy_ecc{mm}(:,2),'method','approximate');
    fprintf(['ALL VOXELS R2>TH - Wilcoxon HRF ecc: z = ' ...
        num2str(stats.zval) ', p = ' num2str(p) '\n']);
    fprintf(['Mean:' num2str(nanmean(diff(xy_ecc{mm},1,2))) ', std ' ...
        num2str(nanstd(diff(xy_ecc{mm},1,2))) '\n'])
    
    diffmat2{1}=[];
    for r=1:length(roi)
        SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,...
            ck_GetROIidx(roilabels(r),rois) );
        [n,x] = hist(...
            T(modidx.(MRI_MODELS{mm,1})).mod.rfs(SSS) - ...
            T(modidx.(MRI_MODELS{mm,2})).mod.rfs(SSS), 100);
        f = n./sum(SSS);
        
        dsz = T(modidx.(MRI_MODELS{mm,1})).mod.rfs - ...
            T(modidx.(MRI_MODELS{mm,2})).mod.rfs;
        dsz = dsz(SSS);
        dsz = dsz(isfinite(dsz));
        
        m = mean(dsz);
        sd = std(dsz);
        se = sd ./ sqrt(length(dsz));
        diffmat2{1} = [diffmat2{1}; m sd se size(dsz,1)];
    end
  
end

%% MRI ECC vs PRF Size ====================================================
Rth=5;

f5=figure;
set(f5,'Position',[100 100 800 1200]);
set(f5,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));

for m=1:length(MRI_MODELS)-1
    s_R2 = T(modidx.(MRI_MODELS{m,1})).mod.R2 > Rth;
    EccBin = 0.5:1:16.5;
    
    subplot(3,2,(m-1)*2 +1);hold on;
    for r=1:length(roi)
        SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,...
            ck_GetROIidx(roilabels(r),rois) );
        
        ES{r}=[];
        scatter(T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS),...
            T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS),100,'Marker','.');        
        
        for b=1:length(EccBin)
            bb=[EccBin(b)-0.5 EccBin(b)+0.5];
            PSZ=T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS);
            ECC=T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS);
            ES{r}=[ES{r}; EccBin(b) median(PSZ(ECC>=bb(1) & ECC<=bb(2)))];
        end
    end
    title(['Ecc vs pRF size [' MMS{m} ', R>' num2str(Rth) ']'],...
        'interpreter','none');
    xlabel('Eccentricity');ylabel('pRF size');
    set(gca, 'Box','off', 'xlim', [0 10], 'ylim',[0 10]);
    
    subplot(3,2,(m-1)*2 +2);hold on;
    for r=1:length(roi)
        h=plot(ES{r}(:,1),ES{r}(:,2),'o');
        set(h,'MarkerSize',6,'markerfacecolor', get(h, 'color'));
    end
    title(['Ecc vs pRF size [' MMS{m} ', R>' num2str(Rth) ']'],...
        'interpreter','none'); xlabel('Eccentricity');ylabel('pRF size');
    set(gca, 'Box','off', 'xlim', [0 10], 'ylim',[0 10]);
end
if SaveFigs
    saveas(f5,fullfile(pngfld, 'MRI_Ecc_vs_Size.png'));
    saveas(f5,fullfile(svgfld, 'MRI_Ecc_vs_Size.svg'));
end
if CloseFigs; close(f5); end

%% Fig 5 & SFig 5: Ecc-Sz for CSS per area including CI ===================
SaveFigs=false;

for Rth=[5 3]
    
    clear tbl lm;
    
    max_ecc = 100;
    bins=[0:18; 2:20]';
    
    f5css=figure;
    set(f5css,'Position',[100 100 2000 1500]);
    set(f5css,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
    sgtitle(['Voxels & fit, R>' num2str(Rth)])
    
    m=3; % css
    s_R2 = T(modidx.(MRI_MODELS{m,1})).mod.R2 > Rth & ...
        T(modidx.(MRI_MODELS{m,1})).mod.ecc < max_ecc;
    
    nsp=length(roi);
    nrc=ceil(sqrt(nsp));
    
    nvox_roi = []; binned_data=[];
    for r=1:length(roi)
        ES{r} = []; ES_m{r} = [];
        ESb{r} = []; ESf{r} = [];
        
        SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,...
            ck_GetROIidx(roilabels(r),rois) );
        
        s_Monkey = strcmp(T(modidx.(MRI_MODELS{m,1})).mod.Monkey,'M01');
        SSS_m1 = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,...
            ck_GetROIidx(roilabels(r),rois)) & s_Monkey;
        
        s_Monkey = strcmp(T(modidx.(MRI_MODELS{m,1})).mod.Monkey,'M02');
        SSS_m2 = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,...
            ck_GetROIidx(roilabels(r),rois)) & s_Monkey;
        
        nvox_roi = [nvox_roi; sum(SSS_m1) sum(SSS_m2)];
        roi_n{r,1} = rois{r};
        roi_n{r,2} = sum(SSS_m1);
        roi_n{r,3} = sum(SSS_m2);
        
        ES{r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS) ...
            T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS) ];
        ES{r}( ES{r}(:,2) > 50,: ) = []; % remove insanely large pRF's
        
        ES_m{1,r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS_m1) ...
            T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS_m1) ];
        ES_m{2,r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS_m2) ...
            T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS_m2) ];
        
        %subplot(nrc,nrc,r); hold on;
        subplot(5,8,r); hold on;
        
        scatter(ES_m{1,r}(:,1),ES_m{1,r}(:,2),10,'Marker','o',...
            'MarkerFaceColor','b','MarkerFaceAlpha',.2,...
            'MarkerEdgeColor','none');
        scatter(ES_m{2,r}(:,1),ES_m{2,r}(:,2),10,'Marker','o',...
            'MarkerFaceColor','k','MarkerFaceAlpha',.2,...
            'MarkerEdgeColor','none');
        
        % make table
        tbl = table(ES{r}(:,1),ES{r}(:,2),'VariableNames',{'ecc','sz'});
        if ~isempty(tbl)
            % remove inf's
            tbl(isinf(tbl.sz),:)=[];
            % fit
            lm{r} = fitlm(tbl,'sz ~ 1 + ecc');
            % CI
            lmCI{r}=coefCI(lm{r},0.05);
            if ~isnan(lm{r}.Coefficients.pValue(1))
                % Prediction
                x=[0;20];
                [ypred,yci] = predict(lm{r},x);
                yci2=yci(:,2)-ypred;
                if lm{r}.Coefficients.pValue(2) < 0.05
                    boundedline(x,ypred',yci2,'r','alpha')
                else
                    boundedline(x,ypred',yci2,'k','alpha')
                end
                %fprintf([roilabels{r} ' n = ' num2str(size(tbl,1)) '\n'])
                %fprintf([roilabels{r} ' p = ' num2str(lm{r}.Coefficients.pValue(2)) '\n'])
                fprintf([roilabels{r} ' slope = ' num2str(lm{r}.Coefficients.Estimate(2)) '\n'])
            end
            
            % binned
            for b=1:size(bins,1)
                seldata = (tbl.ecc>bins(b,1) & tbl.ecc<=bins(b,2)) ;
                if sum(seldata)>1
                    binned_data=[binned_data;...
                        r b mean(tbl.sz(seldata)) std(tbl.sz(seldata))./sqrt(sum(seldata))];
                end
            end
        end
        
        title(['Ecc vs pRF size [' roilabels{r} ', R>' num2str(Rth) ']'],...
            'interpreter','none');
        title(roilabels{r});
        xlabel('Eccentricity');ylabel('pRF size');
        set(gca, 'Box','off', 'xlim', [0 20], 'ylim',[0 25]);
    end
    
    % binned
    f5css_binned=figure;
    set(f5css_binned,'Position',[100 100 2000 1500]);
    set(f5css_binned,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
    sgtitle(['Binned & voxel-based fit, R>' num2str(Rth)])
    
    for r=1:length(roi)
        subplot(5,8,r); hold on;
        errorbar(...
            binned_data(binned_data(:,1)==r,2),...
            binned_data(binned_data(:,1)==r,3),...
            binned_data(binned_data(:,1)==r,4),...
            'LineStyle','none','Marker','o','MarkerSize',6,...
            'Color','k','MarkerFaceColor',[.3 .3 .3]);
        plot([8 8],[0 25],'k--')
        
        if ~isempty(lm{r}) && ~isnan(lm{r}.Coefficients.pValue(1))
            x=[0;20];
            [ypred,yci] = predict(lm{r},x);
            yci2=yci(:,2)-ypred;
            if lm{r}.Coefficients.pValue(2) < 0.05
                boundedline(x,ypred',yci2,'r','alpha')
            else
                boundedline(x,ypred',yci2,'k','alpha')
            end
        end
        
        title([roilabels{r} ', R>' num2str(Rth) ', nVox [' ...
            num2str(roi_n{r,2}) ' ' num2str(roi_n{r,3}) ']'],...
            'interpreter','none');
        %title(roilabels{r});
        xlabel('Ecc (dva)');ylabel('Sz (dva)');
        set(gca, 'Box','off', 'xlim', [0 20], 'ylim',[0 25]);
    end
    
    if SaveFigs
        saveas(f5css,fullfile(pngfld, ['MRI_Ecc_vs_Size_R' num2str(Rth) '.png']));
        saveas(f5css,fullfile(svgfld, ['MRI_Ecc_vs_Size_R' num2str(Rth) '.svg']));
        saveas(f5css_binned,fullfile(pngfld, ['MRI_Ecc_vs_Size_binned_R' num2str(Rth) '.png']));
        saveas(f5css_binned,fullfile(svgfld, ['MRI_Ecc_vs_Size_binned_R' num2str(Rth) '.svg']));
    end
    if CloseFigs; close(f5css_R5); close(f5css_binned_R5); end
    
    % Select areas with data in both monkeys =================================
    if false
        minVox_both=25;
        % which ROIs have a minimum number of voxels in both animals?
        s_roi = mean(nvox_roi>=minVox_both,2)==1;
        fprintf(['ROIs with more than ' num2str(minVox_both) ' sign. voxels in BOTH:\n'])
        rois(s_roi,1)
        
        minVox_s1=minVox_both;
        s1_roi = nvox_roi(:,1)>minVox_s1;
        fprintf(['ROIs with more than ' num2str(minVox_s1) ' sign. voxels in S1:\n'])
        rois(s1_roi,1)
        
        minVox_s2=minVox_both;
        s2_roi = nvox_roi(:,2)>minVox_s2;
        fprintf(['ROIs with more than ' num2str(minVox_s2) ' sign. voxels in S2:\n'])
        rois(s2_roi,1)
        
        
        figure; hold on;
        lidx=1; clear L
        for idx=find(s_roi==1)'
            x=[0;20];
            [ypred,yci] = predict(lm{idx},x);
            yci2=yci(:,2)-ypred;
            plot(x,ypred','LineWidth',3)
            L{lidx}=roilabels{idx};
            lidx=lidx+1;
        end
        legend(L)
    end
    
    % NEW Fig 5: Ecc-Sz for CSS, selected ROIs ===============================
    A_idx = [1:5 8:9];
    B_idx = [10:13 23 25:26];
    C_idx = [19 18 20];
    D_idx = [29:31 33 37:38];
    
    
    % precreate color matrices for plotting
    A_cols = brewermap(length(A_idx),def_cmap);
    B_cols = brewermap(length(B_idx),def_cmap);
    C_cols = brewermap(length(C_idx),def_cmap);
    D_cols = brewermap(length(D_idx),def_cmap);
    
    f_eccsz = figure;
    set(f_eccsz,'Position',[100 100 2400 2000]);
    sgtitle(['Ecc vs Sz, R>' num2str(Rth)]) 
    
    subplot(2,2,1); hold on;
    lidx=1; clear L
    for r=A_idx
        errorbar(...
            binned_data(binned_data(:,1)==r,2),...
            binned_data(binned_data(:,1)==r,3),...
            binned_data(binned_data(:,1)==r,4),...
            'LineStyle','none','Marker','o','MarkerSize',6,...
            'Color',A_cols(lidx,:),'MarkerFaceColor',A_cols(lidx,:));
        lidx=lidx+1;
    end
    lidx=1;
    for r=A_idx
        x=[0;20];
        [ypred,yci] = predict(lm{r},x);
        yci2=yci(:,2)-ypred;
        %boundedline(x,ypred',yci2,'cmap',A_cols(lidx,:),'alpha')
        plot(x,ypred','-','LineWidth',2,'Color',A_cols(lidx,:))
        L{lidx}=roilabels{r};
        lidx=lidx+1;
    end
    legend(L,'Location','NorthWest')
    plot([8 8],[0 25],'k--')
    set(gca,'ylim',[0 15], 'TickDir','out');
    
    subplot(2,2,2); hold on;
    lidx=1; clear L
    for r=B_idx
        errorbar(...
            binned_data(binned_data(:,1)==r,2),...
            binned_data(binned_data(:,1)==r,3),...
            binned_data(binned_data(:,1)==r,4),...
            'LineStyle','none','Marker','o','MarkerSize',6,...
            'Color',B_cols(lidx,:),'MarkerFaceColor',B_cols(lidx,:));
        lidx=lidx+1;
    end
    lidx=1;
    for r=B_idx
        x=[0;20];
        [ypred,yci] = predict(lm{r},x);
        yci2=yci(:,2)-ypred;
        %boundedline(x,ypred',yci2,'cmap',B_cols(lidx,:),'alpha')
        plot(x,ypred','-','LineWidth',2,'Color',B_cols(lidx,:))
        L{lidx}=roilabels{r};
        lidx=lidx+1;
    end
    legend(L,'Location','NorthWest')
    plot([8 8],[0 25],'k--')
    set(gca,'ylim',[0 20], 'TickDir','out');
    
    C_idx = [19 18 20];
    D_idx = [29:31 33 37 38];
    
    subplot(2,2,3); hold on;
    lidx=1; clear L
    for r=C_idx
        errorbar(...
            binned_data(binned_data(:,1)==r,2),...
            binned_data(binned_data(:,1)==r,3),...
            binned_data(binned_data(:,1)==r,4),...
            'LineStyle','none','Marker','o','MarkerSize',6,...
            'Color',C_cols(lidx,:),'MarkerFaceColor',C_cols(lidx,:));
        lidx=lidx+1;
    end
    lidx=1;
    for r=C_idx
        x=[0;20];
        [ypred,yci] = predict(lm{r},x);
        yci2=yci(:,2)-ypred;
        %boundedline(x,ypred',yci2,'cmap',C_cols(lidx,:),'alpha')
        plot(x,ypred','-','LineWidth',2,'Color',C_cols(lidx,:))
        L{lidx}=roilabels{r};
        lidx=lidx+1;
    end
    legend(L,'Location','NorthWest')
    plot([8 8],[0 25],'k--')
    set(gca,'ylim',[0 20], 'TickDir','out');
    
    subplot(2,2,4); hold on;
    lidx=1; clear L
    for r=D_idx
        errorbar(...
            binned_data(binned_data(:,1)==r,2),...
            binned_data(binned_data(:,1)==r,3),...
            binned_data(binned_data(:,1)==r,4),...
            'LineStyle','none','Marker','o','MarkerSize',6,...
            'Color',D_cols(lidx,:),'MarkerFaceColor',D_cols(lidx,:));
        lidx=lidx+1;
    end
    lidx=1;
    for r=D_idx
        x=[0;20];
        [ypred,yci] = predict(lm{r},x);
        yci2=yci(:,2)-ypred;
        %boundedline(x,ypred',yci2,'cmap',D_cols(lidx,:),'alpha')
        plot(x,ypred','-','LineWidth',2,'Color',D_cols(lidx,:))
        L{lidx}=roilabels{r};
        lidx=lidx+1;
    end
    legend(L,'Location','NorthWest')
    plot([8 8],[0 25],'k--')
    set(gca,'ylim',[0 20], 'TickDir','out');
    
    if SaveFigs
        saveas(f_eccsz,fullfile(pngfld, ['MRI_Ecc_vs_Size_ROIs_R' num2str(Rth) '.png']));
        saveas(f_eccsz,fullfile(svgfld, ['MRI_Ecc_vs_Size_ROIs_R' num2str(Rth) '.svg']));
    end
    if CloseFigs; close(f_eccsz); end
    
end

%% SFig 5A/B: Ecc-Sz per area including CI and binned data with SEM =======
reroi = [19 18 20 1:5 7:14 21:38];
for Rth=[5 3]
    
    clear tbl lm;
    bins=[0:18; 2:20]';
    
    sf5css=figure;
    set(sf5css,'Position',[100 100 1600 1200]);
    set(sf5css,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap));
    sgtitle(['R>' num2str(Rth)])
    
    m=3; % css
    s_R2 = T(modidx.(MRI_MODELS{m,1})).mod.R2 > Rth;
    
    nsp=length(reroi);
   
    nvox_roi = []; binned_data=[];
    for r=1:length(reroi)
        rr=reroi(r);
        ES2{r} = []; ES2_m{r} = [];
        ES2b{r} = []; ES2f{r} = [];
        
        SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,...
            ck_GetROIidx(roilabels(rr),rois) );
        
        s_Monkey = strcmp(T(modidx.(MRI_MODELS{m,1})).mod.Monkey,'M01');
        SSS_m1 = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,...
            ck_GetROIidx(roilabels(rr),rois)) & s_Monkey;
        
        s_Monkey = strcmp(T(modidx.(MRI_MODELS{m,1})).mod.Monkey,'M02');
        SSS_m2 = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,...
            ck_GetROIidx(roilabels(rr),rois)) & s_Monkey;
        
        nvox_roi = [nvox_roi; sum(SSS_m1) sum(SSS_m2)];
        roi_n{r,1} = rois{rr};
        roi_n{r,2} = sum(SSS_m1);
        roi_n{r,3} = sum(SSS_m2);
        
        ES2{r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS) ...
            T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS) ];
        ES2{r}( ES2{r}(:,2) > 50,: ) = []; % remove insanely large pRF's
        
        ES2_m{1,r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS_m1) ...
            T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS_m1) ];
        ES2_m{2,r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS_m2) ...
            T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS_m2) ];
        
        subplot(5,7,r); hold on;
        
        scatter(ES2_m{1,r}(:,1),ES2_m{1,r}(:,2),20,'Marker','o',...
            'MarkerFaceColor','b','MarkerFaceAlpha',.2,...
            'MarkerEdgeColor','none');
        scatter(ES2_m{2,r}(:,1),ES2_m{2,r}(:,2),20,'Marker','o',...
            'MarkerFaceColor','k','MarkerFaceAlpha',.2,...
            'MarkerEdgeColor','none');
        
        % make table
        tbl2 = table(ES2{r}(:,1),ES2{r}(:,2),'VariableNames',{'ecc','sz'});
        if ~isempty(tbl2)
            % remove inf's
            tbl2(isinf(tbl2.sz),:)=[];
            % fit
            lm2{r} = fitlm(tbl2,'sz ~ 1 + ecc');
            % CI
            lmCI2{r}=coefCI(lm2{r},0.05);
            if ~isnan(lm2{r}.Coefficients.pValue(1))
                % Prediction
                x=[0;20];
                [ypred,yci] = predict(lm2{r},x);
                yci2=yci(:,2)-ypred;
                if lm2{r}.Coefficients.pValue(2) < 0.05
                    boundedline(x,ypred',yci2,'r','alpha')
                else
                    boundedline(x,ypred',yci2,'k','alpha')
                end
                fprintf([roilabels{rr} ' n = ' num2str(size(tbl2,1)) '\n'])
                fprintf([roilabels{rr} ' p = ' num2str(lm2{r}.Coefficients.pValue(2)) '\n'])
                fprintf([roilabels{rr} ' slope = ' num2str(lm2{r}.Coefficients.Estimate(2)) '\n'])
            end
            
            % binned
            for b=1:size(bins,1)
                seldata = (tbl2.ecc>bins(b,1) & tbl2.ecc<=bins(b,2)) ;
                if sum(seldata)>1
                    binned_data=[binned_data;...
                        rr b mean(tbl2.sz(seldata)) std(tbl2.sz(seldata))./sqrt(sum(seldata))];
                end
            end
        end
        
        errorbar(...
            binned_data(binned_data(:,1)==rr,2),...
            binned_data(binned_data(:,1)==rr,3),...
            binned_data(binned_data(:,1)==rr,4),...
            'LineStyle','none','Marker','o','MarkerSize',6,...
            'Color','k','MarkerFaceColor',[.8 .8 .8]);
        plot([8 8],[0 25],'k--')
        
        %xlabel('Eccentricity');ylabel('pRF size');
        set(gca, 'Box','off', 'xlim', [0 20], 'ylim',[0 25],...
            'FontName','Myriad Pro','FontSize',11,'TitleHorizontalAlignment','left');
        title(['Ecc vs pRF size [' roilabels{rr} ', R>' num2str(Rth) ']'],...
            'interpreter','none');
        title(roilabels{rr});
        text(2,20,sprintf(['n = [' num2str(roi_n{r,2}) ' ' num2str(roi_n{r,3}) ']\n'...
            'p = ' num2str(lm2{r}.Coefficients.pValue(2)) '\n'...
            'slope = ' num2str(lm2{r}.Coefficients.Estimate(2))]));
    end
        
    if SaveFigs
        saveas(sf5css,fullfile(pngfld, ['SF5_MRI_Ecc_vs_Size_R' num2str(Rth) '.png']));
        saveas(sf5css,fullfile(svgfld, ['SF5_MRI_Ecc_vs_Size_R' num2str(Rth) '.svg']));
    end
    if CloseFigs; close(sf5css); end    
end



%% ========================================================================
%  EPHYS PRF ANALYSIS -----------------------------------------------------
%  ========================================================================

%% Fig 9: EPHYS ECC vs PRF Size ===========================================
Rth=25; 
SNRth=3;

m=unique(tMUA.Model);
R2=[];
for i=2%1:length(m)
    R2 = [R2 tMUA.R2(strcmp(tMUA.Model,m{i}))];
end

subs = tMUA.Monkey;
for i=1:length(subs)
    if strcmp(subs{i},'M03')
        subs{i} = 1;
    elseif strcmp(subs{i},'M04')
        subs{i} = 2;
    end
end
subs=cell2mat(subs); names={'M03','M04'};

mm=unique(tMUA.Model);
for m=2%1:length(mm)
    % collect
    if ~strcmp(mm{m},'classicRF')
        eccsz = [tMUA.R2(strcmp(tMUA.Model,mm{m})) ...
            tMUA.ecc(strcmp(tMUA.Model,mm{m})) ...
            tMUA.rfs(strcmp(tMUA.Model,mm{m})) ...
            tMUA.Area(strcmp(tMUA.Model,mm{m})) ...
            tMUA.Array(strcmp(tMUA.Model,mm{m})) ...
            subs(strcmp(tMUA.Model,mm{m})) ...
            tMUA.gain(strcmp(tMUA.Model,mm{m})) ...
            ];
    else 
        eccsz = [tMUA.SNR(strcmp(tMUA.Model,mm{m})) ...
            tMUA.ecc(strcmp(tMUA.Model,mm{m})) ...
            tMUA.rfs(strcmp(tMUA.Model,mm{m})) ...
            tMUA.Area(strcmp(tMUA.Model,mm{m})) ...
            tMUA.Array(strcmp(tMUA.Model,mm{m})) ...
            subs(strcmp(tMUA.Model,mm{m})) ...
            ];
    end
    
    f_eccsz_arr = figure;
    set(f_eccsz_arr,'Position',[100 100 1800 600]);
    for ss=1:2
        a=[1 4];
        for aa=1:length(a)
            leg={};
            if strcmp(mm{m},'classicRF')
                sel = eccsz(:,1)>SNRth & eccsz(:,4)==a(aa) & eccsz(:,6)==ss;
            else
                sel = eccsz(:,1)>Rth & eccsz(:,4)==a(aa) & eccsz(:,6)==ss;
            end
            es1=eccsz(sel,:);
            subplot(2,2,(ss-1)*2 + aa); 
            hold on;
            arr=unique(es1(:,5));
            for aaa = 1:length(arr)
                if arr(aaa)~=99
                    sel2=es1(:,5)==arr(aaa);
                    scatter(es1(sel2,2),es1(sel2,3));
                    leg=[leg {num2str(arr(aaa))}];
                end
            end
            if a(aa)==1
                set(gca,'xlim',[0 8],'ylim',[0 3]);
            else
                set(gca,'xlim',[0 8],'ylim',[0 5]);
            end
            legend(leg,'Location','NorthEastOutside');
            title([names{ss} ' V' num2str(a(aa))]);
            xlabel('Ecc (dva)'); ylabel('RF-size (sd, dva)')
        end
    end
    sgtitle(['MODEL: ' mm{m} ', R2th: ' num2str(Rth)],...
        'interpreter','none');
    if SaveFigs
        saveas(f_eccsz_arr,fullfile(svgfld, ...
            ['EPHYS_MUA_Ecc_vs_Size_' mm{m} '.svg']));
        saveas(f_eccsz_arr,fullfile(pngfld, ...
            ['EPHYS_MUA_Ecc_vs_Size_' mm{m} '.png']));
    end
    if CloseFigs; close(f_eccsz_arr); end
    
    f_eccsz_v1 = figure;
    set(f_eccsz_v1,'Position',[100 100 1600 400]);
    for ss=1:2
        leg={};
        if strcmp(mm{m},'classicRF')
            sel = eccsz(:,1)>SNRth & eccsz(:,4)==1 & eccsz(:,6)==ss;
        else
            sel = eccsz(:,1)>Rth & eccsz(:,4)==1 & eccsz(:,6)==ss;
        end
        es1=eccsz(sel,:);
        subplot(1,2,ss); hold on;
        arr=unique(es1(:,5));
        for aaa = 1:length(arr)
            sel2=es1(:,5)==arr(aaa);
            scat=scatter(es1(sel2,2),es1(sel2,3));
            %if (ss==1 && (arr(aaa)==11 || arr(aaa)==13)) || ...
            %        (ss==2 && (arr(aaa)==10 || arr(aaa)==12))
            if (ss==1 && arr(aaa)==11) || ...
                    (ss==2 && arr(aaa)==10)
                set(scat,'MarkerFaceColor','r','MarkerEdgeColor','r');
            elseif (ss==1 && arr(aaa)==13) || ...
                    (ss==2 && arr(aaa)==12)
                set(scat,'MarkerFaceColor','b','MarkerEdgeColor','b');
                
            else
                set(scat,'MarkerFaceColor','k','MarkerEdgeColor','k');
            end
            leg=[leg {num2str(arr(aaa))}];
        end
        title([names{ss} ' V1']);
        xlabel('Ecc (dva)'); ylabel('RF-size (sd, dva)')
        set(gca,'xlim',[0 8],'ylim',[0 3]);
        %legend(leg,'Location','NorthEastOutside');
    end
    sgtitle(['MODEL: ' mm{m} ', R2th: ' num2str(Rth)],...
        'interpreter','none');
    if SaveFigs
        saveas(f_eccsz_v1,fullfile(pngfld, ...
            ['EPHYS_MUA_Ecc_vs_Size_' mm{m} '.png']));
        saveas(f_eccsz_v1,fullfile(svgfld, ...
            ['EPHYS_MUA_Ecc_vs_Size_' mm{m} '.svg']));
    end
    if CloseFigs; close(f_eccsz_v1); end
      
end
%close all;

%% Fig 6 / SFig 6 EPHYS VFC MUA (scatter and heatmap) =====================
% scatter MUA =============================================================
R2TH = 50;
s=tMUA.R2>R2TH;
model='css_ephys_cv1';

fm=figure; set(fm,'Position',[100 100 1400 1000]);
[cmap,~,~] = brewermap(length(unique(tMUA.Array)),'spectral');
set(fm,'DefaultAxesColorOrder',cmap);
msz = 5;

% M03 V1 ----
m=strcmp(tMUA.Monkey,'M03') & strcmp(tMUA.Model,model);
v=tMUA.Area==1;

subplot(2,2,1);hold on;
plot([-10 20],[0 0],'k'); plot([0 0],[-30 10],'k');
lt={'meridian','meridian'};
nelec=0;
for r=unique(tMUA.Array)'
    a=tMUA.Array==r;
    currcol = cmap(r,:);
    if sum(s & m & a & v) > 0
        nelec=nelec+sum(s & m & a & v);
        p=plot(tMUA.X(s & m & a & v),...
            tMUA.Y(s & m & a & v),'o','Color',currcol,...
            'LineStyle','none','LineWidth',1,...
            'MarkerSize',msz,'MarkerFaceColor',currcol);
        lt=[lt num2str(r)];
    end
end
set(gca, 'Box','off', 'xlim', [-1 8], 'ylim',[-6 1]);
title(['M03 V1 MUA (n = ' num2str(nelec) ')'])
l=legend(lt); set(l,'Location','NorthEastOutside'); clear lt;

% M03 V4 ----
v=tMUA.Area==4;

subplot(2,2,3);hold on;
plot([-10 20],[0 0],'k'); plot([0 0],[-30 10],'k');
lt={'meridian','meridian'};
nelec=0;
for r=unique(tMUA.Array)'
    a=tMUA.Array==r;
    currcol = cmap(r,:);
    if sum(s & m & a & v) > 0
        nelec=nelec+sum(s & m & a & v);
        p=plot(tMUA.X(s & m & a & v),...
            tMUA.Y(s & m & a & v),'o','Color',currcol,...
            'LineStyle','none','LineWidth',1,...
            'MarkerSize',msz,'MarkerFaceColor',currcol);
        lt=[lt num2str(r)];
    end
end
set(gca, 'Box','off', 'xlim', [-1 8], 'ylim',[-6 1]);
% set(gca, 'Box','off', 'xlim', [-2 25], 'ylim',[-25 2]);
title(['M03 V4 MUA (n = ' num2str(nelec) ')'])
l=legend(lt); set(l,'Location','NorthEastOutside'); clear lt;

% M04 V1 ----
m=strcmp(tMUA.Monkey,'M04') & strcmp(tMUA.Model,model);
v=tMUA.Area==1;

subplot(2,2,2);hold on;
plot([-10 20],[0 0],'k'); plot([0 0],[-30 10],'k');
lt={'meridian','meridian'};
nelec=0;
for r=unique(tMUA.Array)'
    a=tMUA.Array==r;
    currcol = cmap(r,:);
    if sum(s & m & a & v) > 0
        nelec=nelec+sum(s & m & a & v);
        p=plot(tMUA.X(s & m & a & v),...
            tMUA.Y(s & m & a & v),'o','Color',currcol,...
            'LineStyle','none','LineWidth',1,...
            'MarkerSize',msz,'MarkerFaceColor',currcol);
        lt=[lt num2str(r)];
    end
end
set(gca, 'Box','off', 'xlim', [-1 8], 'ylim',[-6 1]);
title(['M04 V1 MUA (n = ' num2str(nelec) ')'])
l=legend(lt); set(l,'Location','NorthEastOutside'); clear lt;

% M04 V4 ----
v=tMUA.Area==4;

subplot(2,2,4);hold on;
plot([-10 30],[0 0],'k'); plot([0 0],[-30 10],'k');
lt={'meridian','meridian'};
nelec=0;
for r=unique(tMUA.Array)'
    a=tMUA.Array==r;
    currcol = cmap(r,:);
    if sum(s & m & a & v) > 0
        nelec=nelec+sum(s & m & a & v);
        p=plot(tMUA.X(s & m & a & v),...
            tMUA.Y(s & m & a & v),'o','Color',currcol,...
            'LineStyle','none','LineWidth',1,...
            'MarkerSize',msz,'MarkerFaceColor',currcol);
        lt=[lt num2str(r)];
    end
end
set(gca, 'Box','off', 'xlim', [-2 25], 'ylim',[-25 2]);
title(['M04 V4 MUA (n = ' num2str(nelec) ')'])
l=legend(lt); set(l,'Location','NorthEastOutside'); clear lt;

if SaveFigs
    saveas(fm,fullfile(pngfld, 'MUA_VFC.png'));
    saveas(fm,fullfile(svgfld, 'MUA_VFC.svg'));
end
if CloseFigs; close(fm); end


% Heatmap MUA =============================================================
fhm=figure; set(fhm,'Position',[100 100 1800 600]);
settings.PixPerDeg = 29.5032;
settings.meshsize = 2000;   
colormap(inferno) 

% M03 V1 ----
m=strcmp(tMUA.Monkey,'M03') & strcmp(tMUA.Model,model);
v=tMUA.Area==1;

subplot(2,4,1);hold on;
allprf=[];
for r=unique(tMUA.Array)'
    a=tMUA.Array==r;
    if sum(s & m & a & v) > 0
        allprf = [allprf; ...
            tMUA.X(s & m & a & v) ...
            tMUA.Y(s & m & a & v) ...
            tMUA.rfs(s & m & a & v)];
    end
end

res = ck_2dPRF_ephys(allprf(:,1),allprf(:,2),allprf(:,3));
img=flipud(res.img); res.ymesh2 = fliplr(res.ymesh);
% zoom in on plot
xrange = [-3 8]; yrange = [-6 3];
xrange_idx = [...
    find(res.xmesh >= xrange(1),1,'first') find(res.xmesh <= xrange(2),1,'last')];
yrange_idx = [...
    find(res.ymesh2 >= yrange(1),1,'first') find(res.ymesh2 <= yrange(2),1,'last')];
xrr=res.xmesh(xrange_idx); yrr=res.ymesh2(yrange_idx);
% plot
sumimg=sum(img,3);
imagesc(sumimg);
set(gca,'xlim',xrange_idx,'ylim',...
    yrange_idx,'Color','k','xtick',[],'ytick',[])
colorbar;
subplot(2,4,5); hold on;
plot(1.2*res.xr,[0 0],'w'); plot([0 0],1.2*res.yr,'w');
set(gca,'xlim',xrr,'ylim',yrr,'Color','k')
colorbar; clear('res','img','sumimg');
title('M03 V1 MUA')

% M03 V4 ----
v=tMUA.Area==4;

subplot(2,4,2);hold on;
allprf=[];
for r=unique(tMUA.Array)'
    a=tMUA.Array==r;
    if sum(s & m & a & v) > 0
        allprf = [allprf; ...
            tMUA.X(s & m & a & v) ...
            tMUA.Y(s & m & a & v) ...
            tMUA.rfs(s & m & a & v)];
    end
end

res = ck_2dPRF_ephys(allprf(:,1),allprf(:,2),allprf(:,3));
img=flipud(res.img); res.ymesh2 = fliplr(res.ymesh);
% zoom in on plot
xrange = [-3 8]; yrange = [-6 3];
xrange_idx = [...
    find(res.xmesh >= xrange(1),1,'first') find(res.xmesh <= xrange(2),1,'last')];
yrange_idx = [...
    find(res.ymesh2 >= yrange(1),1,'first') find(res.ymesh2 <= yrange(2),1,'last')];
xrr=res.xmesh(xrange_idx); yrr=res.ymesh2(yrange_idx);
% plot
sumimg=sum(img,3);
imagesc(sumimg);
set(gca,'xlim',xrange_idx,'ylim',...
    yrange_idx,'Color','k','xtick',[],'ytick',[])
colorbar;
subplot(2,4,6); hold on;
plot(1.2*res.xr,[0 0],'w'); plot([0 0],1.2*res.yr,'w');
set(gca,'xlim',xrr,'ylim',yrr,'Color','k')
colorbar; clear('res','img','sumimg');
title('M03 V4 MUA')

% M04 V1 ----
m=strcmp(tMUA.Monkey,'M04') & strcmp(tMUA.Model,model);
v=tMUA.Area==1;

subplot(2,4,3);hold on;
allprf=[];
for r=unique(tMUA.Array)'
    a=tMUA.Array==r;
    if sum(s & m & a & v) > 0
        allprf = [allprf; ...
            tMUA.X(s & m & a & v) ...
            tMUA.Y(s & m & a & v) ...
            tMUA.rfs(s & m & a & v)];
    end
end

res = ck_2dPRF_ephys(allprf(:,1),allprf(:,2),allprf(:,3));
img=flipud(res.img); res.ymesh2 = fliplr(res.ymesh);
% zoom in on plot
xrange = [-3 8]; yrange = [-6 3];
xrange_idx = [...
    find(res.xmesh >= xrange(1),1,'first') find(res.xmesh <= xrange(2),1,'last')];
yrange_idx = [...
    find(res.ymesh2 >= yrange(1),1,'first') find(res.ymesh2 <= yrange(2),1,'last')];
xrr=res.xmesh(xrange_idx); yrr=res.ymesh2(yrange_idx);
% plot
sumimg=sum(img,3);
imagesc(sumimg);
set(gca,'xlim',xrange_idx,'ylim',...
    yrange_idx,'Color','k','xtick',[],'ytick',[])
colorbar;
subplot(2,4,7); hold on;
plot(1.2*res.xr,[0 0],'w'); plot([0 0],1.2*res.yr,'w');
set(gca,'xlim',xrr,'ylim',yrr,'Color','k')
colorbar; clear('res','img','sumimg');
title('M04 V1 MUA')

% M04 V4 ----
v=tMUA.Area==4;

subplot(2,4,4);hold on;
allprf=[];
for r=unique(tMUA.Array)'
    a=tMUA.Array==r;
    if sum(s & m & a & v) > 0
        allprf = [allprf; ...
            tMUA.X(s & m & a & v) ...
            tMUA.Y(s & m & a & v) ...
            tMUA.rfs(s & m & a & v)];
    end
end

res = ck_2dPRF_ephys(allprf(:,1),allprf(:,2),allprf(:,3));
img=flipud(res.img); res.ymesh2 = fliplr(res.ymesh);
% zoom in on plot
xrange = [-5 25]; yrange = [-25 5];
xrange_idx = [...
    find(res.xmesh >= xrange(1),1,'first') find(res.xmesh <= xrange(2),1,'last')];
yrange_idx = [...
    find(res.ymesh2 >= yrange(1),1,'first') find(res.ymesh2 <= yrange(2),1,'last')];
xrr=res.xmesh(xrange_idx); yrr=res.ymesh2(yrange_idx);
% plot
sumimg=sum(img,3);
imagesc(sumimg);
set(gca,'xlim',xrange_idx,'ylim',...
    yrange_idx,'Color','k','xtick',[],'ytick',[])
colorbar;
subplot(2,4,8); hold on;
plot(1.2*res.xr,[0 0],'w'); plot([0 0],1.2*res.yr,'w');
set(gca,'xlim',xrr,'ylim',yrr,'Color','k')
colorbar; clear('res','img','sumimg');
title('M04 V4 MUA')

if SaveFigs
    saveas(fhm,fullfile(pngfld, 'MUA_VFC_hm.png'));
    saveas(fhm,fullfile(svgfld, 'MUA_VFC_hm.svg'));
end
if CloseFigs; close(fhm); end

%% Ephys location difference & size difference ============================
rth=25; snrth=3; c_rth=25;

ephys_MMS = MMS(:,1);
clear C R2m SZ
for m=1:length(ephys_MOD)
 
    C{m}=[];R2m{m}=[];SZ{m}=[];
    
    model=ephys_MOD{m};
    s = strcmp(tMUA.Model,model);
    C{m}=[C{m} tMUA.R2(s) tMUA.X(s) tMUA.Y(s) tMUA.rfs(s)];
    R2m{m}=[R2m{m} tMUA.R2(s)];
    SZ{m}=[SZ{m} tMUA.R2(s) tMUA.rfs(s)];
    
    PRF_EST(m,1).sig = 'MUA';
    PRF_EST(m,1).R2 = tMUA.R2(s);
    PRF_EST(m,1).X = tMUA.X(s);
    PRF_EST(m,1).Y = tMUA.Y(s);
    PRF_EST(m,1).S = tMUA.rfs(s);
    PRF_EST(m,1).A = tMUA.Area(s);
    PRF_EST(m,1).G = tMUA.gain(s);
    PRF_EST(m,1).M = tMUA.Monkey(s);
    PRF_EST(m,1).ARR = tMUA.Array(s);
    
    s = strcmp(tMUA.Model,'classicRF');
    %C{m}=[C{m} tMUA.X(s)./668.745 tMUA.Y(s)./668.745 tMUA.rfs(s)./2];
    
    PRF_EST(m,2).sig = 'MUACLASSIC';
    PRF_EST(m,2).R2 = tMUA.R2(s);
    PRF_EST(m,2).X = tMUA.X(s);
    PRF_EST(m,2).Y = tMUA.Y(s);
    PRF_EST(m,2).S = tMUA.rfs(s);
    PRF_EST(m,2).A = tMUA.Area(s);
    PRF_EST(m,2).G = tMUA.gain(s);
    PRF_EST(m,2).M = tMUA.Monkey(s);
    PRF_EST(m,2).SNR = tMUA.SNR(s);
    PRF_EST(m,2).ARR = tMUA.Array(s);

    s = strcmp(tLFP.Model,model);
    sig=unique(tLFP.SigType);
    lfp_order = [3 1 2 5 4];
    cnt=1;
    for i=lfp_order
        b=strcmp(tLFP.SigType,sig{i});
        C{m}=[C{m} tLFP.R2(s & b) tLFP.X(s & b) tLFP.Y(s & b) tLFP.rfs(s & b)];
        R2m{m}=[R2m{m} tLFP.R2(s & b)];
        SZ{m}=[SZ{m} tLFP.R2(s & b) tLFP.rfs(s & b)];
        
        PRF_EST(m,2+cnt).sig = sig{i};
        PRF_EST(m,2+cnt).R2 = tLFP.R2(s & b);
        PRF_EST(m,2+cnt).X = tLFP.X(s & b);
        PRF_EST(m,2+cnt).Y = tLFP.Y(s & b);
        PRF_EST(m,2+cnt).S =  tLFP.rfs(s & b);
        PRF_EST(m,2+cnt).A =  tLFP.Area(s & b);
        PRF_EST(m,2+cnt).G =  tLFP.gain(s & b);
        
        PRF_EST(m,2+cnt).M = tLFP.Monkey(s & b);
        
        cnt=cnt+1;
    end
    
    s= sum(R2m{m}>rth,2)==size(R2m{m},2);
    distRF{m} = [...
        sqrt(((C{m}(s,2)-C{m}(s,5)).^2) + ((C{m}(s,3)-C{m}(s,6)).^2)) ...
        sqrt(((C{m}(s,2)-C{m}(s,9)).^2) + ((C{m}(s,3)-C{m}(s,10)).^2)) ...
        sqrt(((C{m}(s,2)-C{m}(s,13)).^2) + ((C{m}(s,3)-C{m}(s,14)).^2)) ...
        sqrt(((C{m}(s,2)-C{m}(s,17)).^2) + ((C{m}(s,3)-C{m}(s,18)).^2)) ...
        sqrt(((C{m}(s,2)-C{m}(s,21)).^2) + ((C{m}(s,3)-C{m}(s,22)).^2)) ];
    
    distSZ{m} = [...
        C{m}(s,4)-C{m}(s,7) ...
        C{m}(s,4)-C{m}(s,11) ...
        C{m}(s,4)-C{m}(s,15) ...
        C{m}(s,4)-C{m}(s,19) ...
        C{m}(s,4)-C{m}(s,23) ];
    
    % normalize by MUA-pRF
    normSz{m} = [...
        C{m}(s,7)./C{m}(s,4) ...
        C{m}(s,11)./C{m}(s,4) ...
        C{m}(s,15)./C{m}(s,4) ...
        C{m}(s,19)./C{m}(s,4) ...
        C{m}(s,23)./C{m}(s,4) ];
    
%     % normalize by cRF
%     normSz{m} = [...
%         C{m}(s,7)./PRF_EST(m,2).S(s,:) ...
%         C{m}(s,11)./PRF_EST(m,2).S(s,:) ...
%         C{m}(s,15)./PRF_EST(m,2).S(s,:) ...
%         C{m}(s,19)./PRF_EST(m,2).S(s,:) ...
%         C{m}(s,23)./PRF_EST(m,2).S(s,:) ];
    
end
 
%% Fig 7: compare MUA with classic ========================================
for m=[1 3] 
    mlabel=ephys_MMS{m};
    cRF_coll = [];
    pRF_coll = [];

    f_cRFpRF=figure;
    set(f_cRFpRF,'Position',[100 100 1500 600]);
    
    pn=0;
    for area = [1 4]
        pn=pn+1;
        fprintf(['=== AREA V' num2str(area) ' ===\n'])
        % location
        sel = PRF_EST(m,1).A == area & ...
            PRF_EST(m,1).R2 > rth & ...
            PRF_EST(m,2).R2 >= c_rth & ...
            PRF_EST(m,2).SNR >= snrth & ...
            ( strcmp(PRF_EST(m,1).M,'M03') & (PRF_EST(m,1).ARR ~= 11 & PRF_EST(m,1).ARR ~= 13) | ... 
            strcmp(PRF_EST(m,1).M,'M04') & (PRF_EST(m,1).ARR ~= 10 & PRF_EST(m,1).ARR ~= 12));
        dLOCATION = sqrt(...
            (PRF_EST(m,1).X(sel)-PRF_EST(m,2).X(sel)).^2 + ...
            (PRF_EST(m,1).Y(sel)-PRF_EST(m,2).Y(sel)).^2);
        
        mselect = PRF_EST(m,1).M(sel);
        m3idx = strcmp(mselect,'M03'); m4idx = strcmp(mselect,'M04');
        
        fprintf(['MODEL ' ephys_MOD{m} ', MUA vs CLASSIC distance ----\n'])
        fprintf(['Mean ' num2str(nanmean(dLOCATION)) ', STD ' num2str(nanstd(dLOCATION)) '\n'])
        fprintf(['Median ' num2str(nanmedian(dLOCATION)) ', IQR ' num2str(iqr(dLOCATION)) '\n'])

        % size
        SIZE_MUA = [PRF_EST(m,1).S(sel) PRF_EST(m,2).S(sel)];
        rmINF = isinf(SIZE_MUA(:,1));
        SIZE_MUA(rmINF,:)=[];
        m3idx(rmINF,:)=[];
        m4idx(rmINF,:)=[];
        
        [p,h,stats] = signrank(SIZE_MUA(:,1),SIZE_MUA(:,2));
        fprintf(['MODEL ' ephys_MOD{m} ', MUA vs CLASSIC size ----\n'])
        fprintf(['dSZ: Wilcoxon z = ' ...
            num2str(stats.zval) ', p = ' num2str(p) '\n']);
        fprintf(['Mean (mua-classic) ' num2str(nanmean(SIZE_MUA(:,1)-SIZE_MUA(:,2))) ...
            ', STD ' num2str(nanstd(SIZE_MUA(:,1)-SIZE_MUA(:,2))) '\n'])
        fprintf(['Median ' num2str(nanmedian(SIZE_MUA(:,1)-SIZE_MUA(:,2))) ...
            ', IQR ' num2str(iqr(SIZE_MUA(:,1)-SIZE_MUA(:,2))) '\n'])
        
        subplot(1,2,pn);hold on;
        splim=[3 8];
        plot([0 10],[0 10],'-r');

        % split by monkey
        pRF3=SIZE_MUA(m3idx,1); cRF3=SIZE_MUA(m3idx,2);
        scatter(pRF3,cRF3,50,'Marker','o',...
            'MarkerEdgeColor','none',...
            'MarkerFaceColor',[.3 .3 .3],'MarkerFaceAlpha',0.3);      
        idx = ~isnan(pRF3) & ~isnan(cRF3);
        polyfit(pRF3(idx),cRF3(idx),1)
        
        pRF4=SIZE_MUA(m4idx,1); cRF4=SIZE_MUA(m4idx,2);
        scatter(pRF4,cRF4,50,'Marker','o',...
            'MarkerEdgeColor','none',...
            'MarkerFaceColor',[.0 .0 .5],'MarkerFaceAlpha',0.3);
        idx = ~isnan(pRF4) & ~isnan(cRF4);
        polyfit(pRF4(idx),cRF4(idx),1)
        
        pRF=[pRF3;pRF4];cRF=[cRF3;cRF4];
        n_electrodes = [length(pRF3) length(pRF4)]
        plot([0 10],[0 10],'-r');
        set(gca, 'xlim',[0 splim(pn)],'ylim',[0 splim(pn)]);
        legend({'x=y','M3','M4'})
        title(['Area V' num2str(area)])
        xlabel('MUA pRF size'); ylabel('Classic RF size')
        
        % Stats per area
        pRF_both = [pRF3;pRF4];
        cRF_both = [cRF3;cRF4];
        [R,p] = corr(cRF_both,pRF_both,'rows','complete');
        fprintf(['Pearson corr for V' num2str(area) ': R=' num2str(R) ', p=' num2str(p) '\n'])
        
        cRF_coll = [cRF_coll; ...
            cRF3 3*ones(size(cRF3)) area*ones(size(cRF3));...
            cRF4 4*ones(size(cRF4)) area*ones(size(cRF4))];
        pRF_coll = [pRF_coll;  ...
            pRF3 3*ones(size(cRF3)) area*ones(size(pRF3));...
            pRF4 4*ones(size(cRF4)) area*ones(size(pRF4))];
    end
    [R,p] = corr(cRF_coll(:,1),pRF_coll(:,1),'rows','complete');
    fprintf(['Pearson corr: R=' num2str(R) ', p=' num2str(p) '\n'])
    
    %%
    f_cRFpRF2 = figure;
    set(f_cRFpRF2,'Position',[0 0 2000 400]);
    
    subplot(1,3,1); hold on
    plot([0 10],[0 10],'-r');
    
    % M3 V1
    s = pRF_coll(:,2) == 3 & pRF_coll(:,3) == 1;
    scatter(pRF_coll(s,1),cRF_coll(s,1),55,'Marker','o',...
        'MarkerEdgeColor','none',...
        'MarkerFaceColor','k','MarkerFaceAlpha',0.3);
    % M4 V1
    s = pRF_coll(:,2) == 4 & pRF_coll(:,3) == 1;
    scatter(pRF_coll(s,1),cRF_coll(s,1),55,'Marker','o',...
        'MarkerEdgeColor','none',...
        'MarkerFaceColor',[.0 .0 .5],'MarkerFaceAlpha',0.3);
    % M3 V4
    s = pRF_coll(:,2) == 3 & pRF_coll(:,3) == 4;
    scatter(pRF_coll(s,1),cRF_coll(s,1),50,'Marker','o',...
        'MarkerFaceColor',[1 1 1],'LineWidth',1,...
        'MarkerEdgeColor','k','MarkerFaceAlpha',0.3);
    % M4 V4
    s = pRF_coll(:,2) == 4 & pRF_coll(:,3) == 4;
    scatter(pRF_coll(s,1),cRF_coll(s,1),50,'Marker','o',...
        'MarkerFaceColor',[1 1 1],'LineWidth',1,...
        'MarkerEdgeColor',[.0 .0 .5],'MarkerFaceAlpha',0.3);
    
    set(gca, 'xlim',[0 splim(pn)],'ylim',[0 splim(pn)]);
    legend({'x=y','M3 V1','M4 V1','M3 V4','M4 V4'})
    xlabel('MUA pRF size'); ylabel('Classic RF size')
    ppp=pRF_coll(:,1);ccc=cRF_coll(:,1);
    title(['pRF model: ' mlabel]);
    
    
    subplot(1,3,2); hold on;
    histogram(ppp-ccc,-5.05:0.1:5.05,'FaceColor',[.5 .5 .5],'LineStyle','none')
    xlabel('pRF size - cRF size');
    ylabel('n channels');
    med_diff = nanmedian(ppp-ccc);
    iqr_diff = [med_diff-(iqr(ppp-ccc)/2) med_diff+(iqr(ppp-ccc)/2)];  
    m_diff = nanmean(ppp-ccc);
    sem_diff = nanstd(ppp-ccc)./sqrt(sum(~isnan(ppp.*ccc)));
    [p,h,stats] = signrank(ppp-ccc,0);
    fprintf(['Difference pRF-cRF: Wilcoxon z = ' ...
        num2str(stats.zval) ', p = ' num2str(p) '\n']);
    
    plot([med_diff med_diff],[0 100],'r');
    plot([m_diff m_diff],[0 100],'r--');
    plot([0 0],[0 100],'k');
    legend({'histogram', ['median ' num2str(med_diff)], ['mean ' num2str(m_diff)]})
    title(['IQR ' num2str(iqr_diff(1)) ' ' num2str(iqr_diff(2)) ...
        '; SEM ' num2str(sem_diff)])
    
    
    subplot(1,3,3); hold on;
    histogram(ppp./ccc,-0.05:0.1:10.05,'FaceColor',[.5 .5 .5],'LineStyle','none')
    xlabel('pRF size / cRF size');
    ylabel('n channels');
    med_rat = nanmedian(ppp./ccc);
    iqr_rat = [med_rat-(iqr(ppp./ccc)/2) med_rat+(iqr(ppp./ccc)/2)];  
    m_rat = nanmean(ppp./ccc);
    sem_rat = nanstd(ppp./ccc)./sqrt(sum(~isnan(ppp.*ccc)));
    [p,h,stats] = signrank(ppp./ccc,1);
    fprintf(['Ratio pRF/cRF: Wilcoxon z = ' ...
        num2str(stats.zval) ', p = ' num2str(p) '\n']);

    plot([med_rat med_rat],[0 100],'r');
    plot([m_rat m_rat],[0 100],'r--');
    plot([1 1],[0 100],'k');
    legend({'histogram', ['median ' num2str(med_rat)], ['mean ' num2str(m_rat)]})
    title(['IQR ' num2str(iqr_rat(1)) ' ' num2str(iqr_rat(2)) ...
        '; SEM ' num2str(sem_rat)])
    
    if SaveFigs
        saveas(f_cRFpRF2,fullfile(pngfld, ['MUA_cRF-pRF_' mlabel '.png']));
        saveas(f_cRFpRF2,fullfile(svgfld, ['MUA_cRF-pRF_' mlabel '.svg']));
        saveas(f_cRFpRF,fullfile(pngfld, ['MUA_cRF-pRF' mlabel '.png']));
        saveas(f_cRFpRF,fullfile(svgfld, ['MUA_cRF-pRF' mlabel '.svg']));
    end
    if CloseFigs; close(f_cRFpRF); close(f_cRFpRF2); end
end

%% Fig 8 MUA model comparison =============================================
m=unique(tMUA.Model);
R2=[];
RTH = 0;

for i=1:length(m)
    R2 = [R2 tMUA.R2(strcmp(tMUA.Model,m{i}))];
end

v1=tMUA.Area(strcmp(tMUA.Model,m{1}))==1;
v4=tMUA.Area(strcmp(tMUA.Model,m{1}))==4;

sub = 'both';
if strcmp(sub,'M03') || strcmp(sub,'M04')
    v1=v1 & strcmp(tMUA.Monkey(strcmp(tMUA.Model,m{1})),sub);
    v4=v4 & strcmp(tMUA.Monkey(strcmp(tMUA.Model,m{1})),sub);
end

f6=figure;
msz=15;
set(f6,'Position',[100 100 1200 1200]);
for row=1:4
    for column=1:4
        subplot(4,4,((row-1)*4)+column); hold on;
        plot([0 100],[0 100],'k');
        
        XY = [R2(v1,row+1), R2(v1,column+1)];
        XYs = XY(XY(:,1)>RTH & XY(:,2)>RTH,:);
        scatter(XYs(:,1), XYs(:,2),msz,'Marker','o',...
            'MarkerEdgeColor',[.3 .3 .3],'MarkerFaceColor',[.3 .3 .3],...
            'MarkerFaceAlpha',0.5);
        XY = [R2(v4,row+1), R2(v4,column+1)];
        XYs = XY(XY(:,1)>RTH & XY(:,2)>RTH,:);
        scatter(XYs(:,1), XYs(:,2),msz,'Marker','o',...
            'MarkerEdgeColor',[.3 .8 .3],'MarkerFaceColor',[.3 .8 .3],...
            'MarkerFaceAlpha',0.5);
        set(gca, 'Box','off', 'xlim', [-2 100], 'ylim',[-2 100]);
        xlabel(m{row+1},'interpreter','none'); 
        ylabel(m{column+1},'interpreter','none');
        title('MUA');
        if row==1 && column==1
            legend({'','V1','V4'},'location','NorthWest');
        end
    end
end
if SaveFigs
    saveas(f6,fullfile(pngfld, 'EPHYS_MUA_ModelComparison.png'));
    saveas(f6,fullfile(svgfld, 'EPHYS_MUA_ModelComparison.svg'));
end
if CloseFigs; close(f6); end

%% Stats model comparison R2 MUA ==========================================
for RTH = [0 25]
    % stats V1
    MR2=R2(v1,2:5);
    sel=logical(sum(MR2>RTH,2));
    [p,tbl,stats] = kruskalwallis(MR2(sel,:),{'css','dog','p-lin','u-lin'});
    [c,m,h,gnames] = multcompare(stats);
    for i=1:size(c,1)
        fprintf(['V1 RTH-' num2str(RTH) ': ' gnames{c(i,1)} ' vs ' gnames{c(i,2)} ...
            ', p = ' num2str(c(i,6))  '\n'])
    end
    % stats V4
    MR2=R2(v4,2:5);
    sel=logical(sum(MR2>RTH,2));
    [p,tbl,stats] = kruskalwallis(MR2(sel,:),{'css','dog','p-lin','u-lin'});
    [c,m,h,gnames] = multcompare(stats);
    for i=1:size(c,1)
        fprintf(['V4 RTH-' num2str(RTH) ': ' gnames{c(i,1)} ' vs ' gnames{c(i,2)} ...
            ', p = ' num2str(c(i,6))  '\n'])
    end
end

%% SFig 7 & SFig 8: LFP model comparison ==================================
f7=figure; 
set(f7,'Position',[100 100 1600 1200]);
msz=15;

m=unique(tMUA.Model);
v1=tMUA.Area(strcmp(tMUA.Model,m{1}))==1;
v4=tMUA.Area(strcmp(tMUA.Model,m{1}))==4;

sub = 'both';
if strcmp(sub,'M03') || strcmp(sub,'M04')
    v1=v1 & strcmp(tMUA.Monkey(strcmp(tMUA.Model,m{1})),sub);
    v4=v4 & strcmp(tMUA.Monkey(strcmp(tMUA.Model,m{1})),sub);
end

% scatter dots
sig=unique(tLFP.SigType);
lfp_order = [3 1 2 5 4];
spn=1; fbn=1;

for fb=lfp_order
    lfpmods{fbn,1}=[];
    lfpmods{fbn,2}=[];

    m=unique(tLFP.Model);
    % reorder ---------
    m=m([3 4 1 2]);
    modname={'P-LIN','U-LIN','CSS','DoG'};
    % -----------------
    R2=[];
    for i=1:length(m)
        R2 = [R2 tLFP.R2(...
            strcmp(tLFP.Model,m{i}) & ...
            strcmp(tLFP.SigType,sig{fb}))];
    end
        
    for m1=1:4
        lfpmods{fbn,1}=[lfpmods{fbn,1} R2(v1,m1)];
        lfpmods{fbn,2}=[lfpmods{fbn,2} R2(v4,m1)];
        for m2=m1+1:4
            subplot(length(sig),6,spn); hold on;
            plot([-5 100],[-5 100],'k');

%             M1R2=R2(v1,m1); 
%             M2R2=R2(v1,m2);           
%             M1R2=M1R2(M1R2>0 & M2R2>0);
%             M2R2=M2R2(M1R2>0 & M2R2>0);
            
            scatter(R2(v1,m1), R2(v1,m2),msz,'Marker','o',...
                'MarkerEdgeColor',[.3 .3 .3],'MarkerEdgeAlpha',0,...
                'MarkerFaceColor',[.3 .3 .3],'MarkerFaceAlpha',0.25);
            scatter(R2(v4,m1), R2(v4,m2),msz,'Marker','o',...
                'MarkerEdgeColor',[.3 .8 .3],'MarkerEdgeAlpha',0,...
                'MarkerFaceColor',[.3 .8 .3],'MarkerFaceAlpha',0.25);

%             M1R2=R2(v4,m1); M1R2=M1R2(M1R2>0 & M2R2>0);
%             M2R2=R2(v4,m1); M1R2=M1R2(M1R2>0 & M2R2>0);          
                        
%             scatter(R2(v4,m1), R2(v4,m2),msz,'Marker','o',...
%                 'MarkerEdgeColor',[.3 .3 .3],'MarkerEdgeAlpha',0,...
%                 'MarkerFaceColor',[.3 .3 .3],'MarkerFaceAlpha',0.25);
            set(gca, 'Box','off', 'xlim', [-5 100], 'ylim',[-5 100],...
                'xticklabels',{},'yticklabels',{},'TickDir','out');
            xlabel(modname{m1},'interpreter','none');ylabel(modname{m2},'interpreter','none')
            title(sig{fb})
            if spn==1
                %legend({'','V1','V4'},'location','SouthEast');
            end
            spn=spn+1;
        end
    end
    fbn=fbn+1;
end
if SaveFigs
    saveas(f7,fullfile(pngfld, 'EPHYS_LFP_ModelComparison.png'));
    saveas(f7,fullfile(svgfld, 'EPHYS_LFP_ModelComparison.svg'));
end
if CloseFigs; close(f7); end

%% stats ==================================================================
FreqNames = {'Theta','Alpha','Beta','Gamma-low','Gamma-high'};

for fb=1:5
    % stats V1
    fprintf(['\n= V1 LFP-' FreqNames{fb} ' ================\n'])
    [p,tbl,stats] = kruskalwallis(lfpmods{fb,1},modname);
    fprintf(['H =' num2str(tbl{2,5}(1)) ', df = ' num2str(tbl{2,3}(1)) ', p = ' num2str(tbl{2,6}(1)) '\n'])
    [c,m,h,gnames] = multcompare(stats);
    for i=1:size(c,1)
        fprintf([gnames{c(i,1)} ' vs ' gnames{c(i,2)} ...
            ', p = ' num2str(c(i,6))  '\n'])
    end
    % stats V4
    fprintf(['\n= V4 LFP-' FreqNames{fb} ' ================\n'])
    [p,tbl,stats] = kruskalwallis(lfpmods{fb,2},modname);
    fprintf(['H =' num2str(tbl{2,5}(1)) ', df = ' num2str(tbl{2,3}(1)) ', p = ' num2str(tbl{2,6}(1)) '\n'])
    [c,m,h,gnames] = multcompare(stats);
    for i=1:size(c,1)
        fprintf([gnames{c(i,1)} ' vs ' gnames{c(i,2)} ...
            ', p = ' num2str(c(i,6))  '\n'])
    end
end

%% SFig 9: R2 for different ephys signals =================================
r2th=0;

ephys_MOD={'linear_ephys_cv1','linear_ephys_cv1_neggain',...
    'css_ephys_cv1','dog_ephys_cv1'};
ephys_MMS = MMS(:,1);

for m=1:length(ephys_MOD)
    RR=[];

    fprintf(['\n============ ' ephys_MOD{m} ' ===========\n'])
    
    model=ephys_MOD{m};
    s = strcmp(tMUA.Model,model);
    RR=[RR tMUA.R2(s)];
    
    s = strcmp(tLFP.Model,model);
    sig=unique(tLFP.SigType);
    lfp_order = [3 1 2 5 4];
    for i=lfp_order
        b=strcmp(tLFP.SigType,sig{i});
        RR=[RR tLFP.R2(s & b)];
    end
    LAB=['MUA';sig(lfp_order)];
    
    f8=figure; 
    set(f8,'Position',[100 100 1900 1350]);
    sgtitle(['R2 per Model: ' model],'interpreter','none');
    
    c=0;d=0;
    for ref=1:6
        c=c+1;
        for fb=1:6
            d=d+1;
            s=(RR(:,ref)>r2th & RR(:,fb)>r2th);
            subplot(6,6,d); hold on; 
            %scatter(RR(s,ref),RR(s,fb),120,[0.3 0.3 0.3],'Marker','.');
            binscatter(RR(s,ref),RR(s,fb),25,...
                'XLimits', [0 100],...
                'YLimits', [0 100],...
                'ShowEmptyBins', 'off')
            colorbar; 
            set(gca,'ColorScale','log');
            colormap(inferno)
            caxis([1 256]) 
            plot([0 100],[0 100],'Color',[.7 .7 .7],'LineWidth',2);
            xlabel(LAB{ref});ylabel(LAB{fb});
            %title(['R2 ' model],'interpreter','none');
            set(gca,'xlim',[0 100],'ylim',[0 100]);
            set(gca,'TickDir','out','xtick',[0 50 100],'ytick',[0 50 100]);
        end
    end
    if SaveFigs
        saveas(f8,fullfile(pngfld, ['EPHYS_MUA_R2_' ephys_MOD{m} '.png']));
        saveas(f8,fullfile(svgfld, ['EPHYS_MUA_R2_' ephys_MOD{m} '.svg']));
    end
    if CloseFigs; close(f8); end
    
    % Stats
    [p,tbl,stats] = kruskalwallis(RR,LAB');
    fprintf(['H =' num2str(tbl{2,5}(1)) ', df = ' num2str(tbl{2,3}(1)) ', p = ' num2str(tbl{2,6}(1)) '\n'])
    [c,m,h,gnames] = multcompare(stats);
    for i=1:size(c,1)
        fprintf([gnames{c(i,1)} ' (' num2str(m(c(i,1),1)) ' +/- ' num2str(m(c(i,1),2)) ') vs ' ...
            gnames{c(i,2)} ' (' num2str(m(c(i,2),1)) ' +/- ' num2str(m(c(i,2),2)) '), p = ' num2str(c(i,6))  '\n'])
    end    
    close all
end

%% Fig 11A,B: location and size across ephys channels =====================
RTH=25; SNRTH = 3;

% CSS -----
midx = 3; % 2 = U-LIN, 3 = CSS
fprintf(['MODEL ' ephys_MOD{midx} '\n']);
% 1 = MUA, 2 = ClasRF, 3 = Theta, 4 = Alpha, 5 = Beta, 6 = Gamma-low, 7 =
% Gamma-high

SigCompNames = {};
SigComp_nElec =[];
SigCompDist = [];
SigCompDist_columns = {'mean','std','median','iqr'};
SigCompSz = [];
SigCompSz_columns = {'mean_nS1','std_nS1','median_nS1','iqr_nS1',...
    'mean_nS2','std_nS2','median_nS2','iqr_nS2',...
    'mean_nS2/nS1','std_nS2/nS1','median_nS2/nS1','iqr_nS2/nS1'};

sub = 'both';

rown=1;
for sigidx1 = 1:7
    for sigidx2 = 1:7 
        SigCompNames{rown,1} = PRF_EST(midx,sigidx1).sig;
        SigCompNames{rown,2} = PRF_EST(midx,sigidx2).sig;
        
        % threshold
        if strcmp(sub,'M03') || strcmp(sub,'M04')
            elec = PRF_EST(midx,sigidx1).R2 > RTH & ...
                PRF_EST(midx,sigidx2).R2 > RTH & ...
                strcmp(PRF_EST(midx,sigidx2).M,sub);
        else
            elec = PRF_EST(midx,sigidx1).R2 > RTH & PRF_EST(midx,sigidx2).R2 > RTH;
        end
        
        % calculate distance
        DIST = sqrt((PRF_EST(midx,sigidx1).X(elec) - PRF_EST(midx,sigidx2).X(elec)).^2 + ...
            (PRF_EST(midx,sigidx1).Y(elec) - PRF_EST(midx,sigidx2).Y(elec)).^2);
        SigCompDist = [SigCompDist; nanmean(DIST) nanstd(DIST) nanmedian(DIST) iqr(DIST)];
        
        % calculate normalized size
        % normalize by MUA pRF
%         nS1 = PRF_EST(midx,sigidx1).S(elec)./PRF_EST(midx,1).S(elec);
%         nS2 = PRF_EST(midx,sigidx2).S(elec)./PRF_EST(midx,1).S(elec);
        % normalize by cRF
        nS1 = PRF_EST(midx,sigidx1).S(elec)./PRF_EST(midx,2).S(elec);
        nS2 = PRF_EST(midx,sigidx2).S(elec)./PRF_EST(midx,2).S(elec);
        
        SigCompSz = [SigCompSz;...
            nanmean(nS1(~isinf(nS1))) nanstd(nS1(~isinf(nS1))) ...
            nanmedian(nS1(~isinf(nS1))) iqr(nS1(~isinf(nS1))) ...
            nanmean(nS2(~isinf(nS2))) nanstd(nS2(~isinf(nS2))) ...
            nanmedian(nS2(~isinf(nS2))) iqr(nS2(~isinf(nS2))) ...
            nanmean(nS2(~isinf(nS1) & ~isinf(nS2))./nS1(~isinf(nS1) & ~isinf(nS2))) ...
            nanstd(nS2(~isinf(nS1) & ~isinf(nS2))./nS1(~isinf(nS1) & ~isinf(nS2))) ...
            nanmedian(nS2(~isinf(nS1) & ~isinf(nS2))./nS1(~isinf(nS1) & ~isinf(nS2))) ...
            iqr(nS2(~isinf(nS1) & ~isinf(nS2))./nS1(~isinf(nS1) & ~isinf(nS2)))];
        
        SigComp_nElec = [SigComp_nElec;sum(elec)];
        
        rown = rown+1;
    end
end

[sorted_names,sidx] = sortrows(SigCompNames);
sorted_n = SigComp_nElec(sidx,:);
sorted_dist = SigCompDist(sidx,:);
sorted_sz = SigCompSz(sidx,:);

uSig = unique(SigCompNames);
nSig = size(unique(SigCompNames),1);

DistMat_median = zeros(nSig);
DistMat_iqr = zeros(nSig);
DistMat_n = zeros(nSig);
SzMat_median = zeros(nSig);
SzMat_iqr = zeros(nSig);

for i=1:nSig
    ii=((i-1)*nSig)+1;
    DistMat_median(:,i) = sorted_dist(ii:ii+nSig-1,3);
    DistMat_iqr(:,i) = sorted_dist(ii:ii+nSig-1,4)./2;
    DistMat_n(:,i) = sorted_n(ii:ii+nSig-1);
    SzMat_median(:,i) = sorted_sz(ii:ii+nSig-1,11);
    SzMat_iqr(:,i) = sorted_sz(ii:ii+nSig-1,12)./2;
end

% re-order for plot
order=[4 3 5 1 2 7 6];
DistMat_median = DistMat_median(order,:); 
DistMat_median = DistMat_median(:,order); 
DistMat_iqr = DistMat_iqr(order,:); 
DistMat_iqr = DistMat_iqr(:,order); 
SzMat_median = SzMat_median(order,:); 
SzMat_median = SzMat_median(:,order); 
SzMat_iqr = SzMat_iqr(order,:); 
SzMat_iqr = SzMat_iqr(:,order); 
DistMat_n = DistMat_n(order,:);
DistMat_n = DistMat_n(:,order); 
uSig = uSig(order);

%
fcss=figure;
set(fcss,'Position',[10 10 2200 1200],'Renderer','painters');
colormap(viridis)
%colormap(brewermap([],'RdBu'));

subplot(2,3,1);
imagesc(DistMat_median)
set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,...
    'XTickLabelRotation',45)
caxis([0 1.5]);colorbar; 
title('pRF distance median');

subplot(2,3,2);
imagesc(DistMat_iqr)
set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,...
    'XTickLabelRotation',45)
caxis([0 0.8]);colorbar; 
title('pRF distance iqr');

subplot(2,3,3);
imagesc(DistMat_n)
set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,...
    'XTickLabelRotation',45,'ColorScale','log');
colorbar; 
caxis([1 1700]) 
title('n');

subplot(2,3,4);
imagesc(SzMat_median)
set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,...
    'XTickLabelRotation',45)
caxis([0 2]);colorbar; 
title('pRF relative sz median');

subplot(2,3,5);
imagesc(SzMat_iqr)
set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,...
    'XTickLabelRotation',45)
caxis([0 2]);colorbar; 
title('pRF relative sz iqr');

subplot(2,3,6);
nMask = DistMat_n >= 10;
imagesc(nMask)
set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,...
    'XTickLabelRotation',45);
colorbar; 
caxis([0 1]) 
title('mask n > 10');

sgtitle(['MODEL ' ephys_MOD{midx}], 'interpreter','none')

if SaveFigs
    saveas(fcss,fullfile(pngfld, 'EPHYS_xsig_dist.png'));
    saveas(fcss,fullfile(svgfld, 'EPHYS_xsig_dist.svg'));
end
if CloseFigs; close(fcss); end

fcss_sz=figure;hold on;
%errorbar(1:7,SzMat_median(:,2),SzMat_iqr(:,2),'ko','LineStyle','none'); % norm by MUA pRF
errorbar(1:7,SzMat_median(:,1),SzMat_iqr(:,1),'ko','LineStyle','none'); % norm by MUA cRF
ylabel('Normalized pRF size')
set(gca,'xlim',[0.5 7.5], 'xtick',1:7, 'xticklabels',...
    {'cRF-MUA','pRF-MUA','Theta','Alpha', 'Beta', 'Gam-low', 'Gam-high'},...
    'XTickLabelRotation',45);

if SaveFigs
    saveas(fcss_sz,fullfile(pngfld, 'EPHYS_xsig_sz.png'));
    saveas(fcss_sz,fullfile(svgfld, 'EPHYS_xsig_sz.svg'));
end
if CloseFigs; close(fcss_sz); end

%% Fig 10: Split Alpha and Beta LFP by positive negative gain =============

R2th = 20; % minimum R2
R2enh = 5; % R2 improvement
fb = {'Alpha','Beta'};

%% ALPHA Fig 10, SFig 11 ==================================================
fidx = 1;

DoG = tLFP(...
    strcmp(tLFP.Model,'dog_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:);
lin_n = tLFP(...
    strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,fb{fidx}),:);
lin = tLFP(...
    strcmp(tLFP.Model,'linear_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:);
css = tLFP(...
    strcmp(tLFP.Model,'css_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:);

lin_nGAM1 = tLFP(...
    strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,'lGamma'),:);
lin_nGAM2 = tLFP(...
    strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,'hGamma'),:);


% U-LIN ----
f_neg2 = figure; 
roi=1; % only do this for V1 channels
set(f_neg2,'Position',[10 10 900 800]);
chan_sel = lin_n.Area==roi & lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh;
chan_sel2 = lin_n.Area==roi & lin_n.R2>R2th;

chan_pgain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain>0;
chan_ngain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain<0;

chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th;
chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th;

subplot(2,2,1); hold on;
% gain alpha U-LIN
histogram(lin_n.gain(chan_sel2),-2000:50:2000,'FaceColor','k','FaceAlpha',0.5);
histogram(lin_n.gain(chan_ngain),-2000:50:2000,'FaceColor','r','FaceAlpha',0.5);
histogram(lin_n.gain(chan_pgain),-2000:50:2000,'FaceColor','b','FaceAlpha',0.5);

xlabel('gain U-LIN - ALL ELEC');ylabel('nChannels');
set(gca,'xlim',[-800 1800],'TickDir','out');
MM=median(lin_n.gain(chan_sel2));
yy=get(gca,'ylim');
plot([MM MM], [0 yy(2)+40],'k','Linewidth',5)
set(gca,'ylim',[0 yy(2)+10]);
title('Gain')

fprintf(['UNSELECTED - ALPHA MEDIAN GAIN: ' num2str(MM) ', IQR ' num2str(iqr(lin_n.gain(chan_sel))) '\n'])

subplot(2,2,2); hold on;
% gain alpha U-LIN
histogram(lin_n.gain(chan_sel),-1000:50:1000,'FaceColor','k','FaceAlpha',0.5);
xlabel('gain LIN-POSNEG');ylabel('nChannels');
set(gca,'xlim',[-800 1700],'TickDir','out');
MM=median(lin_n.gain(chan_sel));
yy=get(gca,'ylim');
plot([MM MM], [0 yy(2)+40],'k','Linewidth',5)
set(gca,'ylim',[0 yy(2)+10]);
title('Gain selected channels (based on R2 ULIN>PLIN')

fprintf(['ALPHA MEDIAN GAIN: ' num2str(MM) ', IQR ' num2str(iqr(lin_n.gain(chan_sel))) '\n'])

% Wilcoxon 1-tailed < 1
[p,h,stats] = signrank(lin_n.gain(chan_sel),0,'tail','left');
fprintf(['Gain < 0: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

subplot(2,2,3); hold on;
bb = [lin_n.ecc(chan_sel) lin.ecc(chan_sel)];

mEcc = [mean(lin_n.ecc(chan_ngain)) mean(lin_n.ecc(chan_pgain))];
sdEcc = [std(lin_n.ecc(chan_ngain)) std(lin_n.ecc(chan_pgain))];

mdEcc = [median(lin_n.ecc(chan_ngain)) median(lin_n.ecc(chan_pgain)) ];
iqrEcc = [iqr(lin_n.ecc(chan_ngain)) iqr(lin_n.ecc(chan_pgain)) ]./2;

errorbar(1:2,mdEcc,iqrEcc,...
    'ko','MarkerSize',4,'MarkerFaceColor','k','Linewidth',2)
set(gca,'xtick',1:2,'xticklabels',{'ULIN-N','ULIN-P'},...
    'xlim',[0.8 2.2],'TickDir','out')

ylabel('Eccentricity');
title('Ecc')

% Ecc difference
fprintf('-- STATS Ecc --\n');
fprintf(['mEcc nGain U-LIN: ' num2str(mEcc(1)) ', STD ' num2str(sdEcc(1)) '\n'])
fprintf(['mEcc pGain U-LIN: ' num2str(mEcc(2)) ', STD ' num2str(sdEcc(2)) '\n'])
fprintf(['mdEcc nGain U-LIN: ' num2str(mdEcc(1)) ', IQR ' num2str(iqrEcc(1)) '\n'])
fprintf(['mdEcc pGain U-LIN: ' num2str(mdEcc(2)) ', IQR ' num2str(iqrEcc(2)) '\n'])
[p,h,stats] = ranksum(lin_n.ecc(chan_ngain),lin_n.ecc(chan_pgain));
fprintf(['Ecc difference pos gain U-LIN vs neg gain U-LIN: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);


subplot(2,2,4); hold on;
bb = [lin_n.rfs(chan_sel) lin.rfs(chan_sel)];

mSz = [mean(lin_n.rfs(chan_ngain)) mean(lin_n.rfs(chan_pgain)) ];
sdSz = [std(lin_n.rfs(chan_ngain)) std(lin_n.rfs(chan_pgain)) ];

mdSz = [median(lin_n.rfs(chan_ngain)) median(lin_n.rfs(chan_pgain))];
iqrSz = [iqr(lin_n.rfs(chan_ngain)) iqr(lin_n.rfs(chan_pgain)) ]./2;

errorbar(1:2,mdSz,iqrSz,'ko','MarkerSize',4,'MarkerFaceColor','k','Linewidth',2)
set(gca,'xtick',1:2,'xticklabels',{'ULIN-N','ULIN-P'},...
    'xlim',[0.8 2.2],'TickDir','out')
ylabel('Size');
title('Size')

% Sz difference
fprintf('-- STATS Sz --\n');
fprintf(['mSz nGain U-LIN: ' num2str(mSz(1)) ', STD ' num2str(sdSz(1)) '\n'])
fprintf(['mSz pGain U-LIN: ' num2str(mSz(2)) ', STD ' num2str(sdSz(2)) '\n'])
fprintf(['mdSz nGain U-LIN: ' num2str(mdSz(1)) ', IQR ' num2str(iqrSz(1)) '\n'])
fprintf(['mdSz pGain U-LIN: ' num2str(mdSz(2)) ', IQR ' num2str(iqrSz(2)) '\n'])
[p,h,stats] = ranksum(lin_n.rfs(chan_ngain),lin_n.rfs(chan_pgain));
fprintf(['Sz difference pos gain U-LIN vs neg gain U-LIN: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);


if SaveFigs
    saveas(f_neg2,fullfile(pngfld, 'ALPHA_posneg_ecc_sz.png'));
    saveas(f_neg2,fullfile(svgfld, 'ALPHA_posneg_ecc_sz.svg'));
end
if CloseFigs; close(f_neg2); end

%% Compare eccentricities between alpha and gamma prfs ====================
fextra=figure;
set(fextra,'Position',[100 100 600 1500]);

dE=[];

subplot(4,2,1);hold on;
plot([0 15],[0 15]);
scatter(lin_n.ecc(chan_ngain),lin_nGAM1.ecc(chan_ngain));
set(gca,'xlim',[0 15],'ylim',[0 15]);
xlabel('A- CHAN: Ecc LOW GAM');ylabel('A- CHAN: Ecc ALPHA');

fprintf('Are Negative ALPHA pRF shifted relative to fixation from Low GAMMA?\n');
[p,h,stats] = signrank(lin_n.ecc(chan_ngain),lin_nGAM1.ecc(chan_ngain));
fprintf(['Median dEcc is ' ...
    num2str(median(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))) ...
    ' IQR ' ...
    num2str(iqr(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2) ...
    ' (POS --> towards fixation) \n']);
fprintf(['p = ' num2str(p) '\n']);
dE=[dE; ...
    median(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain)) ...
    iqr(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2];

subplot(4,2,2);hold on;
plot([0 15],[0 15]);
scatter(lin_n.ecc(chan_pgain),lin_nGAM1.ecc(chan_pgain));
set(gca,'xlim',[0 15],'ylim',[0 15]);
xlabel('A+ CHAN: Ecc LOW GAM');ylabel('A+ CHAN: Ecc ALPHA');

fprintf('Are Positive ALPHA pRF shifted relative to fixation from Low GAMMA?\n');
[p,h,stats] = signrank(lin_n.ecc(chan_pgain),lin_nGAM1.ecc(chan_pgain));
fprintf(['Median dEcc is ' ...
    num2str(median(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))) ...
    ' IQR ' ...
    num2str(iqr(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2) ...
    ' (POS --> towards fixation) \n']);
fprintf(['p = ' num2str(p) '\n']);
dE=[dE; ...
    median(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain)) ...
    iqr(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2];

subplot(4,2,3);hold on;
plot([0 15],[0 15]);
scatter(lin_n.ecc(chan_ngain),lin_nGAM2.ecc(chan_ngain));
set(gca,'xlim',[0 15],'ylim',[0 15]);
xlabel('A- CHAN: Ecc HIGH GAM');ylabel('A- CHAN: Ecc ALPHA');

fprintf('Are Negative ALPHA pRF shifted relative to fixation from High GAMMA?\n');
[p,h,stats] = signrank(lin_n.ecc(chan_ngain),lin_nGAM2.ecc(chan_ngain));
fprintf(['Median dEcc is ' ...
    num2str(median(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))) ...
    ' IQR ' ...
    num2str(iqr(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2) ...
    ' (POS --> towards fixation) \n']);
fprintf(['p = ' num2str(p) '\n']);
dE=[dE; ...
    median(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain)) ...
    iqr(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2];

subplot(4,2,4);hold on;
plot([0 15],[0 15]);
scatter(lin_n.ecc(chan_pgain),lin_nGAM2.ecc(chan_pgain));
set(gca,'xlim',[0 15],'ylim',[0 15]);
xlabel('A+ CHAN: Ecc HIGH GAM');ylabel('A+ CHAN: Ecc ALPHA');

fprintf('Are Positive ALPHA pRF shifted relative to fixation from High GAMMA?\n');
[p,h,stats] = signrank(lin_n.ecc(chan_pgain),lin_nGAM2.ecc(chan_pgain));
fprintf(['Median dEcc is ' ...
    num2str(median(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))) ...
    ' IQR ' ...
    num2str(iqr(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2) ...
    ' (POS --> towards fixation) \n']);
fprintf(['p = ' num2str(p) '\n']);
dE=[dE; ...
    median(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain)) ...
    iqr(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2];

subplot(4,2,5:8);hold on;
hold on;
bar(1:4, dE(:,1));
errorbar(1:4, dE(:,1), dE(:,2), 'k','Linestyle', 'none');
ylabel('dEcc Alpha vs Gam')
set(gca,'xtick',1:4,'xticklabels',{'A- v lG','A+ v lG','A- v hG','A+ v hG'})

%% plot the diff in within channel pRF location for low/high LFP freq =====
fshiftprf = figure;
set(fshiftprf,'Position',[100 100 1600 1600]);

Xa1 = lin_n.X(chan_ngain); Ya1 = lin_n.Y(chan_ngain);
Xlg1 = lin_nGAM1.X(chan_ngain); Ylg1 = lin_nGAM1.Y(chan_ngain);
Xa2 = lin_n.X(chan_pgain); Ya2 = lin_n.Y(chan_pgain);
Xlg2 = lin_nGAM1.X(chan_pgain); Ylg2 = lin_nGAM1.Y(chan_pgain);

% only prfs within 8 deg for display
csel1 = Xa1<8 & Ya1<2 & Xlg1<8 & Ylg1<2 & ...
    Xa1>-2 & Ya1>-8 & Xlg1>-2 & Ylg1>-8;
csel2 = Xa2<8 & Ya2<2 & Xlg2<8 & Ylg2<2 & ...
    Xa2>-2 & Ya2>-8 & Xlg2>-2 & Ylg2>-8;

subplot(2,2,1);hold on;
plot([0 0],[-10 10],'--','LineWidth',2,'Color',[0.5 0.5 0.5]);
plot([-10 10],[0 0],'--','LineWidth',2,'Color',[0.5 0.5 0.5]);

quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),...
    'Color',[0.8 0.8 0.8],...
    'AutoScale',0,'ShowArrowHead',0); % lines from gam to alpha pRF
scatter(Xa1(csel1),Ya1(csel1),20,'ko','MarkerFaceColor','r'); % alpha pRF center
scatter(Xlg1(csel1),Ylg1(csel1),20,'ko','MarkerFaceColor','b'); % lgam prf center
% quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),...
%     'k','LineWidth',1,'AutoScale',1,'ShowArrowHead',1); % arrows from gam to alpha pRF
plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g'); % fovea
set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,...
    'FontName','Myriad Pro','TickDir','out');
xlabel('Horizontal dva');ylabel('Vertical dva');

shiftsz = sqrt((Xa1(csel1)-Xlg1(csel1)).^2 + (Ya1(csel1)-Ylg1(csel1)).^2);
m_ss = mean(shiftsz); se_ss = std(shiftsz)./sqrt(length(shiftsz)); 
fprintf(['Average negALPHA shift size: ' num2str(m_ss) ' +/- ' num2str(se_ss) ' SEM\n']);

subplot(2,2,2);hold on;
quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),...
    'k','AutoScale',1,'ShowArrowHead',1)
plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g')
set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,...
    'FontName','Myriad Pro','TickDir','out');
xlabel('Horizontal dva');ylabel('Vertical dva');

subplot(2,2,3);hold on;
plot([0 0],[-10 10],'--','LineWidth',2,'Color',[0.5 0.5 0.5]);
plot([-10 10],[0 0],'--','LineWidth',2,'Color',[0.5 0.5 0.5]);

quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),...
    'Color',[0.8 0.8 0.8],...
    'AutoScale',0,'ShowArrowHead',0); % lines from gam to alpha pRF
scatter(Xa2(csel2),Ya2(csel2),20,'ko','MarkerFaceColor','r'); % alpha pRF center
scatter(Xlg2(csel2),Ylg2(csel2),20,'ko','MarkerFaceColor','b'); % lgam prf center
% quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),...
%     'k','LineWidth',1,'AutoScale',1,'ShowArrowHead',1); % arrows from gam to alpha pRF
plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g'); % fovea
set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,...
    'FontName','Myriad Pro','TickDir','out');
xlabel('Horizontal dva');ylabel('Vertical dva');

shiftsz = sqrt((Xa2(csel2)-Xlg2(csel2)).^2 + (Ya2(csel2)-Ylg2(csel2)).^2);
m_ss = mean(shiftsz); se_ss = std(shiftsz)./sqrt(length(shiftsz)); 
fprintf(['Average posALPHA shift size: ' num2str(m_ss) ' +/- ' num2str(se_ss) ' SEM\n']);

subplot(2,2,4);hold on;
quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),...
    'k','AutoScale',1,'ShowArrowHead',1)
plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g')
set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,...
    'FontName','Myriad Pro','TickDir','out');
xlabel('Horizontal dva');ylabel('Vertical dva');



if SaveFigs
    saveas(fextra,fullfile(pngfld, 'ALPHAvGAM_posneg_ecc_sz.png'));
    saveas(fextra,fullfile(svgfld, 'ALPHAvGAM_posneg_ecc_sz.svg'));
end
if CloseFigs; close(fextra); end

%% ALPHA vs GAMMA pRF separation index ====================================
lin_n = tLFP(...
    strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,'Alpha'),:);
DoG = tMUA(...
    strcmp(tMUA.Model,'dog_ephys_cv1'),:);

roi=1; % only do this for V1 channels
chan_sel = lin_n.Area==roi & lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh;
chan_sel2 = lin_n.Area==roi & lin_n.R2>R2th;

chan_pgain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain>0;
chan_ngain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain<0;

chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th;
chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th;
[mean(sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2)) ...
    std(sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2))./sqrt(sum(chan_ngain))]
[mean(sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2)) ...
    std(sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2))./sqrt(sum(chan_pgain))]

[mean(sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2)) ...
    std(sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM2.Y(chan_ngain)).^2))./sqrt(sum(chan_ngain))]
[mean(sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2)) ...
    std(sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2))./sqrt(sum(chan_pgain))]

fprintf('=== ALPHA ===\n')

% size
fprintf('Size diff (a-lG)\n')
fprintf('NGAIN\n')
[mean(lin_n.rfs(chan_ngain)-lin_nGAM1.rfs(chan_ngain)) ...
    std(lin_n.rfs(chan_ngain)-lin_nGAM1.rfs(chan_ngain))./sqrt(sum(chan_ngain))]
fprintf('PGAIN\n')
[mean(lin_n.rfs(chan_pgain)-lin_nGAM1.rfs(chan_pgain)) ...
    std(lin_n.rfs(chan_pgain)-lin_nGAM1.rfs(chan_pgain))./sqrt(sum(chan_pgain))]

% normalized size
fprintf('Norm Size diff (a/lG)\n')
fprintf('NGAIN\n')
[mean(lin_n.rfs(chan_ngain)./lin_nGAM1.rfs(chan_ngain)) ...
    std(lin_n.rfs(chan_ngain)./lin_nGAM1.rfs(chan_ngain))./sqrt(sum(chan_ngain))]
fprintf('PGAIN\n')
[mean(lin_n.rfs(chan_pgain)./lin_nGAM1.rfs(chan_pgain)) ...
    std(lin_n.rfs(chan_pgain)./lin_nGAM1.rfs(chan_pgain))./sqrt(sum(chan_pgain))]
fprintf('Norm Size diff (a/hG)\n')
fprintf('NGAIN\n')
[mean(lin_n.rfs(chan_ngain)./lin_nGAM2.rfs(chan_ngain)) ...
    std(lin_n.rfs(chan_ngain)./lin_nGAM2.rfs(chan_ngain))./sqrt(sum(chan_ngain))]
fprintf('PGAIN\n')
[mean(lin_n.rfs(chan_pgain)./lin_nGAM2.rfs(chan_pgain)) ...
    std(lin_n.rfs(chan_pgain)./lin_nGAM2.rfs(chan_pgain))./sqrt(sum(chan_pgain))]

% distance in relation to size
distn = sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2);
distp = sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2);
sz2n=lin_n.rfs(chan_ngain)+lin_nGAM1.rfs(chan_ngain);
sz2p=lin_n.rfs(chan_pgain)+lin_nGAM1.rfs(chan_pgain);


ff=figure;
set(ff,'Position',[10 10 600 1500]);

fprintf('DIST REL/Sz (dist/(s1+s2))\n')

subplot(4,2,1);hold on;
scatter(distn,sz2n);
scatter(distp,sz2p);
plot([0 20],[0 20]);
xlabel('distance A-G pRF');ylabel('sumsize A-G pRF')
set(gca,'xlim',[0 20],'ylim',[0 20]);
legend({'negative', 'positive'});

subplot(4,2,2);hold on;
fprintf('lG NEG\n')
fprintf('mean std median iqr\n')
[mean(distn./sz2n) std(distn./sz2n)./sqrt(sum(chan_ngain))...
    median(distn./sz2n) iqr(distn./sz2n)./2]
histogram(distn./sz2n,100,'Normalization','probability')
fprintf('lG POS\n')
[mean(distp./sz2p) std(distp./sz2p)./sqrt(sum(chan_pgain))...
    median(distp./sz2p) iqr(distp./sz2p)./2]
histogram(distp./sz2p,100,'Normalization','probability')
title('Low Gamma'); xlabel('Separation Index');

BC=[median(distn./sz2n) median(distp./sz2p)];
BCE=[iqr(distn./sz2n)/2 iqr(distp./sz2p)/2];

[p,h,stats] = ranksum(distn./sz2n, distp./sz2p);
fprintf('SEPARATION INDEX [Alpha-lGam] - Neg vs Pos channels\n');
fprintf(['p = ' num2str(p) '\n']);


distn = sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM2.Y(chan_ngain)).^2);
distp = sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2);
sz2n=lin_n.rfs(chan_ngain)+lin_nGAM2.rfs(chan_ngain);
sz2p=lin_n.rfs(chan_pgain)+lin_nGAM2.rfs(chan_pgain);

subplot(4,2,3);hold on;
scatter(distn,sz2n);
scatter(distp,sz2p);
plot([0 20],[0 20]);
xlabel('distance A-G pRF');ylabel('sumsize A-G pRF')
set(gca,'xlim',[0 20],'ylim',[0 20]);
legend({'negative', 'positive'});

subplot(4,2,4);hold on;
fprintf('hG NEG\n')
fprintf('mean std median iqr\n')
[mean(distn./sz2n) std(distn./sz2n)./sqrt(sum(chan_ngain))...
    median(distn./sz2n) iqr(distn./sz2n)./2]
histogram(distn./sz2n,100,'Normalization','probability')
fprintf('hG POS\n')
[mean(distp./sz2p) std(distp./sz2p)./sqrt(sum(chan_pgain))...
    median(distp./sz2p) iqr(distp./sz2p)./2]
histogram(distp./sz2p,100,'Normalization','probability')
title('High Gamma');  xlabel('Separation Index');

BC=[BC median(distn./sz2n) median(distp./sz2p)];
BCE=[BCE iqr(distn./sz2n)/2 iqr(distp./sz2p)/2];

[p,h,stats] = ranksum(distn./sz2n, distp./sz2p);
fprintf('SEPARATION INDEX [Alpha-hGam] - Neg vs Pos channels\n');
fprintf(['p = ' num2str(p) '\n']);

subplot(4,2,5:8); hold on;
bar(1:4, BC);
errorbar(1:4, BC, BCE, 'k', 'Linestyle','none')
ylabel('Separation Index')
set(gca, 'xtick',1:4,'xticklabels',{'lG-','lG+','hG-','hG+'})

if SaveFigs
    saveas(ff,fullfile(pngfld, 'ALPHAvGAM_sepidx.png'));
    saveas(ff,fullfile(svgfld, 'ALPHAvGAM_sepidx.svg'));
end
if CloseFigs; close(ff); end

%% ALPHA DoG ==============================================================
chan_sel = DoG.R2>R2th & DoG.R2>lin.R2+R2enh & ...
    DoG.normamp~=0 & DoG.ecc<16;

MM=median(DoG.normamp(chan_sel));
fprintf(['ALPHA MEDIAN NAMP: ' num2str(MM) ', IQR ' num2str(iqr(DoG.normamp(chan_sel))) '\n'])

% Wilcoxon 1-tailed < 1
[p,h,stats] = signrank(lin_n.gain(chan_sel),0,'tail','left');
fprintf(['Gain < 0: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

bb = [DoG.ecc(chan_sel) lin.ecc(chan_sel)];
% Wilcoxon 
[p,h,stats] = signrank(bb(:,1),bb(:,2));
fprintf(['ECC diff Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

MM=median(bb(:,2)-bb(:,1));
fprintf(['ALPHA MEDIAN ECC DIFF: ' num2str(median(lin.ecc(chan_sel)-DoG.ecc(chan_sel))) ...
    ', IQR ' num2str(iqr(lin.ecc(chan_sel)-DoG.ecc(chan_sel))) '\n'])

%% BETA SFig 10, 11 =======================================================
fidx = 2;

DoG = tLFP(...
    strcmp(tLFP.Model,'dog_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:);
lin_n = tLFP(...
    strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,fb{fidx}),:);
lin = tLFP(...
    strcmp(tLFP.Model,'linear_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:);
css = tLFP(...
    strcmp(tLFP.Model,'css_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:);

% U-LIN ----
f_neg2 = figure;
roi=1; % only do this for V1 channels
set(f_neg2,'Position',[10 10 900 1200]);
chan_sel = lin_n.Area==roi & lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh;
chan_sel2 = lin_n.Area==roi & lin_n.R2>R2th;

chan_pgain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain>0;
chan_ngain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain<0;

subplot(2,2,1); hold on;
% gain alpha U-LIN
histogram(lin_n.gain(chan_sel2),-2000:50:2000,'FaceColor','k','FaceAlpha',0.5);
histogram(lin_n.gain(chan_ngain),-2000:50:2000,'FaceColor','r','FaceAlpha',0.5);
histogram(lin_n.gain(chan_pgain),-2000:50:2000,'FaceColor','b','FaceAlpha',0.5);

xlabel('gain U-LIN - ALL ELEC');ylabel('nChannels');
set(gca,'xlim',[-800 1800],'TickDir','out');
MM=median(lin_n.gain(chan_sel2));
yy=get(gca,'ylim');
plot([MM MM], [0 yy(2)+40],'k','Linewidth',5)
set(gca,'ylim',[0 yy(2)+10]);
title('Gain')

fprintf(['UNSELECTED - BETA MEDIAN GAIN: ' num2str(MM) ', IQR ' num2str(iqr(lin_n.gain(chan_sel))) '\n'])

subplot(2,2,2); hold on;
% gain alpha U-LIN
histogram(lin_n.gain(chan_sel),-1000:50:1000,'FaceColor','k','FaceAlpha',0.5);
xlabel('gain LIN-POSNEG');ylabel('nChannels');
set(gca,'xlim',[-800 1700],'TickDir','out');
MM=median(lin_n.gain(chan_sel));
yy=get(gca,'ylim');
plot([MM MM], [0 yy(2)+40],'k','Linewidth',5)
set(gca,'ylim',[0 yy(2)+10]);
title('Gain selected channels (based on R2 ULIN>PLIN')

fprintf(['BETA MEDIAN GAIN: ' num2str(MM) ', IQR ' num2str(iqr(lin_n.gain(chan_sel))) '\n'])

% Wilcoxon 1-tailed < 1
[p,h,stats] = signrank(lin_n.gain(chan_sel),0,'tail','left');
fprintf(['Gain < 0: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

subplot(2,2,3); hold on;
bb = [lin_n.ecc(chan_sel) lin.ecc(chan_sel)];

mEcc = [mean(lin_n.ecc(chan_ngain)) mean(lin_n.ecc(chan_pgain))];
sdEcc = [std(lin_n.ecc(chan_ngain)) std(lin_n.ecc(chan_pgain))];

mdEcc = [median(lin_n.ecc(chan_ngain)) median(lin_n.ecc(chan_pgain))];
iqrEcc = [iqr(lin_n.ecc(chan_ngain)) iqr(lin_n.ecc(chan_pgain))]./2;

errorbar(1:2,mdEcc,iqrEcc,...
    'ko','MarkerSize',4,'MarkerFaceColor','k','Linewidth',2)
set(gca,'xtick',1:2,'xticklabels',{'ULIN-N','ULIN-P'},...
    'xlim',[0.8 2.2],'TickDir','out')

ylabel('Eccentricity');
title('Ecc')

% Ecc difference
fprintf('-- STATS Ecc --\n');
fprintf(['mEcc nGain U-LIN: ' num2str(mEcc(1)) ', STD ' num2str(sdEcc(1)) '\n'])
fprintf(['mEcc pGain U-LIN: ' num2str(mEcc(2)) ', STD ' num2str(sdEcc(2)) '\n'])
fprintf(['mdEcc nGain U-LIN: ' num2str(mdEcc(1)) ', IQR ' num2str(iqrEcc(1)) '\n'])
fprintf(['mdEcc pGain U-LIN: ' num2str(mdEcc(2)) ', IQR ' num2str(iqrEcc(2)) '\n'])
[p,h,stats] = ranksum(lin_n.ecc(chan_ngain),lin_n.ecc(chan_pgain));
fprintf(['Ecc difference pos gain U-LIN vs neg gain U-LIN: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

subplot(2,2,4); hold on;
bb = [lin_n.rfs(chan_sel) lin.rfs(chan_sel)];

mSz = [mean(lin_n.rfs(chan_ngain)) mean(lin_n.rfs(chan_pgain))];
sdSz = [std(lin_n.rfs(chan_ngain)) std(lin_n.rfs(chan_pgain))];

mdSz = [median(lin_n.rfs(chan_ngain)) median(lin_n.rfs(chan_pgain))];
iqrSz = [iqr(lin_n.rfs(chan_ngain)) iqr(lin_n.rfs(chan_pgain))]./2;

errorbar(1:2,mdSz,iqrSz,'ko','MarkerSize',4,'MarkerFaceColor','k','Linewidth',2)
set(gca,'xtick',1:4,'xticklabels',{'ULIN-N','ULIN-P'},...
    'xlim',[0.8 2.2],'TickDir','out')
ylabel('Size');
title('Size')


% Sz difference
fprintf('-- STATS Sz --\n');
fprintf(['mSz nGain U-LIN: ' num2str(mSz(1)) ', STD ' num2str(sdSz(1)) '\n'])
fprintf(['mSz pGain U-LIN: ' num2str(mSz(2)) ', STD ' num2str(sdSz(2)) '\n'])
fprintf(['mdSz nGain U-LIN: ' num2str(mdSz(1)) ', IQR ' num2str(iqrSz(1)) '\n'])
fprintf(['mdSz pGain U-LIN: ' num2str(mdSz(2)) ', IQR ' num2str(iqrSz(2)) '\n'])
[p,h,stats] = ranksum(lin_n.rfs(chan_ngain),lin_n.rfs(chan_pgain));
fprintf(['Sz difference pos gain U-LIN vs neg gain U-LIN: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

if SaveFigs
    saveas(f_neg2,fullfile(pngfld, 'BETA_posneg_ecc_sz.png'));
    saveas(f_neg2,fullfile(svgfld, 'BETA_posneg_ecc_sz.svg'));
end
if CloseFigs; close(f_neg2); end

%% Compare eccentricities between beta and gamma prfs =====================
lin_n = tLFP(...
    strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,'Beta'),:);

roi=1; % only do this for V1 channels
chan_sel = lin_n.Area==roi & lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh;
chan_sel2 = lin_n.Area==roi & lin_n.R2>R2th;

chan_pgain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain>0;
chan_ngain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain<0;

chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th;
chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th;

chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th;
chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th;

fextra=figure;
set(fextra,'Position',[10 10 600 1500]);

dE=[];

subplot(4,2,1);hold on;
plot([0 15],[0 15]);
scatter(lin_n.ecc(chan_ngain),lin_nGAM1.ecc(chan_ngain));
set(gca,'xlim',[0 15],'ylim',[0 15]);
xlabel('B- CHAN: Ecc LOW GAM');ylabel('B- CHAN: Ecc BETA');

fprintf('Are Negative BETA pRF shifted relative to fixation from Low GAMMA?\n');
[p,h,stats] = signrank(lin_n.ecc(chan_ngain),lin_nGAM1.ecc(chan_ngain));
fprintf(['Median dEcc is ' ...
    num2str(median(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))) ...
    ' IQR ' ...
    num2str(iqr(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2) ...
    ' (POS --> towards fixation) \n']);
fprintf(['p = ' num2str(p) '\n']);
dE=[dE; ...
    median(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain)) ...
    iqr(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2];


subplot(4,2,2);hold on;
plot([0 15],[0 15]);
scatter(lin_n.ecc(chan_pgain),lin_nGAM1.ecc(chan_pgain));
set(gca,'xlim',[0 15],'ylim',[0 15]);
xlabel('B+ CHAN: Ecc LOW GAM');ylabel('B+ CHAN: Ecc BETA');

fprintf('Are Positive BETA pRF shifted relative to fixation from Low GAMMA?\n');
[p,h,stats] = signrank(lin_n.ecc(chan_pgain),lin_nGAM1.ecc(chan_pgain));
fprintf(['Median dEcc is ' ...
    num2str(median(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))) ...
    ' IQR ' ...
    num2str(iqr(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2) ...
    ' (POS --> towards fixation) \n']);
fprintf(['p = ' num2str(p) '\n']);
dE=[dE; ...
    median(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain)) ...
    iqr(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2];

subplot(4,2,3);hold on;
plot([0 15],[0 15]);
scatter(lin_n.ecc(chan_ngain),lin_nGAM2.ecc(chan_ngain));
set(gca,'xlim',[0 15],'ylim',[0 15]);
xlabel('B- CHAN: Ecc HIGH GAM');ylabel('B- CHAN: Ecc BETA');

fprintf('Are Negative BETA pRF shifted relative to fixation from High GAMMA?\n');
[p,h,stats] = signrank(lin_n.ecc(chan_ngain),lin_nGAM2.ecc(chan_ngain));
fprintf(['Median dEcc is ' ...
    num2str(median(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))) ...
    ' IQR ' ...
    num2str(iqr(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2) ...
    ' (POS --> towards fixation) \n']);
fprintf(['p = ' num2str(p) '\n']);
dE=[dE; ...
    median(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain)) ...
    iqr(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2];

subplot(4,2,4);hold on;
plot([0 15],[0 15]);
scatter(lin_n.ecc(chan_pgain),lin_nGAM2.ecc(chan_pgain));
set(gca,'xlim',[0 15],'ylim',[0 15]);
xlabel('B+ CHAN: Ecc HIGH GAM');ylabel('B+ CHAN: Ecc BETA');

fprintf('Are Positive BETA pRF shifted relative to fixation from High GAMMA?\n');
[p,h,stats] = signrank(lin_n.ecc(chan_pgain),lin_nGAM2.ecc(chan_pgain));
fprintf(['Median dEcc is ' ...
    num2str(median(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))) ...
    ' IQR ' ...
    num2str(iqr(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2) ...
    ' (POS --> towards fixation) \n']);
fprintf(['p = ' num2str(p) '\n']);
dE=[dE; ...
    median(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain)) ...
    iqr(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2];

subplot(4,2,5:8);hold on;
hold on;
bar(1:4, dE(:,1));
errorbar(1:4, dE(:,1), dE(:,2), 'k','Linestyle', 'none');
ylabel('dEcc Beta vs Gam')
set(gca,'xtick',1:4,'xticklabels',{'B- v lG','B+ v lG','B- v hG','B+ v hG'})

%% plot the diff in within channel pRF location for low/high LFP freq =====
fshiftprf = figure;
set(fshiftprf,'Position',[100 100 1600 1600]);

Xa1 = lin_n.X(chan_ngain); Ya1 = lin_n.Y(chan_ngain);
Xlg1 = lin_nGAM1.X(chan_ngain); Ylg1 = lin_nGAM1.Y(chan_ngain);
Xa2 = lin_n.X(chan_pgain); Ya2 = lin_n.Y(chan_pgain);
Xlg2 = lin_nGAM1.X(chan_pgain); Ylg2 = lin_nGAM1.Y(chan_pgain);

% only prfs within 8 deg for display
csel1 = Xa1<8 & Ya1<2 & Xlg1<8 & Ylg1<2 & ...
    Xa1>-2 & Ya1>-8 & Xlg1>-2 & Ylg1>-8;
csel2 = Xa2<8 & Ya2<2 & Xlg2<8 & Ylg2<2 & ...
    Xa2>-2 & Ya2>-8 & Xlg2>-2 & Ylg2>-8;

subplot(2,2,1);hold on;
plot([0 0],[-10 10],'--','LineWidth',2,'Color',[0.5 0.5 0.5]);
plot([-10 10],[0 0],'--','LineWidth',2,'Color',[0.5 0.5 0.5]);

quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),...
    'Color',[0.8 0.8 0.8],...
    'AutoScale',0,'ShowArrowHead',0); % lines from gam to alpha pRF
scatter(Xa1(csel1),Ya1(csel1),20,'ko','MarkerFaceColor','r'); % alpha pRF center
scatter(Xlg1(csel1),Ylg1(csel1),20,'ko','MarkerFaceColor','b'); % lgam prf center
% quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),...
%     'k','LineWidth',1,'AutoScale',1,'ShowArrowHead',1); % arrows from gam to alpha pRF
plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g'); % fovea
set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,...
    'FontName','Myriad Pro','TickDir','out');
xlabel('Horizontal dva');ylabel('Vertical dva');

shiftsz = sqrt((Xa1(csel1)-Xlg1(csel1)).^2 + (Ya1(csel1)-Ylg1(csel1)).^2);
m_ss = mean(shiftsz); se_ss = std(shiftsz)./sqrt(length(shiftsz)); 
fprintf(['Average negALPHA shift size: ' num2str(m_ss) ' +/- ' num2str(se_ss) ' SEM\n']);

subplot(2,2,2);hold on;
quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),...
    'k','AutoScale',1,'ShowArrowHead',1)
plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g')
set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,...
    'FontName','Myriad Pro','TickDir','out');
xlabel('Horizontal dva');ylabel('Vertical dva');

subplot(2,2,3);hold on;
plot([0 0],[-10 10],'--','LineWidth',2,'Color',[0.5 0.5 0.5]);
plot([-10 10],[0 0],'--','LineWidth',2,'Color',[0.5 0.5 0.5]);

quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),...
    'Color',[0.8 0.8 0.8],...
    'AutoScale',0,'ShowArrowHead',0); % lines from gam to alpha pRF
scatter(Xa2(csel2),Ya2(csel2),20,'ko','MarkerFaceColor','r'); % alpha pRF center
scatter(Xlg2(csel2),Ylg2(csel2),20,'ko','MarkerFaceColor','b'); % lgam prf center
% quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),...
%     'k','LineWidth',1,'AutoScale',1,'ShowArrowHead',1); % arrows from gam to alpha pRF
plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g'); % fovea
set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,...
    'FontName','Myriad Pro','TickDir','out');
xlabel('Horizontal dva');ylabel('Vertical dva');

shiftsz = sqrt((Xa2(csel2)-Xlg2(csel2)).^2 + (Ya2(csel2)-Ylg2(csel2)).^2);
m_ss = mean(shiftsz); se_ss = std(shiftsz)./sqrt(length(shiftsz)); 
fprintf(['Average posALPHA shift size: ' num2str(m_ss) ' +/- ' num2str(se_ss) ' SEM\n']);

subplot(2,2,4);hold on;
quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),...
    'k','AutoScale',1,'ShowArrowHead',1)
plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g')
set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,...
    'FontName','Myriad Pro','TickDir','out');
xlabel('Horizontal dva');ylabel('Vertical dva');

if SaveFigs
    saveas(fextra,fullfile(pngfld, 'BETAvGAM_posneg_ecc_sz.png'));
    saveas(fextra,fullfile(svgfld, 'BETAvGAM_posneg_ecc_sz.svg'));
end
if CloseFigs; close(fextra); end

%% Beta vs Gamma pRF separation index =====================================
lin_n = tLFP(...
    strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,'Beta'),:);

roi=1; % only do this for V1 channels
chan_sel = lin_n.Area==roi & lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh;
chan_sel2 = lin_n.Area==roi & lin_n.R2>R2th;

chan_pgain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain>0;
chan_ngain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain<0;

chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th;
chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th;

chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th;
chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th;
[mean(sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2)) ...
    std(sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2))./sqrt(sum(chan_ngain))]
[mean(sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2)) ...
    std(sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2))./sqrt(sum(chan_pgain))]

[mean(sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2)) ...
    std(sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM2.Y(chan_ngain)).^2))./sqrt(sum(chan_ngain))]
[mean(sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2)) ...
    std(sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2))./sqrt(sum(chan_pgain))]

fprintf('=== BETA ===\n')

% size
fprintf('Size diff (a-lG)\n')
fprintf('NGAIN\n')
[mean(lin_n.rfs(chan_ngain)-lin_nGAM1.rfs(chan_ngain)) ...
    std(lin_n.rfs(chan_ngain)-lin_nGAM1.rfs(chan_ngain))./sqrt(sum(chan_ngain))]
fprintf('PGAIN\n')
[mean(lin_n.rfs(chan_pgain)-lin_nGAM1.rfs(chan_pgain)) ...
    std(lin_n.rfs(chan_pgain)-lin_nGAM1.rfs(chan_pgain))./sqrt(sum(chan_pgain))]

% normalized size
fprintf('Norm Size diff (a/lG)\n')
fprintf('NGAIN\n')
[mean(lin_n.rfs(chan_ngain)./lin_nGAM1.rfs(chan_ngain)) ...
    std(lin_n.rfs(chan_ngain)./lin_nGAM1.rfs(chan_ngain))./sqrt(sum(chan_ngain))]
fprintf('PGAIN\n')
[mean(lin_n.rfs(chan_pgain)./lin_nGAM1.rfs(chan_pgain)) ...
    std(lin_n.rfs(chan_pgain)./lin_nGAM1.rfs(chan_pgain))./sqrt(sum(chan_pgain))]
fprintf('Norm Size diff (a/hG)\n')
fprintf('NGAIN\n')
[mean(lin_n.rfs(chan_ngain)./lin_nGAM2.rfs(chan_ngain)) ...
    std(lin_n.rfs(chan_ngain)./lin_nGAM2.rfs(chan_ngain))./sqrt(sum(chan_ngain))]
fprintf('PGAIN\n')
[mean(lin_n.rfs(chan_pgain)./lin_nGAM2.rfs(chan_pgain)) ...
    std(lin_n.rfs(chan_pgain)./lin_nGAM2.rfs(chan_pgain))./sqrt(sum(chan_pgain))]

% distance in relation to size
distn = sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2);
distp = sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2);
sz2n=lin_n.rfs(chan_ngain)+lin_nGAM1.rfs(chan_ngain);
sz2p=lin_n.rfs(chan_pgain)+lin_nGAM1.rfs(chan_pgain);


fextra=figure;
set(fextra,'Position',[10 10 600 1500]);

subplot(4,2,1);hold on;
scatter(distn,sz2n);
scatter(distp,sz2p);
plot([0 20],[0 20]);
set(gca,'xlim',[0 20],'ylim',[0 20]);
xlabel('distance B-G pRF');ylabel('sumsize B-G pRF')
legend({'negative', 'positive'});

fprintf('DIST REL/Sz (dist/(s1+s2))\n')

subplot(4,2,2);hold on;
fprintf('lG NEG\n')
fprintf('mean std median iqr\n')
[mean(distn./sz2n) std(distn./sz2n)./sqrt(sum(chan_ngain))...
    median(distn./sz2n) iqr(distn./sz2n)./2]
histogram(distn./sz2n,100,'Normalization','probability')
fprintf('lG POS\n')
[mean(distp./sz2p) std(distp./sz2p)./sqrt(sum(chan_pgain))...
    median(distp./sz2p) iqr(distp./sz2p)./2]
histogram(distp./sz2p,100,'Normalization','probability')
title('Low Gamma'); xlabel('Separation Index');

BC=[median(distn./sz2n) median(distp./sz2p)];
BCE=[iqr(distn./sz2n)/2 iqr(distp./sz2p)/2];

distn = sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ...
    (lin_n.Y(chan_ngain)-lin_nGAM2.Y(chan_ngain)).^2);
distp = sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ...
    (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2);
sz2n=lin_n.rfs(chan_ngain)+lin_nGAM2.rfs(chan_ngain);
sz2p=lin_n.rfs(chan_pgain)+lin_nGAM2.rfs(chan_pgain);

subplot(4,2,3);hold on;
scatter(distn,sz2n);
scatter(distp,sz2p);
plot([0 20],[0 20]);
xlabel('distance B-G pRF');ylabel('sumsize B-G pRF')
set(gca,'xlim',[0 20],'ylim',[0 20]);
legend({'negative', 'positive'});

subplot(4,2,4);hold on;
fprintf('hG NEG\n')
fprintf('mean std median iqr\n')
[mean(distn./sz2n) std(distn./sz2n)./sqrt(sum(chan_ngain))...
    median(distn./sz2n) iqr(distn./sz2n)./2]
histogram(distn./sz2n,100,'Normalization','probability')
fprintf('hG POS\n')
[mean(distp./sz2p) std(distp./sz2p)./sqrt(sum(chan_pgain))...
    median(distp./sz2p) iqr(distp./sz2p)./2]
histogram(distp./sz2p,100,'Normalization','probability')
title('High Gamma');  xlabel('Separation Index'); 

BC=[BC median(distn./sz2n) median(distp./sz2p)];
BCE=[BCE iqr(distn./sz2n)/2 iqr(distp./sz2p)/2];

subplot(4,2,5:8); hold on;
bar(1:4, BC);
ylabel('Separation Index')
errorbar(1:4, BC, BCE, 'k', 'Linestyle','none')
set(gca, 'xtick',1:4,'xticklabels',{'lG-','lG+','hG-','hG+'})

if SaveFigs
    saveas(fextra,fullfile(pngfld, 'BETAvGAM_sepidx.png'));
    saveas(fextra,fullfile(svgfld, 'BETAvGAM_sepidx.svg'));
end
if CloseFigs; close(fextra); end

%% Beta DoG ===============================================================
chan_sel = DoG.R2>R2th & DoG.R2>lin.R2+R2enh & ...
    DoG.normamp~=0 & DoG.ecc<16;

MM=median(DoG.normamp(chan_sel));
fprintf(['BETA MEDIAN NAMP: ' num2str(MM) ', IQR ' num2str(iqr(DoG.normamp(chan_sel))) '\n'])

% Wilcoxon 1-tailed < 1
[p,h,stats] = signrank(lin_n.gain(chan_sel),0,'tail','left');
fprintf(['Gain < 0: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

bb = [DoG.ecc(chan_sel) lin.ecc(chan_sel)];
% Wilcoxon 
[p,h,stats] = signrank(bb(:,1),bb(:,2));
fprintf(['ECC diff Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

MM=median(bb(:,2)-bb(:,1));
fprintf(['BETA MEDIAN ECC DIFF: ' num2str(median(lin.ecc(chan_sel)-DoG.ecc(chan_sel))) ...
    ', IQR ' num2str(iqr(lin.ecc(chan_sel)-DoG.ecc(chan_sel))) '\n'])



%% ========================================================================
%  FMRI-EPHYS PRF ANALYSIS ------------------------------------------------
%  ========================================================================

%% Fig 12 & SFig 12: Correlate MRI-ephys REGRESSION BASED =================
% based on the ecc vs size correlation
% PER SIGNAL SOURCE
% - calculate linear regression on random (filtered) subset of data
% - repeat
% - filter results
% ACROSS SIGNAL SOURCES
% - correlate the bootstrapped fit-result parameters
% - optional >> do this in a bootstrapped way
rng(1); % seed the random number generator

% adjust these to check robustness ----------------------------------------
Rth_mri = 5; % R2 threshold MRI
Rth_ephys = 50; % R2 threshold ephys
selVFC = 1; % if true, only use lower right visual quadrant
selANAT = 0; % if true, only use approximate area of electrodes
ephysSUB = 'Both';

% -------------------------------------------------------------------------

MaxECC = [10 30]; % max ecc to use for fitting [V1 V4]
MODS = {...
    'linhrf_cv1_mhrf','linear_ephys_cv1';...
    'linhrf_cv1_mhrf_neggain','linear_ephys_cv1_neggain';...
    'csshrf_cv1_mhrf','css_ephys_cv1';...
    'doghrf_cv1_mhrf','dog_ephys_cv1';...
    };
MRI_MODEL = MODS(:,1); EPHYS_MODEL = MODS(:,2);
MMS={'linear','linear_ng','css','dog'};

warning off;
cmROI = {'V1','V4'};

fprintf('=======================\n');
for m = 3 % select which model
    fprintf(['\nCrossmodal Correlation for Model: ' MODS{m} '\n']);
    s_R2 = T(modidx.(MRI_MODEL{m})).mod.R2 > Rth_mri & ...
        T(modidx.(MRI_MODEL{m})).mod.rfs < 1000;
    
    % collect mri prfs
    for r = 1:size(cmROI,2)
        if ~selANAT
            SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{m})).mod.ROI,...
                ck_GetROIidx(cmROI(r),rois) );
        else
            if strcmp(cmROI{r},'V1')
                SSS = s_R2 & T(modidx.(MRI_MODEL{m})).mod.ELEC_V1 > 0 & ...
                    ismember( T(modidx.(MRI_MODELS{m})).mod.ROI,...
                ck_GetROIidx(cmROI(r),rois) );
            elseif strcmp(cmROI{r},'V4')
                SSS = s_R2 & T(modidx.(MRI_MODEL{m})).mod.ELEC_V4 > 0 & ...
                    ismember( T(modidx.(MRI_MODELS{m})).mod.ROI,...
                ck_GetROIidx(cmROI(r),rois) );
            end
        end
        if selVFC
            VFSEL = T(modidx.(MRI_MODEL{m})).mod.X >= 0 & ...
                T(modidx.(MRI_MODEL{m})).mod.Y <=0;
            SSS = SSS & VFSEL;
        end
        
%         if selANAT && strcmp(cmROI{r},'V1')
%             ANATSEL = T(modidx.(MRI_MODEL{m})).mod.ELEC_V1 > 0; 
%             SSS = SSS & ANATSEL;
%         elseif selANAT && strcmp(cmROI{r},'V4')
%             ANATSEL = T(modidx.(MRI_MODEL{m})).mod.ELEC_V4 > 0; 
%             SSS = SSS & ANATSEL;
%         end

        if strcmp(cmROI{r},'V1') % V1
            mri1(m).ECC = T(modidx.(MRI_MODEL{m})).mod.ecc(SSS);
            mri1(m).S = T(modidx.(MRI_MODEL{m})).mod.rfs(SSS);
            mri1(m).G = T(modidx.(MRI_MODEL{m})).mod.gain(SSS);
        elseif strcmp(cmROI{r},'V4') % V4
            mri4(m).ECC = T(modidx.(MRI_MODEL{m})).mod.ecc(SSS);
            mri4(m).S = T(modidx.(MRI_MODEL{m})).mod.rfs(SSS);
            mri4(m).G = T(modidx.(MRI_MODEL{m})).mod.gain(SSS);
        end
    end
    
    % collect ephys prfs
    % MUA V1
    if strcmp(ephysSUB,'M03') || strcmp(ephysSUB,'M04')
        s = strcmp(tMUA.Model,EPHYS_MODEL{m}) & ...
             strcmp(tMUA.Monkey,ephysSUB) & ...
            tMUA.Area == 1 & tMUA.R2 > Rth_ephys & tMUA.rfs < 1000;
    else 
        s = strcmp(tMUA.Model,EPHYS_MODEL{m}) & ...
            tMUA.Area == 1 & tMUA.R2 > Rth_ephys & tMUA.rfs < 1000;
    end
    % Exclude 2 outlier V1 arrays in both monkeys
    s2 = ( strcmp(tMUA.Monkey,'M03') & (tMUA.Array == 11 | tMUA.Array == 13)) | ...
        ( strcmp(tMUA.Monkey,'M04') & (tMUA.Array == 10 | tMUA.Array == 12));
    s(s2) = false;   
    mua1(m).ECC = tMUA.ecc(s);
    mua1(m).S = tMUA.rfs(s);
    mua1(m).G = tMUA.gain(s);
    % MUA V4
    s = strcmp(tMUA.Model,EPHYS_MODEL{m}) & ...
        tMUA.Area == 4 & tMUA.R2 > Rth_ephys & tMUA.rfs < 1000;
    mua4(m).ECC = tMUA.ecc(s);
    mua4(m).S = tMUA.rfs(s);
    mua4(m).G = tMUA.gain(s);
    
    % LFP
    freqband=unique(tLFP.SigType);
    for fb = 1: length(freqband)
        % V1
        if strcmp(ephysSUB,'M03') || strcmp(ephysSUB,'M04')
            s = strcmp(tLFP.Model,EPHYS_MODEL{m}) & ...
                strcmp(tLFP.Monkey,ephysSUB) & ...
                tLFP.Area == 1 & tLFP.R2 > Rth_ephys & tLFP.rfs < 1000 & ...
                strcmp(tLFP.SigType, freqband{fb});
        else
            s = strcmp(tLFP.Model,EPHYS_MODEL{m}) & ...
                tLFP.Area == 1 & tLFP.R2 > Rth_ephys & tLFP.rfs < 1000 & ...
                strcmp(tLFP.SigType, freqband{fb});
        end
        s2 = ( strcmp(tMUA.Monkey,'M03') & (tMUA.Array == 11 | tMUA.Array == 13)) | ...
            ( strcmp(tMUA.Monkey,'M04') & (tMUA.Array == 10 | tMUA.Array == 12));
        s(s2) = false;
        
        lfp1(fb,m).freqband =  freqband{fb};
        lfp1(fb,m).ECC =  tLFP.ecc(s);
        lfp1(fb,m).S =  tLFP.rfs(s);
        lfp1(fb,m).G =  tLFP.gain(s);
        
        % V4
        s = strcmp(tLFP.Model,EPHYS_MODEL{m}) & ...
            tLFP.Area == 4 & tLFP.R2 > Rth_ephys & tLFP.rfs < 1000 & ...
            strcmp(tLFP.SigType, freqband{fb});
        lfp4(fb,m).freqband =  freqband{fb};
        lfp4(fb,m).ECC =  tLFP.ecc(s);
        lfp4(fb,m).S =  tLFP.rfs(s);
        lfp4(fb,m).G =  tLFP.gain(s);
    end
    
    % Calculate  & bootstrap linear regressions
    % =============================================
    fprintf('Performing regressions on ALL data\n')
    MRI_stats1 = regstats(mri1(m).S,mri1(m).ECC);
    MRI_stats4 = regstats(mri4(m).S,mri4(m).ECC);
    MUA_stats1 = regstats(mua1(m).S,mua1(m).ECC);
    MUA_stats4 = regstats(mua4(m).S,mua4(m).ECC);
    
    for fb=1:length(freqband)
        try
            LFP_stats1{fb} = regstats(lfp1(fb,m).S,lfp1(fb,m).ECC);
        catch
            LFP_stats1{fb} = [];
        end
        try
            LFP_stats4{fb} = regstats(lfp4(fb,m).S,lfp4(fb,m).ECC);
        catch
            LFP_stats4{fb} = [];
        end
    end
    
    % Do statistics on model comparisons with fitlm
    ck_xmod_stats;
    XMOD(m).model = MMS{m}; XMOD(m).stats = stats;
    
    if m==2
        ck_xmod_stats_lowfreqlfp;
    end
end
warning on;

%% Fig 11C: Cross-signal comparison exponential parameter =================
RTHRES = 25;
clear exptvals_MUA exptvals_LFP exptvals_MRI kw
rois2 = [1 4];
f=figure;set(f,'Position',[10 10 1200 500]);
ephysSub='M04';

for r=1:length(rois2)
    kw{r}=[];
    exptvals_MRI{r} = [exptv(1).roi{rois2(r)};exptv(2).roi{rois2(r)}];  
    kw{r}=[kw{r}; exptvals_MRI{r} ones(size(exptvals_MRI{r}))];
    
    if strcmp(ephysSub,'M03') || strcmp(ephysSub,'M04')
        exptvals_MUA{r} = tMUA.expt(strcmp(tMUA.Model,'css_ephys_cv1') & ...
            strcmp(tMUA.Monkey,ephysSub) & strcmp(tMUA.SigType,'MUA') & ...
            tMUA.Area==rois2(r) & tMUA.R2 > RTHRES );
        fprintf(['Only ' ephysSub '\n'])
    else
        exptvals_MUA{r} = tMUA.expt(strcmp(tMUA.Model,'css_ephys_cv1') & ...
            strcmp(tMUA.SigType,'MUA') & tMUA.Area==rois2(r) & tMUA.R2 > RTHRES );
    end
    
    kw{r}=[kw{r}; exptvals_MUA{r} 2*ones(size(exptvals_MUA{r}))];
    sig=unique(tLFP.SigType);
    lfp_order = [3 1 2 5 4]; spn=1;
    
    cl=2;
    for fb=lfp_order
        cl=cl+1;
        if strcmp(ephysSub,'M03') || strcmp(ephysSub,'M04')
            exptvals_LFP{r,spn} = tLFP.expt(strcmp(tLFP.Model,'css_ephys_cv1') & ...
                strcmp(tLFP.Monkey,ephysSub) & strcmp(tLFP.SigType,sig{fb}) & ...
                tLFP.Area==rois2(r) & tLFP.R2 > RTHRES );
            fprintf(['Only ' ephysSub '\n'])
        else
            exptvals_LFP{r,spn} = tLFP.expt(strcmp(tLFP.Model,'css_ephys_cv1') & ...
                strcmp(tLFP.SigType,sig{fb}) & tLFP.Area==rois2(r) & tLFP.R2 > RTHRES );
        end
        if size(exptvals_LFP{r,spn},1)>3
            kw{r}=[kw{r}; exptvals_LFP{r,spn} cl*ones(size(exptvals_LFP{r,spn}))];
        end
        spn=spn+1;
    end
    
    figure(f);
    subplot(1,2,r); hold on;
    bar(1:7,[mean(exptvals_MRI{r}) ...
        mean(exptvals_MUA{r}) ...
        mean(exptvals_LFP{r,1}) ...
        mean(exptvals_LFP{r,2}) ...
        mean(exptvals_LFP{r,3}) ...
        mean(exptvals_LFP{r,4}) ...
        mean(exptvals_LFP{r,5})]);
    errorbar(1:7,[mean(exptvals_MRI{r}) ...
        mean(exptvals_MUA{r}) ...
        mean(exptvals_LFP{r,1}) ...
        mean(exptvals_LFP{r,2}) ...
        mean(exptvals_LFP{r,3}) ...
        mean(exptvals_LFP{r,4}) ...
        mean(exptvals_LFP{r,5})],...
        [std(exptvals_MRI{r}) ...
        std(exptvals_MUA{r}) ...
        std(exptvals_LFP{r,1}) ...
        std(exptvals_LFP{r,2}) ...
        std(exptvals_LFP{r,3}) ...
        std(exptvals_LFP{r,4}) ...
        std(exptvals_LFP{r,5})],'ko','Linestyle','none')
    set(gca,'xlim',[.5 7.5],'ylim',[0 1.1],'xtick',1:7,'xticklabels',...
        {'MRI','MUA', 'THETA','ALPHA','BETA','lGAM','hGAM'})
    title(['V' num2str(rois2(r))])
    ylabel('pRF exponent')
    
    
    fprintf('-- Compare < 1 --\n');
    fprintf(['AREA V' num2str(rois2(r)) '\n']);
    
    [p,h,stats] = signrank(exptvals_MRI{r},1,'tail','left');
    fprintf(['MRI, Wilcoxon EXPT < 1 : z = ' num2str(stats.zval) ', p = ' num2str(p) '\n']);
    [p,h,stats] = signrank(exptvals_MUA{r},1,'tail','left');
    fprintf(['MUA, Wilcoxon EXPT < 1 : z = ' num2str(stats.zval) ', p = ' num2str(p) '\n']);
    for i=1:5
        if size(exptvals_LFP{r,i},1) > 5
            [p,h,stats] = signrank(exptvals_LFP{r,i},1,'tail','left');
            if isfield(stats,'zval')
                fprintf(['LFP-' num2str(i) ', Wilcoxon EXPT < 1 : z = ' num2str(stats.zval) ', p = ' num2str(p) '\n']);
            end
        end
    end
    
    fprintf('-- Compare across signals --\n');
    fprintf(['AREA V' num2str(rois2(r)) '\n']);
    [expcss(r).p,expcss(r).tbl,expcss(r).stats] = kruskalwallis(kw{r}(:,1), kw{r}(:,2));
    [expcss(r).c,expcss(r).m,expcss(r).h,expcss(r).gnames] = multcompare(expcss(r).stats);
end

if SaveFigs
    saveas(f,fullfile(pngfld, 'xsig_prfexp.png'));
    saveas(f,fullfile(svgfld, 'xsig_prfexp.svg'));
end
if CloseFigs; close(f); end



%% ========================================================================
%  CLEAN UP ---------------------------------------------------------------
%  ========================================================================

%% clear and close  
close all; clc;