V1-V4-LFPs-and-Visual-Attention/Code/process_cGC_within_V1V4_MultiMethods_PreFirstDim.m

addpath(genpath('./support_routines'));
addpath(genpath('./support_tools'));

clear all; close all; clc;

sSubjectName='monkey 2';
sVisualArea='V4';

reComputeGCs=0;
reComputeGCsShuffling=0;

plotSinglePenGCs=0;
plotPenAvgGCs=1;
reportSignifThreshold=1;
plotPenAvgGCsSingChRes=0;

lamLayerTag={'Supra','Granr','Infra'};
gratCondTag={'RF','OUTrp','OUT1','OUT2'};
spectWinTag={'Full','Alpha','Beta','Gamma'};

%% cGC implementation settings
numInformChs=-2;
% numInformChs = [-Inf -2 -1 0 1 2 Inf];
% number of informative channels used as regressors in the cGC strategy
%   0 = Unconditional GC
%  >0 = cGC to most informative channels
%  <0 = cGC to m. inf. outside laminar compartments of each X,Y channel pair
% Inf = cGC to ALL channles (ignores most informative computation)
compMethod='Mvgc';
% compMethod = 'Mvgc','Geweke','MatPart'; choose method to compute cGC
%   Mvgc  = (Barnett & Seth, J. Neurosci. Methods, 2014)
%  Geweke = (Geweke, J. of the American Stat. Assoc., 1984)
% MatPart = (Chen et al., J. Neurosci Methods 2006)

%% DATA SAMPLING SETTINGS
SAMPLING_FREQ=1017.375;
NUM_TIME_SAMPLES_EXTRACTED=1024;
timeAxis0=linspace(-NUM_TIME_SAMPLES_EXTRACTED/SAMPLING_FREQ,0,NUM_TIME_SAMPLES_EXTRACTED)*1e3;
NUM_TIME_SAMPLES=512;
timeAxisPreFirstDim=timeAxis0(NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end);
freqAxisPreFirstDim=linspace(0,SAMPLING_FREQ/2,NUM_TIME_SAMPLES);
DOWNSAMPLING_FACTOR=4;
NUM_TIME_DWSAMPLES=NUM_TIME_SAMPLES/DOWNSAMPLING_FACTOR;
DWSAMPLING_FREQ=SAMPLING_FREQ/DOWNSAMPLING_FACTOR;
timeAxisPreFirstDimDwSamp=linspace(timeAxisPreFirstDim(1),timeAxisPreFirstDim(end),NUM_TIME_DWSAMPLES);
freqAxisPreFirstDimDwSamp=linspace(0,DWSAMPLING_FREQ/2,NUM_TIME_DWSAMPLES+1);
freqIndices2Plot=find(freqAxisPreFirstDimDwSamp>2);


if reComputeGCs
    subjectInfos=getSubjectInfos(sSubjectName);
    subjectData=getSubjectData(sSubjectName,sVisualArea,'ALL','LFPs','ALL');
    
    % REMOVE DATA IN TIME WINDOWS OTHER THAN PREFIRSTDIM TO SAVE MEMORY
    for pp=1:length(subjectData.LfpStruct)
        if ~isempty(subjectData.LfpStruct(pp).Sorted)
            rmfn=fieldnames(subjectData.LfpStruct(pp).Sorted.Full.Supra.RF);
            rmfn(cellfun(@(fnm) strcmpi(fnm,'PreFirstDimDataBiZSc'),rmfn))=[];
            for ll=1:3
                subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF=rmfield(subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF,rmfn);
                subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1=rmfield(subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1,rmfn);
                subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2=rmfield(subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2,rmfn);
            end
        end
    end
    
    %% RANDPERM FOR AVG SUBSAMPLING OF TRIALS AT OUTSIDE LOCATIONS
    rng(0); % will pick exactly the same trials at every run,
    % Note: varying this seed might give slightly different cGC magnitudes
    % (different trials are selected) but results are robust to permutation
    for pp=1:length(subjectData.LfpStruct)
        if ~isempty(subjectData.LfpStruct(pp).Sorted)
            currNtrOUT1=size(subjectData.LfpStruct(pp).Sorted.Full.Supra.OUT1.PreFirstDimDataBiZSc,2);
            currRandPermTrials=randperm(2*currNtrOUT1,currNtrOUT1);
            for ll=1:3
                currOUTcat=cat(2,subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1.PreFirstDimDataBiZSc,...
                    subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2.PreFirstDimDataBiZSc);
                subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUTrp.PreFirstDimDataBiZSc=currOUTcat(:,currRandPermTrials,:);
            end
        end
    end
    eval(clearoutmem('curr*'));
    
    %% Information Toolbox settings (support_tools) - DOI: 10.1186/1471-2202-10-81
    numQuantBins=4;
    MImethod='dr';
    MIbias='qe';
    MIbtsp=0;
    bandWinsWidth=2;
    bandWinsStart=.01; % Starts at 1 Hz
    bandWinsEnd=100;   % Stops at ~100 Hz
    bandWinsEnd=min(bandWinsEnd,SAMPLING_FREQ/2);
    NUM_BAND_WINS=floor((bandWinsEnd-bandWinsStart)/bandWinsWidth);
    bandWinsCutoffs=[linspace(bandWinsStart,bandWinsEnd-bandWinsWidth,NUM_BAND_WINS);...
        linspace(bandWinsStart,bandWinsEnd-bandWinsWidth,NUM_BAND_WINS)+bandWinsWidth]';
    freqAxisBandWins=mean(bandWinsCutoffs,2);
    filterOrder=2;
    filterBCoeffs=nan(NUM_BAND_WINS,2*filterOrder+1);
    filterACoeffs=nan(NUM_BAND_WINS,2*filterOrder+1);
    filterFreqResp=nan(NUM_BAND_WINS,512);
    filterImpResp=nan(NUM_BAND_WINS,2048);
    filterFreqAxis=nan(1,512);
    filterTimeAxis=nan(1,512);
    for bw=1:NUM_BAND_WINS
        [filterBCoeffs(bw,:),filterACoeffs(bw,:)]=butter(filterOrder,bandWinsCutoffs(bw,:)./(SAMPLING_FREQ/2),'bandpass');
        [filterFreqResp(bw,:),filterFreqAxis]=freqz(filterBCoeffs(bw,:),filterACoeffs(bw,:),512,SAMPLING_FREQ);
        [filterImpResp(bw,:),filterTimeAxis]=impz(filterBCoeffs(bw,:),filterACoeffs(bw,:),2048,SAMPLING_FREQ);
    end
    fresEstim=1/4*NUM_TIME_DWSAMPLES;
    lambdafreqAxis=linspace(0,freqAxisPreFirstDimDwSamp(end),fresEstim+1);
    
    %% MVGC Toolbox settings (support_tools) - DOI: 10.1016/j.jneumeth.2013.10.018
    regmode='LWR';
    morder=10;     % set up by optimization of AIC and R² (see Suppl. Figure S7)
    acmaxlags=[];  % autocovariance max lags
    acdectol=[];   % autocovariance decay tolerance (sets memory of the process)
    NShuffles=100; % number of shuffles used to compute significance threshold
    
    parfor pp=1:length(subjectData.LfpStruct)
        if ~isempty(subjectData.LfpStruct(pp).Sorted)
            for cnd=1:2
                disp(['Processing cGC for ' sSubjectName ' ' sVisualArea ' PEN ' num2str(subjectData.LfpStruct(pp).penNum) ', ' gratCondTag{cnd}]);
                
                currX0=[subjectData.LfpStruct(pp).Sorted.Full.Supra.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end);...
                    subjectData.LfpStruct(pp).Sorted.Full.Granr.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end);...
                    subjectData.LfpStruct(pp).Sorted.Full.Infra.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end)];
                currXp=permute(currX0,[1 3 2]);
                %eval(clearoutmem('currX0'));
                
                currSchs=1:size(subjectData.LfpStruct(pp).Sorted.Full.Supra.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1);
                currGchs=max([currSchs 0])+(1:size(subjectData.LfpStruct(pp).Sorted.Full.Granr.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1));
                currIchs=max([currSchs currGchs 0])+(1:size(subjectData.LfpStruct(pp).Sorted.Full.Infra.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1));
                % very awful but required by parfor transparency constraints
                currSGIchs=nan(3,16); % set to 16, higher than every max_{j=I,G,S}{length(currjchs)}, allows to cat along 1st dim.
                currSGIchs(1,1:length(currSchs))=currSchs;
                currSGIchs(2,1:length(currGchs))=currGchs;
                currSGIchs(3,1:length(currIchs))=currIchs;
                currSGIalgns=nan(3,16);
                currSGIalgns(1,1:length(currSchs))=subjectData.LfpStruct(pp).Sorted.Labels.SupraChsLamAlignsBi;
                currSGIalgns(2,1:length(currGchs))=subjectData.LfpStruct(pp).Sorted.Labels.GranrChsLamAlignsBi;
                currSGIalgns(3,1:length(currIchs))=subjectData.LfpStruct(pp).Sorted.Labels.InfraChsLamAlignsBi;
                % De-trend along trial-dimension may help for stationarity aspects
                currXDetr=mvdetrend(currXp);
                % De-mean along time-dimension
                currXDemn=currXDetr-repmat(mean(currXDetr,2),[1 NUM_TIME_SAMPLES 1]);
                % Downsampling by movmean
                currXDemnp=permute(currXDemn,[2 1 3]);
                currXDownp=movmeandecimatrix(currXDemnp,DOWNSAMPLING_FACTOR);
                currXDwSamp=permute(currXDownp,[2 1 3]);
                % Downsampling by brute decimation
                %currXDwSamp=currXDemn(:,1:DOWNSAMPLING_FACTOR:end,:);
                [currNvars,currNobs,currNtrials]=size(currXDwSamp);
                currfpw=nan(currNvars,currNvars,NUM_TIME_DWSAMPLES+1);
                currfpwSh=nan(NShuffles,currNvars,currNvars,NUM_TIME_DWSAMPLES+1);
                currFFiltHilbPowDown=nan(currNvars,NUM_BAND_WINS,currNtrials,NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR);
                
                % BINNING AND INFORMATION EVALUATION
                paramsMI=[];
                paramsMI.method=MImethod;
                paramsMI.bias=MIbias;
                paramsMI.btsp=MIbtsp;
                currMIMat=nan(NUM_BAND_WINS,NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR,currNvars,currNvars);
                for bw=1:NUM_BAND_WINS % LOOP OVER BAND WINDOWS
                    for cch=1:currNvars
                        % Filter for each band window
                        currBwinFFilt=filtfilt(filterBCoeffs(bw,:),filterACoeffs(bw,:),squeeze(currXDwSamp(cch,:,:)));
                        % Hilbert Transform and abs^2
                        currBwinFFiltHilbPow=(abs(hilbertpad(currBwinFFilt)).^2);
                        % Downsampling by moving average [ch x NUM_BAND_WINS x trials x NUM_TIME_DWSAMAPLES]
                        currFFiltHilbPowDown(cch,bw,:,:)=movmeandecimatrix(currBwinFFiltHilbPow,DOWNSAMPLING_FACTOR)';
                    end
                    % [currNvars x NBWINS x currNtrials x NTIMEDWS]
                    currBwXDwMat=reshape(vectorisen(currFFiltHilbPowDown(:,bw,:,:),[2 3]),[1 currNvars*currNtrials*NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR]);
                    currChsLabels=reshape(repmat(1:currNvars,currNtrials*NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR,1)',[1 currNvars*currNtrials*NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR]);
                    % [1 x (currNtrials*NTIMEDWS) x currNvars]
                    [currRespMatAllChs, currnTrAllChs]= buildr(currChsLabels,currBwXDwMat);
                    currRespMatBinAllChs = binr(currRespMatAllChs, currnTrAllChs, numQuantBins, 'eqpop');
                    currRespMatBinAllChsRes=reshape(currRespMatBinAllChs(1,:,:),[currNtrials NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR currNvars]);
                    paramsMI.nt=currNtrials;%currnTrAllChs(1);
                    for tt=1:(NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR)
                        for ch1=1:currNvars
                            for ch2=ch1:currNvars
                                currMIMat(bw,tt,ch1,ch2)=information(reshape(currRespMatBinAllChsRes(:,tt,[ch1 ch2]),[1 currNtrials 2]), paramsMI, 'I');
                            end
                        end
                    end
                end
                
                currTimeCumulInformMat=squeeze(nansum(vectorisen(currMIMat,[1 2]),2));%squeeze(nansum(currMIMat,1));
                currTimeCumulInformMat2=currTimeCumulInformMat+triu(currTimeCumulInformMat',-1);
                
                currSrtMIMat=nan(currNvars,currNvars);
                currSrtMIChs=nan(currNvars,currNvars);
                currSrtRangeNormMIMat=nan(currNvars,currNvars);
                for ch=1:currNvars
                    [currSrtMIMat(ch,:),currSrtMIChs(ch,:)]=sort(currTimeCumulInformMat2(ch,:),'descend');
                end
                
                for chY=1:currNvars
                    for chX=1:currNvars
                        if chX~=chY %required by parfor, ch2=ch1:nvars is not allowed
                            chsZOut=[];
                            if ~(any(chX==currSchs) || any(chY==currSchs))
                                chsZOut=currSchs;
                            end
                            if  ~(any(chX==currGchs) || any(chY==currGchs))
                                chsZOut=[chsZOut currGchs];
                            end
                            if  ~(any(chX==currIchs) || any(chY==currIchs))
                                chsZOut=[chsZOut currIchs];
                            end
                            
                            if numInformChs==0
                                chsZ=[];
                            elseif numInformChs>0
                                chsZMIwrtX=currSrtMIChs(chX,1:min(numInformChs,end));
                                chsZMIwrtY=currSrtMIChs(chY,1:min(numInformChs,end));
                                chsZ=setdiff(unique([chsZMIwrtX,chsZMIwrtY],'stable'),[chX chY]);
                            else
                                if numInformChs==-inf
                                    chsZ=chsZOut;
                                else
                                    numInformChsAbs=abs(numInformChs);
                                    chsZMIwrtX=intersect(currSrtMIChs(chX,1:min(numInformChsAbs,end)),chsZOut,'stable');
                                    chsZMIwrtY=intersect(currSrtMIChs(chY,1:min(numInformChsAbs,end)),chsZOut,'stable');
                                    chsZ=unique([chsZMIwrtX,chsZMIwrtY],'stable');
                                end
                            end
                            numZvars=length(chsZ);
                            
                            disp([sSubjectName '-' sVisualArea '-PEN-' num2str(subjectData.LfpStruct(pp).penNum) '-' gratCondTag{cnd} '-' compMethod '-cGC_( ' num2str(chY) ' -> ' num2str(chX) ' \ ' num2str(chsZ) ')']);
                            
                            % FIT A VAR(p) MODEL FOR (X,Y,Z) - (Barnett & Seth, J. Neurosci. Methods, 2014)
                            [BxyzMat,ExyzMat]=tsdata_to_var(currXDwSamp([chX chY chsZ],:,:),morder,regmode);
                            assert(~isbad(BxyzMat),'VAR estimation failed');
                            
                            if strcmpi(compMethod,'Mvgc') && (chX < chY)
                                % Autocovariance calculation (<mvgc|A5|>)
                                [GammaMat,GammaInfo] = var_to_autocov(BxyzMat,ExyzMat,acmaxlags,acdectol);
                                var_info(GammaInfo,true); % report results (and bail out on error)
                                
                                % Pairwise Granger causality calculation: frequency domain  (<mvgc|A14|>)
                                f_multichs_fres_21 = autocov_to_smvgc(GammaMat,2,1,fresEstim); % from 1 to 2
                                f_multichs_fres_12 = autocov_to_smvgc(GammaMat,1,2,fresEstim); % from 2 to 1
                                % assert(~isbad(f_multichs_fres,false),'Mvgc spectral GC calculation failed');
                                % if fres=[] the freqAxis ends up being non-uniform
                                % currch1ch2freqAxis=linspace(0,DWSAMPLING_FREQ/2,size(multichsfpw,3));
                                % so interp1 is used to set a uniform output fres
                                
                                f_y_to_x_given_z=interp1(lambdafreqAxis,squeeze(f_multichs_fres_12),freqAxisPreFirstDimDwSamp,'spline');
                                f_x_to_y_given_z=interp1(lambdafreqAxis,squeeze(f_multichs_fres_21),freqAxisPreFirstDimDwSamp,'spline');
                                currfpw(chX,chY,:)=abs(f_y_to_x_given_z);
                                currfpw(chY,chX,:)=abs(f_x_to_y_given_z);
                                
                                if reComputeGCsShuffling %SHUFFLE TRIALS
                                    for nnsh=1:NShuffles
                                        currXDwSampShuf=nan(size(currXDwSamp([chX chY chsZ],:,:)));
                                        for chshuf=1:length([chX chY chsZ])
                                            currXDwSampShuf(chshuf,:,:)=currXDwSamp(chshuf,:,randperm(currNtrials));
                                        end
                                        [BxyzMatShuf,ExyzMatShuf]=tsdata_to_var(currXDwSampShuf,morder,regmode);
                                        assert(~isbad(BxyzMatShuf),'VAR estimation failed');
                                        [GammaMatShuf,GammaInfoShuf] = var_to_autocov(BxyzMatShuf,ExyzMatShuf,acmaxlags,acdectol);
                                        %f_multichs_fres_Sh = autocov_to_spwcgc(GammaMatShuf,fresEstim);
                                        f_multichs_fres_Sh_21 = autocov_to_smvgc(GammaMatShuf,2,1,fresEstim);
                                        f_multichs_fres_Sh_12 = autocov_to_smvgc(GammaMatShuf,1,2,fresEstim);
                                        %assert(~isbad(f_multichs_fres_Sh,false),'Mvgc spectral GC calculation failed');
                                        f_y_to_x_given_z_sh=interp1(lambdafreqAxis,squeeze(f_multichs_fres_Sh_12),freqAxisPreFirstDimDwSamp,'spline');
                                        f_x_to_y_given_z_sh=interp1(lambdafreqAxis,squeeze(f_multichs_fres_Sh_21),freqAxisPreFirstDimDwSamp,'spline');
                                        currfpwSh(nnsh,chX,chY,:)=abs(f_y_to_x_given_z_sh);
                                        currfpwSh(nnsh,chY,chX,:)=abs(f_x_to_y_given_z_sh);
                                        % max(currfpwSh(:,chX,chY,:),[],4) is to be used for prctile
                                    end
                                end
                            end
                            if strcmpi(compMethod,'Geweke')
                                % ORIGINAL GEWEKE SPECTRAL FACTORIZATION METHOD
                                % (Geweke, J. of the American Stat. Assoc., 1984)
                                % FIT A SECOND, REDUCED VAR(p) MODEL FOR (X,Z)
                                [DxzMat,ExzMat]=tsdata_to_var(currXDwSamp([chX chsZ],:,:),morder,regmode);
                                assert(~isbad(DxzMat),'VAR estimation failed');
                                
                                % P MATRICES FOR COVARIANCE NORMALIZATION
                                PxzMat=[1 zeros(1,numZvars); -ExzMat(1,2:end)'/(ExzMat(1,1)) eye(numZvars)];
                                PxyzMat=[1 0 zeros(1,numZvars);     -ExyzMat(1,2)/(ExyzMat(1,1)) 1 zeros(1,numZvars);    -ExyzMat(1,3:end)'/(ExyzMat(1,1)) zeros(numZvars,1) eye(numZvars)]*...
                                    [1 0 zeros(1,numZvars);     0 1 zeros(1,numZvars);      zeros(numZvars,1) -(ExyzMat(3:end,2) - ExyzMat(3:end,1)/ExyzMat(1,1)*ExyzMat(1,2))/(ExyzMat(2,2) - ExyzMat(2,1)/ExyzMat(1,1)*ExyzMat(1,2)) eye(numZvars)];
                                
                                % APPLYING NORMALIZATION BY P-MATRICES MULTIPLICATION IN TIME
                                PExzMat=PxzMat*ExzMat*PxzMat';
                                PExyzMat=PxyzMat*ExyzMat*PxyzMat';
                                PDxzMat=nan(numZvars+1,numZvars+1,morder);  PBxyzMat=nan(numZvars+2,numZvars+2,morder);
                                for mm=1:morder
                                    PDxzMat(:,:,mm)=PxzMat*DxzMat(:,:,mm);
                                    PBxyzMat(:,:,mm)=PxyzMat*BxyzMat(:,:,mm);
                                end
                                [PSxzMat,PGxzMat]=var_to_cpsd(PDxzMat,PExzMat,fresEstim);
                                [PSxyzMat,PHxyzMat]=var_to_cpsd(PBxyzMat,PExyzMat,fresEstim);
                                
                                % APPLYING NORMALIZATION BY P-MATRICES MULTIPLICATION IN FREQUENCY
                                %[SxzMat,GxzMat]=var_to_cpsd(DxzMat,ExzMat,fresEstim);
                                %[SxyzMat,HxyzMat]=var_to_cpsd(BxyzMat,ExyzMat,fresEstim);
                                
                                f_y_to_x_given_z_fres=nan(1,fresEstim+1);
                                for ll=1:fresEstim+1
                                    % NORMALIZED TRANSFER FUNCTIONS P*G(X,Z) and P2*P1*H(X,Y,Z)
                                    %PGxzMat=GxzMat(:,:,lambda)/(PxzMat); %G*inv(P)
                                    %PHxyzMat=HxyzMat(:,:,lambda)/(PxyzMat); %H*inv(P)
                                    % EXPANDED G MATRIX FOR SPECTRAL FACTORIZATION
                                    GxyzMat=[PGxzMat(1,1,ll) 0 PGxzMat(1,2:end,ll);  0 1 zeros(1,numZvars); PGxzMat(2:end,1,ll) zeros(numZvars,1) PGxzMat(2:end,2:end,ll)];
                                    QxyzMat=GxyzMat\PHxyzMat(:,:,ll); %inv(G)*H
                                    %SIGxyz=abs(QxyzMat(1,1)*ExyzMat(1,1)*conj(QxyzMat(1,1)) + QxyzMat(1,2)*ExyzMat(2,2)*conj(QxyzMat(1,2)) + QxyzMat(1,3:end)*ExyzMat(3:end,3:end)*conj(QxyzMat(1,3:end))');
                                    SIGxyz=abs(QxyzMat(1,1)*PExyzMat(1,1)*conj(QxyzMat(1,1)) + QxyzMat(1,2)*PExyzMat(2,2)*conj(QxyzMat(1,2)) + QxyzMat(1,3:end)*PExyzMat(3:end,3:end)*conj(QxyzMat(1,3:end))');
                                    %SIGx=abs(QxyzMat(1,1)*ExyzMat(1,1)*conj(QxyzMat(1,1)));
                                    SIGx=abs(QxyzMat(1,1)*PExyzMat(1,1)*conj(QxyzMat(1,1)));
                                    
                                    if QxyzMat(1,1)*PExyzMat(1,1)*conj(QxyzMat(1,1)) <0
                                        warning('Qxx product is negative');
                                    end
                                    if QxyzMat(1,2)*PExyzMat(2,2)*conj(QxyzMat(1,2)) <0
                                        warning('Qxy product is negative');
                                    end
                                    if QxyzMat(1,3:end)*PExyzMat(3:end,3:end)*conj(QxyzMat(1,3:end)') <0
                                        warning('Qxz product is negative');
                                    end
                                    f_y_to_x_given_z_fres(ll)= log( SIGxyz ./ SIGx );
                                    %f_y_to_x_given_z_fres(lambda)= log( abs(PExzMat(1,1)) ./ SIGx );
                                end
                                f_y_to_x_given_z=interp1(lambdafreqAxis,f_y_to_x_given_z_fres,freqAxisPreFirstDimDwSamp,'spline');
                                currfpw(chX,chY,:)=f_y_to_x_given_z;
                            end
                            if strcmpi(compMethod,'MatPart')
                                % MATRIX PARTITIONING METHOD (Chen et al., J. Neurosci Methods 2006)
                                % P MATRIX FOR COVARIANCE NORMALIZATION
                                PxyzMat=[1 0 zeros(1,numZvars);     -ExyzMat(1,2)/(ExyzMat(1,1)) 1 zeros(1,numZvars);    -ExyzMat(1,3:end)'/(ExyzMat(1,1)) zeros(numZvars,1) eye(numZvars)]*...
                                    [1 0 zeros(1,numZvars);     0 1 zeros(1,numZvars);      zeros(numZvars,1) -(ExyzMat(3:end,2) - ExyzMat(3:end,1)/ExyzMat(1,1)*ExyzMat(1,2))/(ExyzMat(2,2) - ExyzMat(2,1)/ExyzMat(1,1)*ExyzMat(1,2)) eye(numZvars)];
                                
                                % APPLYING NORMALIZATION BY P-MATRICES MULTIPLICATION IN TIME
                                PExyzMat=PxyzMat*ExyzMat*PxyzMat';
                                PBxyzMat=nan(numZvars+2,numZvars+2,morder);
                                for mm=1:morder
                                    PBxyzMat(:,:,mm)=PxyzMat*BxyzMat(:,:,mm);
                                end
                                [PSxyzMat,PHxyzMat]=var_to_cpsd(PBxyzMat,PExyzMat,fresEstim);
                                
                                f_y_to_x_given_z_fres=nan(1,fresEstim+1);
                                for ll=1:fresEstim+1
                                    %PHxyzMat=HxyzMat(:,:,lambda)/(PxyzMat); %H*inv(P)
                                    % MIXED TRANSFER FUNCTIONS HBARxy_zy=[HBARxy; HBARzy]
                                    HBAR_xy_zy=([PHxyzMat(1,1,ll) PHxyzMat(1,3:end,ll); PHxyzMat(3:end,1,ll) PHxyzMat(3:end,3:end,ll)]) \ [PHxyzMat(1,2,ll); PHxyzMat(3:end,2,ll)];
                                    % REDUCED COVARIANCE MATRIX SIGBAR_xz=[SIGBARxx SIGBARxz; SIGBARzx SIGBARzz]
                                    EBAR_xz=[PExyzMat(1,1) PExyzMat(1,3:end); PExyzMat(3:end,1) PExyzMat(3:end,3:end)]+HBAR_xy_zy*[PExyzMat(1,2)' PExyzMat(3:end,2)']+[PExyzMat(1,2)' PExyzMat(3:end,2)']'*conj(HBAR_xy_zy')+PExyzMat(2,2)*HBAR_xy_zy*conj(HBAR_xy_zy');
                                    % REDUCED P MATRIX FOR COVARIANCE NORMALIZATION
                                    PBARMat_xz=[1 zeros(1,numZvars); -EBAR_xz(2:end,1)/EBAR_xz(1,1) eye(numZvars)];
                                    GMat_xz=[PHxyzMat(1,1,ll) PHxyzMat(1,3:end,ll); PHxyzMat(3:end,1,ll) PHxyzMat(3:end,3:end,ll)]/(PBARMat_xz); % G/P = G*inv(P)
                                    % EXPANDED G MATRIX FOR SPECTRAL FACTORIZATION
                                    GMat_xyz=[GMat_xz(1,1) 0 GMat_xz(1,2:end); 0 1 zeros(1,numZvars); GMat_xz(2:end,1) zeros(numZvars,1) GMat_xz(2:end,2:end)];
                                    Q_xyz=(GMat_xyz)\PHxyzMat(:,:,ll); % G\H = inv(G)*H
                                    %SIGxyz=abs( Q_xyz(1,1) * ExyzMat(1,1) * conj(Q_xyz(1,1)) +  Q_xyz(1,2) * ExyzMat(2,2) * conj(Q_xyz(1,2)) + Q_xyz(1,3:end) * ExyzMat(3:end,3:end) * conj(Q_xyz(1,3:end))');
                                    SIGxyz=abs( Q_xyz(1,1) * PExyzMat(1,1) * conj(Q_xyz(1,1)) +  Q_xyz(1,2) * PExyzMat(2,2) * conj(Q_xyz(1,2)) + Q_xyz(1,3:end) * PExyzMat(3:end,3:end) * conj(Q_xyz(1,3:end))');
                                    %SIGx=abs( Q_xyz(1,1) * ExyzMat(1,1) * conj(Q_xyz(1,1)) );
                                    SIGx=abs( Q_xyz(1,1) * PExyzMat(1,1) * conj(Q_xyz(1,1)) );
                                    f_y_to_x_given_z_fres(ll)=log( SIGxyz ./ SIGx );
                                    %f_y_to_x_given_z_fres(lambda)=log( abs(EBAR_xz(1,1)) ./ SIGx );
                                end
                                f_y_to_x_given_z=interp1(lambdafreqAxis,f_y_to_x_given_z_fres,freqAxisPreFirstDimDwSamp,'spline');
                                currfpw(chX,chY,:)=f_y_to_x_given_z;
                            end
                        end
                    end
                end
                currGCsStruct(pp).(gratCondTag{cnd}).pwgcf=currfpw;
                currGCsStruct(pp).(gratCondTag{cnd}).pwgcfSh=currfpwSh;
            end
            currGCsStruct(pp).penNum=subjectData.LfpStruct(pp).penNum;
            currGCsStruct(pp).currSGIchs=currSGIchs;
            currGCsStruct(pp).currSGIalgns=currSGIalgns;
        end
    end
    
    % SORTING BY VISUAL AREA / LAMINAR LAYER / GRATC CONDITION
    for pp=1:length(currGCsStruct)
        if ~isempty(currGCsStruct(pp).currSGIchs)
            currSGIchs=currGCsStruct(pp).currSGIchs;
            currSGIalgns=currGCsStruct(pp).currSGIalgns;
            for cnd=1:2
                for ll1=1:3
                    for ll2=1:3
                        % eval is awful but it saves from extra indexing
                        currll1chs=currSGIchs(ll1,~isnan(currSGIchs(ll1,:)));
                        currll2chs=currSGIchs(ll2,~isnan(currSGIchs(ll2,:)));
                        % NOTE: reverse to get fi->j|[ij] instead of fj->i|[ij]
                        currllsfpw=currGCsStruct(pp).(gratCondTag{cnd}).pwgcf(currll2chs,currll1chs,:);
                        currllsfpwSh=currGCsStruct(pp).(gratCondTag{cnd}).pwgcfSh(:,currll2chs,currll1chs,:);
                        if ~isempty(currllsfpw)
                            % pool chsxchs within ll1xll2
                            currllsfpwRes=reshape(currllsfpw,size(currllsfpw,1)*size(currllsfpw,2),size(currllsfpw,3));
                            currllsfpwResSh=reshape(currllsfpwSh,NShuffles,size(currllsfpwSh,2)*size(currllsfpwSh,3),size(currllsfpwSh,4));
                            % get rid of all-NaNs rowsacdpt[]
                            currllsfpwRes((sum(isnan(currllsfpwRes),2)==size(currllsfpwRes,2)),:)=[];
                            currllsfpwResSh(:,sum(isnan(currllsfpwResSh(1,:,:)),3)==size(currllsfpwResSh,3),:)=[];
                            sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcf=currllsfpwRes;
                            sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcfSh=currllsfpwResSh;
                            %eval(clearoutmem('currlls*'));
                        else
                            sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcf=[];
                            sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcfSh=[];
                        end
                        currll1algns=currSGIalgns(ll1,~isnan(currSGIalgns(ll1,:)));
                        currll2algns=currSGIalgns(ll2,~isnan(currSGIalgns(ll2,:)));
                        if length(currll1chs)~=length(currll1algns) ||  length(currll2chs)~=length(currll2algns); error('Layer labels do not match alignments.'); end
                        
                        for algn1=1:length(currll1algns)
                            for algn2=1:length(currll2algns)
                                currllalgnsfpw=currGCsStruct(pp).(gratCondTag{cnd}).pwgcf(currll2chs(algn2),currll1chs(algn1),:);
                                currllalgnsfpwRes=reshape(currllalgnsfpw,1,size(currllalgnsfpw,3));
                                currllalgnsfpwSh=currGCsStruct(pp).(gratCondTag{cnd}).pwgcfSh(:,currll2chs(algn2),currll1chs(algn1),:);
                                currllalgnsfpwResSh=reshape(currllalgnsfpwSh,NShuffles,1,size(currllalgnsfpwSh,4));
                                
                                if currll1algns(algn1)<0
                                    strcurrll1algn1=['m' num2str(abs(currll1algns(algn1)))];
                                else
                                    strcurrll1algn1=['p' num2str(currll1algns(algn1))];
                                end
                                if currll2algns(algn2)<0
                                    strcurrll2algn2=['m' num2str(abs(currll2algns(algn2)))];
                                else
                                    strcurrll2algn2=['p' num2str(currll2algns(algn2))];
                                end
                                if sum(isnan(currllalgnsfpw),2)~=size(currllalgnsfpw,2)
                                    sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcf'])=currllalgnsfpwRes;
                                    sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcfSh'])=currllalgnsfpwResSh;
                                else
                                    sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcf'])=[];
                                    sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcfSh'])=[];
                                end
                            end
                        end
                    end
                    if cnd==2
                        sGCsStruct(pp).penNum=currGCsStruct(pp).penNum;
                        sGCsStruct(pp).Labels.(['num' lamLayerTag{ll1}(1) 'chs'])=length(currll1chs);
                        sGCsStruct(pp).Labels.([lamLayerTag{ll1} 'ChsAlgns'])=currll1algns;
                    end
                end
            end
        end
    end
    
    save(['./' upper(sSubjectName(1)) sSubjectName(end) sVisualArea '_cGCs_' num2str(numInformChs) '_MostInfChs_' compMethod '_PreFirstDim.mat'],'sGCsStruct','-v7.3');
    
else
    load(['./' upper(sSubjectName(1)) sSubjectName(end) sVisualArea '_cGCs_' num2str(numInformChs) '_MostInfChs_' compMethod '_PreFirstDim.mat']);
end

%% POOL ACROSS PENs GCs
PoolGCs=[];
for cnd=1:2
    PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcf=[];
    PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcfSh=[];
    for ll1=1:3
        for ll2=1:3
            currll1nchs=[];
            PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}])=[];
            for pp=1:length(sGCsStruct)
                if ~isempty(sGCsStruct(pp).RF)
                    currSubFieldNames=fieldnames(sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]));
                    currSubFieldNames=currSubFieldNames(~contains(currSubFieldNames,'Sh'));
                    for sbfs=1:length(currSubFieldNames)
                        if ~isfield(PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]),currSubFieldNames{sbfs})
                            PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs})=[];
                            PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh'])=[];
                        end
                        PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs})=...
                            [PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs});...
                            sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs})];
                        PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh'])=...
                            cat(2,PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh']),...
                            sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh']));
                    end
                    currll1nchs=cat(1,currll1nchs,sGCsStruct(pp).Labels.(['num' lamLayerTag{ll1}(1) 'chs']));
                end
            end
            if cnd==2
                PoolGCs.Labels.(['num' lamLayerTag{ll1}(1) 'chs'])=currll1nchs;
            end
            PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcf=[PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcf;...
                PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcf];
            PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcfSh=[PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcfSh;...
                PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcfSh];
        end
    end
end

%% PLOT PEN AVG GCs Layers x Layers (3x3 subplot)
if plotPenAvgGCs
    ylimplot=[0 0];
    for lltx=1:3
        for llrx=1:3
            if ~isempty(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf)
                ylimplot(1)=min([nanmean(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)-...
                    nansem(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)...
                    nanmean(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)-...
                    nansem(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1) ylimplot(1)]);
                ylimplot(2)=max([nanmean(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)+...
                    nansem(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)...
                    nanmean(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)+...
                    nansem(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1) ylimplot(2)]);
            end
        end
    end
    ylimplot=1*ylimplot; ylimplot(1)=-1e-2;
    hfig=figure('units','normalized','outerposition',[0 0 1 1]);
    for lltx=1:3
        for llrx=1:3
            if ~isempty(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf)
                for ff=1:length(freqAxisPreFirstDimDwSamp)
                    PoolGCs.pValues.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pvalues(ff)=...
                        signrank(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,ff),...
                        PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,ff));
                end
                pBarVec=PoolGCs.pValues.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pvalues(freqIndices2Plot);
                [~,~,~,pBarVecFDR]=fdr_bh(pBarVec);
                
                currLLsPwgcfRFShFreqWise=nan(1,NUM_TIME_DWSAMPLES+1);
                currLLsPwgcfOUTrpShFreqWise=nan(1,NUM_TIME_DWSAMPLES+1);
                
                if reportSignifThreshold
                    for ff=1:(NUM_TIME_DWSAMPLES+1)
                        currLLsFreqPwgcfRFSh=squeeze(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcfSh(:,:,ff));
                        currLLsPwgcfRFShFreqWise(ff)=prctile(nanmean(currLLsFreqPwgcfRFSh,2),95);
                        currLLsFreqPwgcfOUTrpSh=squeeze(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcfSh(:,:,ff));
                        currLLsPwgcfOUTrpShFreqWise(ff)=prctile(nanmean(currLLsFreqPwgcfOUTrpSh,2),95);
                    end
                end
                
                subplot(3,3,(lltx-1)*3+llrx)
                plotcmapdots([1 0 0; 0 0 1;.2 .6 .2]);
                plot([-inf -inf],[-inf -inf],'r:'); plot([-inf -inf],[-inf -inf],'b:');
                plotmsem(freqAxisPreFirstDimDwSamp(freqIndices2Plot),(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot)),'r');%[.73 0 0]);
                plotmsem(freqAxisPreFirstDimDwSamp(freqIndices2Plot),(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot)),'b');%[0 .73 .73]);
                plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),currLLsPwgcfRFShFreqWise(freqIndices2Plot),'r:');
                plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),currLLsPwgcfOUTrpShFreqWise(freqIndices2Plot),'b:');
                plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),double(pBarVecFDR<=0.05)-1+ylimplot(1),'color',[.2 .6 .2],'linewidth',3)
                xlim([0 100]);  ylim(ylimplot);
                title([lamLayerTag{lltx} ' to ' lamLayerTag{llrx}])
                if llrx<=lltx
                    text(100, .9*(ylimplot(2)-ylimplot(1))+ylimplot(1),['tot n=' num2str(size(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf,1))],'color','k');
                end
                if llrx==3 && lltx==1
                    legend({'RF','OUT','p \leq 0.05','95% sh. perc. RF','95% sh. perc. OUT'},'Location','northeast');
                end
                if llrx==1
                    text(100, .8*(ylimplot(2)-ylimplot(1))+ylimplot(1),['avg n=' sprintf('%.2f',mean(PoolGCs.Labels.(['num' lamLayerTag{lltx}(1) 'chs'])))],'color','k');
                end
            end
        end
    end
    if numInformChs==0
        supertitle([sSubjectName ' ' sVisualArea ' Unconditional ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
    elseif numInformChs==Inf
        supertitle([sSubjectName ' ' sVisualArea ' Conditional (to all Complementary Channels) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
    elseif numInformChs>0
        supertitle([sSubjectName ' ' sVisualArea ' Conditional (to n=' num2str(numInformChs) ' Most Inform. Complem. Chs) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
    elseif numInformChs==-Inf
        supertitle([sSubjectName ' ' sVisualArea ' Conditional (to Complementary Chs outside current Layers) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
    else
        supertitle([sSubjectName ' ' sVisualArea ' Conditional (to n=' num2str(abs(numInformChs)) ' Most Inf. Compl. Chs outside current Layers) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
    end
end


%% PLOT PEN AVG GCs Layers x Layers (3x3 subplot)
if plotPenAvgGCs
    ylimplot=[-1e-2 0.12];
    hfig=figure('units','normalized','outerposition',[0 0 1 1]);
    
    for ff=1:length(freqAxisPreFirstDimDwSamp)
        PoolGCs.pValues.AllLayers.pvalues(ff)=signrank(PoolGCs.RF.AllLayers.pwgcf(:,ff),PoolGCs.OUTrp.AllLayers.pwgcf(:,ff));
    end
    pBarVec=PoolGCs.pValues.AllLayers.pvalues(freqIndices2Plot);
    [~,~,~,pBarVecFDR]=fdr_bh(pBarVec);
    
    currLLsPwgcfRFShFreqWise=nan(1,NUM_TIME_DWSAMPLES+1);
    currLLsPwgcfOUTrpShFreqWise=nan(1,NUM_TIME_DWSAMPLES+1);
    
    if reportSignifThreshold
        for ff=1:(NUM_TIME_DWSAMPLES+1)
            currLLsFreqPwgcfRFSh=squeeze(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcfSh(:,:,ff));
            currLLsPwgcfRFShFreqWise(ff)=prctile(nanmean(currLLsFreqPwgcfRFSh,2),95);
            currLLsFreqPwgcfOUTrpSh=squeeze(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcfSh(:,:,ff));
            currLLsPwgcfOUTrpShFreqWise(ff)=prctile(nanmean(currLLsFreqPwgcfOUTrpSh,2),95);
        end
    end
    plotcmapdots([1 0 0; 0 0 1;.2 .6 .2]);
    plot([-inf -inf],[-inf -inf],'r:'); plot([-inf -inf],[-inf -inf],'b:');
    plotmsem(freqAxisPreFirstDimDwSamp(freqIndices2Plot),(PoolGCs.RF.AllLayers.pwgcf(:,freqIndices2Plot)),'r');%[.73 0 0]);
    plotmsem(freqAxisPreFirstDimDwSamp(freqIndices2Plot),(PoolGCs.OUTrp.AllLayers.pwgcf(:,freqIndices2Plot)),'b');%[0 .73 .73]);
    plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),currLLsPwgcfRFShFreqWise(freqIndices2Plot),'r:');
    plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),currLLsPwgcfOUTrpShFreqWise(freqIndices2Plot),'b:');
    plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),double(pBarVecFDR<=0.05)-1+ylimplot(1),'color',[.2 .6 .2],'linewidth',3)
    xlim([0 100]);  ylim(ylimplot);
    text(90, .9*(ylimplot(2)-ylimplot(1))+ylimplot(1),['tot n=' num2str(size(PoolGCs.RF.AllLayers.pwgcf,1))],'color','k');
    if llrx==3 && lltx==1
        legend({'RF','OUT','p \leq 0.05','95% sh. perc. RF','95% sh. perc. OUT'},'Location','northeast');
    end
    if numInformChs==0
        supertitle([sSubjectName ' ' sVisualArea ' Unconditional ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
    elseif numInformChs==Inf
        supertitle([sSubjectName ' ' sVisualArea ' Conditional (to all Complementary Channels) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
    elseif numInformChs>0
        supertitle([sSubjectName ' ' sVisualArea ' Conditional (to n=' num2str(numInformChs) ' Most Inform. Complem. Chs) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
    elseif numInformChs==-Inf
        supertitle([sSubjectName ' ' sVisualArea ' Conditional (to Complementary Chs outside current Layers) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
    else
        supertitle([sSubjectName ' ' sVisualArea ' Conditional (to n=' num2str(abs(numInformChs)) ' Most Inf. Compl. Chs outside current Layers) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
    end
end