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

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

clear all; close all; clc;

sSubjectName='monkey 1';
sVisualArea='V1V4';

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);
    subjectDataV1=getSubjectData(sSubjectName,'V1','ALL','LFPs','ALL');
    subjectDataV4=getSubjectData(sSubjectName,'V4','ALL','LFPs','ALL');
    
    % REMOVE DATA IN TIME WINDOWS OTHER THAN PREFIRSTDIM TO SAVE MEMORY
    rmpp=[];
    for pp=1:length(subjectDataV1.LfpStruct)
        if ~isempty(subjectDataV1.LfpStruct(pp).Sorted)
            rmfn=fieldnames(subjectDataV1.LfpStruct(pp).Sorted.Full.Supra.RF);
            rmfn(cellfun(@(fnm) strcmpi(fnm,'PreFirstDimDataBiZSc'),rmfn))=[];
            for ll=1:3
                subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF=rmfield(subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF,rmfn);
                subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1=rmfield(subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1,rmfn);
                subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2=rmfield(subjectDataV1.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2,rmfn);
            end
        else
            rmpp=cat(1,rmpp,pp);
        end
    end
    subjectDataV1.LfpStruct(rmpp)=[]; subjectDataV1.penInfos(rmpp)=[]; subjectDataV1.penIDs(rmpp)=[];
    rmpp=[];
    for pp=1:length(subjectDataV4.LfpStruct)
        if ~isempty(subjectDataV4.LfpStruct(pp).Sorted)
            rmfn=fieldnames(subjectDataV4.LfpStruct(pp).Sorted.Full.Supra.RF);
            rmfn(cellfun(@(fnm) strcmpi(fnm,'PreFirstDimDataBiZSc'),rmfn))=[];
            for ll=1:3
                subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF=rmfield(subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF,rmfn);
                subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1=rmfield(subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1,rmfn);
                subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2=rmfield(subjectDataV4.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2,rmfn);
            end
        else
            rmpp=cat(1,rmpp,pp);
        end
    end
    subjectDataV4.LfpStruct(rmpp)=[]; subjectDataV4.penInfos(rmpp)=[]; subjectDataV4.penIDs(rmpp)=[];
    
    %% FIND PENs W/SAME TRIALS IN V1 AND V4
    penIDsV1=arrayfun(@(jj) subjectDataV1.LfpStruct(jj).penNum,1:length(subjectDataV1.LfpStruct));
    penIDsV4=arrayfun(@(jj) subjectDataV4.LfpStruct(jj).penNum,1:length(subjectDataV4.LfpStruct));
    penIDsV1V4=intersect(penIDsV1,penIDsV4);% penIDsV1(arrayfun(@(jj) any(penIDsV4==penIDsV1(jj)), 1:length(penIDsV1)));
    penIndicesV1=arrayfun(@(jj) find(penIDsV1==penIDsV1V4(jj)), 1:length(penIDsV1V4));
    penIndicesV4=arrayfun(@(jj) find(penIDsV4==penIDsV1V4(jj)), 1:length(penIDsV1V4));
    
    sLfpStructV1V4=[];
    parfor pp=1:length(penIDsV1V4) % LOOP OVER PENS
        if (subjectDataV1.LfpStruct(penIndicesV1(pp)).penNum==subjectDataV4.LfpStruct(penIndicesV4(pp)).penNum)
            sLfpStructV1V4(pp).penNum=subjectDataV1.LfpStruct(penIndicesV1(pp)).penNum;
        else
            error('PEN Numbers in V1 and V4 do not match.');
        end
        for cnd=1:3 % LOOP OVER ATTENTIONAL CONDITIONS
            % FIND TRIALS WITH SAME EVENT IDs FOR V1 AND V4
            currPenCndTrialIDsV1=subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Labels.selectedNlxEvents(remod(ceil(subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Labels.GratCondition/6),3)==cnd);
            currPenCndTrialIDsV4=subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Labels.selectedNlxEvents(remod(ceil(subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Labels.GratCondition/6),3)==cnd);
            currPenCndTrialIDsV1V4=intersect(currPenCndTrialIDsV1,currPenCndTrialIDsV4);
            currPenCndTrialIndicesV1=arrayfun(@(jj) find(currPenCndTrialIDsV1==currPenCndTrialIDsV1V4(jj)), 1:length(currPenCndTrialIDsV1V4));
            currPenCndTrialIndicesV4=arrayfun(@(jj) find(currPenCndTrialIDsV4==currPenCndTrialIDsV1V4(jj)), 1:length(currPenCndTrialIDsV1V4));
            if ~all(currPenCndTrialIDsV1(currPenCndTrialIndicesV1) == currPenCndTrialIDsV4(currPenCndTrialIndicesV4))
                error('Some of the trials IDs for V1 do not match trial IDs for V4.');
            end
            
            % FIND TRIALS WITH SAME FIRST DIM DELAY FOR V1 AND V4
            currPenCndFirstDimDelayV1=subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Delays.FirstDim(remod(ceil(subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Labels.GratCondition/6),3)==cnd);
            currPenCndFirstDimDelayV4=subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Delays.FirstDim(remod(ceil(subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Labels.GratCondition/6),3)==cnd);
            if ~all(currPenCndFirstDimDelayV1(currPenCndTrialIndicesV1)==currPenCndFirstDimDelayV4(currPenCndTrialIndicesV4))
                error('FirstDim Delays not matching in V1 and V4 trials');
            end
            
            for lls=1:3 % LOOP OVER LAMINAR LAYERS PICK SIMULTANEOUSLY VALID TRIALS
                cnd134=cnd; if cnd>1; cnd134=cnd+1; end
                sLfpStructV1V4(pp).V1.Sorted.Full.(lamLayerTag{lls}).(gratCondTag{cnd134}).PreFirstDimDataBiZSc=...
                    subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Full.(lamLayerTag{lls}).(gratCondTag{cnd134}).PreFirstDimDataBiZSc(:,currPenCndTrialIndicesV1,:);
                sLfpStructV1V4(pp).V4.Sorted.Full.(lamLayerTag{lls}).(gratCondTag{cnd134}).PreFirstDimDataBiZSc=...
                    subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Full.(lamLayerTag{lls}).(gratCondTag{cnd134}).PreFirstDimDataBiZSc(:,currPenCndTrialIndicesV4,:);
            end
        end
        sLfpStructV1V4(pp).V1.Sorted.Labels.SupraChsLamAlignsBi=subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Labels.SupraChsLamAlignsBi;
        sLfpStructV1V4(pp).V1.Sorted.Labels.GranrChsLamAlignsBi=subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Labels.GranrChsLamAlignsBi;
        sLfpStructV1V4(pp).V1.Sorted.Labels.InfraChsLamAlignsBi=subjectDataV1.LfpStruct(penIndicesV1(pp)).Sorted.Labels.InfraChsLamAlignsBi;
        sLfpStructV1V4(pp).V4.Sorted.Labels.SupraChsLamAlignsBi=subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Labels.SupraChsLamAlignsBi;
        sLfpStructV1V4(pp).V4.Sorted.Labels.GranrChsLamAlignsBi=subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Labels.GranrChsLamAlignsBi;
        sLfpStructV1V4(pp).V4.Sorted.Labels.InfraChsLamAlignsBi=subjectDataV4.LfpStruct(penIndicesV4(pp)).Sorted.Labels.InfraChsLamAlignsBi;
    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(penIDsV1V4)
        currNtrRF=size(sLfpStructV1V4(pp).V1.Sorted.Full.Supra.RF.PreFirstDimDataBiZSc,2);
        currNtrOUT1=size(sLfpStructV1V4(pp).V1.Sorted.Full.Supra.OUT1.PreFirstDimDataBiZSc,2);
        currNtrOUT2=size(sLfpStructV1V4(pp).V1.Sorted.Full.Supra.OUT2.PreFirstDimDataBiZSc,2);
        currV1V4OUTrp=randperm(currNtrOUT1+currNtrOUT2,currNtrRF);
        for lls=1:3
            currOUTcatV1=cat(2,sLfpStructV1V4(pp).V1.Sorted.Full.(lamLayerTag{lls}).OUT1.PreFirstDimDataBiZSc,...
                sLfpStructV1V4(pp).V1.Sorted.Full.(lamLayerTag{lls}).OUT2.PreFirstDimDataBiZSc);
            currOUTcatV4=cat(2,sLfpStructV1V4(pp).V4.Sorted.Full.(lamLayerTag{lls}).OUT1.PreFirstDimDataBiZSc,...
                sLfpStructV1V4(pp).V4.Sorted.Full.(lamLayerTag{lls}).OUT2.PreFirstDimDataBiZSc);
            sLfpStructV1V4(pp).V1.Sorted.Full.(lamLayerTag{lls}).OUTrp.PreFirstDimDataBiZSc=currOUTcatV1(:,currV1V4OUTrp,:);
            sLfpStructV1V4(pp).V4.Sorted.Full.(lamLayerTag{lls}).OUTrp.PreFirstDimDataBiZSc=currOUTcatV4(:,currV1V4OUTrp,:);
        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
    
    for pp=1:length(sLfpStructV1V4)
        for cnd=1:2
            disp(['Processing cGC for ' sSubjectName ' ' sVisualArea ' PEN ' num2str(sLfpStructV1V4(pp).penNum) ', ' gratCondTag{cnd}]);
            
            currX0=[sLfpStructV1V4(pp).V1.Sorted.Full.Supra.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end);...
                sLfpStructV1V4(pp).V1.Sorted.Full.Granr.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end);...
                sLfpStructV1V4(pp).V1.Sorted.Full.Infra.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end);...
                sLfpStructV1V4(pp).V4.Sorted.Full.Supra.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end);...
                sLfpStructV1V4(pp).V4.Sorted.Full.Granr.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end);...
                sLfpStructV1V4(pp).V4.Sorted.Full.Infra.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end)];
            %swap time and trial dimensions
            currXp=permute(currX0,[1 3 2]);
            %eval(clearoutmem('currX0'));
            
            currV1Schs=1:size(sLfpStructV1V4(pp).V1.Sorted.Full.Supra.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1);
            currV1Gchs=max([currV1Schs 0])+(1:size(sLfpStructV1V4(pp).V1.Sorted.Full.Granr.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1));
            currV1Ichs=max([currV1Gchs currV1Schs 0])+(1:size(sLfpStructV1V4(pp).V1.Sorted.Full.Infra.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1));
            currV4Schs=max([currV1Ichs currV1Gchs currV1Schs 0])+(1:size(sLfpStructV1V4(pp).V4.Sorted.Full.Supra.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1));
            currV4Gchs=max([currV4Schs currV1Ichs currV1Gchs currV1Schs 0])+(1:size(sLfpStructV1V4(pp).V4.Sorted.Full.Granr.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1));
            currV4Ichs=max([currV4Gchs currV4Schs currV1Ichs currV1Gchs currV1Schs 0])+(1:size(sLfpStructV1V4(pp).V4.Sorted.Full.Infra.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1));
            % required if you run parfor pp=1:length(sLfpStructV1V4) for transparency constraints
            currV1SGIchs=nan(3,16); % set to 16, higher than every max_{j=I,G,S}{length(currjchs)}, allows to cat along 1st dim.
            currV1SGIchs(1,1:length(currV1Schs))=currV1Schs;
            currV1SGIchs(2,1:length(currV1Gchs))=currV1Gchs;
            currV1SGIchs(3,1:length(currV1Ichs))=currV1Ichs;
            currV4SGIchs=nan(3,16);
            currV4SGIchs(1,1:length(currV4Schs))=currV4Schs;
            currV4SGIchs(2,1:length(currV4Gchs))=currV4Gchs;
            currV4SGIchs(3,1:length(currV4Ichs))=currV4Ichs;
            currV1SGIalgns=nan(3,16);
            currV1SGIalgns(1,1:length(currV1Schs))=sLfpStructV1V4(pp).V1.Sorted.Labels.SupraChsLamAlignsBi;
            currV1SGIalgns(2,1:length(currV1Gchs))=sLfpStructV1V4(pp).V1.Sorted.Labels.GranrChsLamAlignsBi;
            currV1SGIalgns(3,1:length(currV1Ichs))=sLfpStructV1V4(pp).V1.Sorted.Labels.InfraChsLamAlignsBi;
            currV4SGIalgns=nan(3,16);
            currV4SGIalgns(1,1:length(currV4Schs))=sLfpStructV1V4(pp).V4.Sorted.Labels.SupraChsLamAlignsBi;
            currV4SGIalgns(2,1:length(currV4Gchs))=sLfpStructV1V4(pp).V4.Sorted.Labels.GranrChsLamAlignsBi;
            currV4SGIalgns(3,1:length(currV4Ichs))=sLfpStructV1V4(pp).V4.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,:,:)));
                    %currBwinFFiltABCtot(:,bw,:)=currBwinFFiltABC;
                    
                    % Hilbert Transform and abs^2
                    currBwinFFiltHilbPow=(abs(hilbertpad(currBwinFFilt)).^2);
                    %currBwinFFiltHilbPowtot(:,bw,:)=currBwinFFiltHilbPow;
                    
                    % 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==currV1Schs) || any(chY==currV1Schs))
                            chsZOut=currV1Schs;
                        end
                        if  ~(any(chX==currV1Gchs) || any(chY==currV1Gchs))
                            chsZOut=[chsZOut currV1Gchs];
                        end
                        if  ~(any(chX==currV1Ichs) || any(chY==currV1Ichs))
                            chsZOut=[chsZOut currV1Ichs];
                        end
                        if ~(any(chX==currV4Schs) || any(chY==currV4Schs))
                            chsZOut=[chsZOut currV4Schs];
                        end
                        if  ~(any(chX==currV4Gchs) || any(chY==currV4Gchs))
                            chsZOut=[chsZOut currV4Gchs];
                        end
                        if  ~(any(chX==currV4Ichs) || any(chY==currV4Ichs))
                            chsZOut=[chsZOut currV4Ichs];
                        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(sLfpStructV1V4(pp).penNum) '-' gratCondTag{cnd} '-' compMethod '-cGC_( ' num2str(chY) ' -> ' num2str(chX) ' \ ' num2str(chsZ) ')']);
                        
                        % FIT A VAR(p) MODEL FOR (X,Y,Z)
                        [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 = autocov_to_spwcgc(GammaMat,fresEstim);
                            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 lambda=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,lambda) 0 PGxzMat(1,2:end,lambda);  0 1 zeros(1,numZvars); PGxzMat(2:end,1,lambda) zeros(numZvars,1) PGxzMat(2:end,2:end,lambda)];
                                QxyzMat=GxyzMat\PHxyzMat(:,:,lambda); %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(lambda)= 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 (Y. Chen et al., JNeurosci 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 lambda=1:fresEstim+1
                                %PHxyzMat=HxyzMat(:,:,lambda)/(PxyzMat); %H*inv(P)
                                % MIXED TRANSFER FUNCTIONS HBARxy_zy=[HBARxy; HBARzy]
                                HBAR_xy_zy=([PHxyzMat(1,1,lambda) PHxyzMat(1,3:end,lambda); PHxyzMat(3:end,1,lambda) PHxyzMat(3:end,3:end,lambda)]) \ [PHxyzMat(1,2,lambda); PHxyzMat(3:end,2,lambda)];
                                % 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,lambda) PHxyzMat(1,3:end,lambda); PHxyzMat(3:end,1,lambda) PHxyzMat(3:end,3:end,lambda)]/(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(:,:,lambda); % 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(lambda)=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
            currGCsStructV1V4(pp).(gratCondTag{cnd}).pwgcf=currfpw;
            currGCsStructV1V4(pp).(gratCondTag{cnd}).pwgcfSh=currfpwSh;
        end
        currGCsStructV1V4(pp).currV1SGIchs=currV1SGIchs;
        currGCsStructV1V4(pp).currV4SGIchs=currV4SGIchs;
        currGCsStructV1V4(pp).currV1SGIalgns=currV1SGIalgns;
        currGCsStructV1V4(pp).currV4SGIalgns=currV4SGIalgns;
        currGCsStructV1V4(pp).penNum=sLfpStructV1V4(pp).penNum;
    end
    
    
    % SORTING BY VISUAL AREA / LAMINAR LAYER / GRATC CONDITION
    for pp=1:length(sLfpStructV1V4)
        currV1SGIchs=currGCsStructV1V4(pp).currV1SGIchs;
        currV4SGIchs=currGCsStructV1V4(pp).currV4SGIchs;
        currV1SGIalgns=currGCsStructV1V4(pp).currV1SGIalgns;
        currV4SGIalgns=currGCsStructV1V4(pp).currV4SGIalgns;
        for cnd=1:2
            for vv1=[1 4]
                for vv2=[1 4]
                    for ll1=1:3
                        for ll2=1:3
                            % eval is awful but it saves from extra indexing
                            currll1chs=eval(['currV' num2str(vv1) 'SGIchs(ll1,~isnan(currV' num2str(vv1) 'SGIchs(ll1,:)))']);
                            currll2chs=eval(['currV' num2str(vv2) 'SGIchs(ll2,~isnan(currV' num2str(vv2) 'SGIchs(ll2,:)))']);
                            % NOTE: reverse to get fi->j|[ij] instead of fj->i|[ij]
                            currllsfpw=currGCsStructV1V4(pp).(gratCondTag{cnd}).pwgcf(currll2chs,currll1chs,:);
                            currllsfpwSh=currGCsStructV1V4(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 rows
                                currllsfpwRes((sum(isnan(currllsfpwRes),2)==size(currllsfpwRes,2)),:)=[];
                                currllsfpwResSh(:,sum(isnan(currllsfpwResSh(1,:,:)),3)==size(currllsfpwResSh,3),:)=[];
                                sGCsStructV1V4(pp).(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcf=currllsfpwRes;
                                sGCsStructV1V4(pp).(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcfSh=currllsfpwResSh;
                                %eval(clearoutmem('currlls*'));
                            else
                                sGCsStructV1V4(pp).(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcf=[];
                                sGCsStructV1V4(pp).(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcfSh=[];
                            end
                            currll1algns=eval(['currV' num2str(vv1) 'SGIalgns(ll1,~isnan(currV' num2str(vv1) 'SGIalgns(ll1,:)))']);
                            currll2algns=eval(['currV' num2str(vv2) 'SGIalgns(ll2,~isnan(currV' num2str(vv2) 'SGIalgns(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=currGCsStructV1V4(pp).(gratCondTag{cnd}).pwgcf(currll2chs(algn2),currll1chs(algn1),:);
                                    currllalgnsfpwRes=reshape(currllalgnsfpw,1,size(currllalgnsfpw,3));
                                    currllalgnsfpwSh=currGCsStructV1V4(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)
                                        sGCsStructV1V4(pp).(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcf'])=currllalgnsfpwRes;
                                        sGCsStructV1V4(pp).(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcfSh'])=currllalgnsfpwResSh;
                                    else
                                        sGCsStructV1V4(pp).(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcf'])=[];
                                        sGCsStructV1V4(pp).(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcfSh'])=[];
                                    end
                                end
                            end
                        end
                        if vv2==4 && cnd==2
                            sGCsStructV1V4(pp).penNum=currGCsStructV1V4(pp).penNum;
                            sGCsStructV1V4(pp).Labels.(['V' num2str(vv1)]).(['num' lamLayerTag{ll1}(1) 'chs'])=length(currll1chs);
                            sGCsStructV1V4(pp).Labels.(['V' num2str(vv1)]).([lamLayerTag{ll1} 'ChsAlgns'])=currll1algns;
                        end
                    end
                end
            end
        end
    end
    
    save(['./' upper(sSubjectName(1)) sSubjectName(end) sVisualArea '_cGCs_' num2str(numInformChs) '_MostInfChs_' compMethod '_PreFirstDim.mat'],'sGCsStructV1V4','-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
    for vv1=[1 4]
        for vv2=[1 4]
            PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcf=[];
            PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcfSh=[];
            for ll1=1:3
                for ll2=1:3
                    currll1nchs=[];
                    PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}])=[];
                    for pp=1:length(sGCsStructV1V4)
                        currSubFieldNames=fieldnames(sGCsStructV1V4(pp).(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]));
                        currSubFieldNames=currSubFieldNames(~contains(currSubFieldNames,'Sh'));
                        for sbfs=1:length(currSubFieldNames)
                            if ~isfield(PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]),currSubFieldNames{sbfs})
                                PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs})=[];
                                PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh'])=[];
                            end
                            PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs})=...
                                [PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs});...
                                sGCsStructV1V4(pp).(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs})];
                            if reportSignifThreshold
                                PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh'])=...
                                    cat(2,PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh']),...
                                    sGCsStructV1V4(pp).(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh']));
                            end
                        end
                        currll1nchs=cat(1,currll1nchs,sGCsStructV1V4(pp).Labels.(['V' num2str(vv1)]).(['num' lamLayerTag{ll1}(1) 'chs']));
                    end
                    if cnd==2
                        PoolGCs.Labels.(['V' num2str(vv1)]).(['num' lamLayerTag{ll1}(1) 'chs'])=currll1nchs;
                    end
                    PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcf=[...
                        PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcf;...
                        PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcf];
                    PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcfSh=[...
                        PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcfSh;...
                        PoolGCs.(gratCondTag{cnd}).(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcfSh];
                end
            end
        end
    end
end

%% PLOT PEN AVG GCs Layers x Layers (3x3 subplot)
if plotPenAvgGCs
    for vv1=[1 4]
        for vv2=[1 4]
            ylimplot=[0 0];
            for lltx=1:3
                for llrx=1:3
                    if ~isempty(PoolGCs.RF.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf)
                        nmRF=nanmean(PoolGCs.RF.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1);
                        nsRF=nansem(PoolGCs.RF.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1);
                        nmOUTrp=nanmean(PoolGCs.OUTrp.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1);
                        nsOUTrp=nansem(PoolGCs.OUTrp.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1);
                        ylimplot(1)=min([nmRF-nsRF nmOUTrp-nsOUTrp ylimplot(1)]);
                        ylimplot(2)=max([nmRF+nsRF nmOUTrp+nsOUTrp ylimplot(2)]);
                    end
                end
            end
            if (strcmpi(sSubjectName,'Taylor') && vv1==4 && vv2==1); ylimplot=2/3*ylimplot; end
            if vv1==vv2; ylimplot(1)=-1e-2; else ylimtoplot(1)=-1e-3; end
            
            hfig=figure('units','normalized','outerposition',[0 0 1 1]);
            for lltx=1:3
                for llrx=1:3
                    if ~isempty(PoolGCs.RF.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf)
                        for ff=1:length(freqAxisPreFirstDimDwSamp)
                            PoolGCs.pValues.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pvalues(ff)=...
                                signrank(PoolGCs.RF.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,ff),...
                                PoolGCs.OUTrp.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,ff));
                        end
                        pBarVec=PoolGCs.pValues.(['V' num2str(vv1) 'toV' num2str(vv2)]).([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.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcfSh(:,:,ff));
                                currLLsPwgcfRFShFreqWise(ff)=prctile(nanmean(currLLsFreqPwgcfRFSh,2),95);
                                currLLsFreqPwgcfOUTrpSh=squeeze(PoolGCs.OUTrp.(['V' num2str(vv1) 'toV' num2str(vv2)]).([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.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot)),'r');
                        plotmsem(freqAxisPreFirstDimDwSamp(freqIndices2Plot),(PoolGCs.OUTrp.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot)),'b');
                        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} ' V' num2str(vv1) ' to ' lamLayerTag{llrx} ' V' num2str(vv2)])
                        text(90, .9*(ylimplot(2)-ylimplot(1))+ylimplot(1),['tot n_' lamLayerTag{lltx}(1) '_,_' lamLayerTag{llrx}(1) '=' num2str(size(PoolGCs.RF.(['V' num2str(vv1) 'toV' num2str(vv2)]).([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf,1))],'color','k');
                        if llrx==3 && lltx==3
                            hleg=legend({'RF','OUT','p \leq .05','95% sh. perc. RF','95% sh. perc. OUT'},'Location','northeast');
                            hleg.Position=hleg.Position.*[.98 .91 1 1]; hleg.EdgeColor=[1 1 1];
                        end
                        if vv1~=vv2 && lltx==1
                            if llrx==1
                                text(90, .7*(ylimplot(2)-ylimplot(1))+ylimplot(1),['avg n_' lamLayerTag{llrx}(1) '_,_V_' num2str(vv2) '=' sprintf('%.2f',mean(PoolGCs.Labels.(['V' num2str(vv2)]).(['num' lamLayerTag{llrx}(1) 'chs'])))],'color','k');
                            else
                                text(90, .8*(ylimplot(2)-ylimplot(1))+ylimplot(1),['avg n_' lamLayerTag{llrx}(1) '_,_V_' num2str(vv2) '=' sprintf('%.2f',mean(PoolGCs.Labels.(['V' num2str(vv2)]).(['num' lamLayerTag{llrx}(1) 'chs'])))],'color','k');
                            end
                        end
                        if llrx==1
                            text(90, .8*(ylimplot(2)-ylimplot(1))+ylimplot(1),['avg n_' lamLayerTag{lltx}(1) '_,_V_' num2str(vv1) '=' sprintf('%.2f',mean(PoolGCs.Labels.(['V' num2str(vv1)]).(['num' lamLayerTag{lltx}(1) 'chs'])))],'color','k');
                        end
                    end
                end
            end
            if numInformChs==0
                supertitle([sSubjectName ' V' num2str(vv1) ' to V' num2str(vv2) ' Unconditional ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
            elseif numInformChs==Inf 
                supertitle([sSubjectName ' V' num2str(vv1) ' to V' num2str(vv2) ' Conditional (to all Complementary Channels) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
            elseif numInformChs>0
                supertitle([sSubjectName ' V' num2str(vv1) ' to V' num2str(vv2) ' 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 ' V' num2str(vv1) ' to V' num2str(vv2) ' Conditional (to Complementary Chs outside current Layers) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']);
            else
                supertitle([sSubjectName ' V' num2str(vv1) ' to V' num2str(vv2) ' 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
    end
end


%% PLOT PEN AVG GCs Layers x Layers (3x3 subplot)
if plotPenAvgGCs
   hfig=figure('units','normalized','outerposition',[0 0 1 1]);
    for vv1=[1 4]
        for vv2=[1 4]
            if vv1==vv2; ylimplot=[-1e-3 0.12]; else; ylimplot=[-1e-3 0.016]; end
            
            if ~isempty(PoolGCs.RF.(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcf)
                for ff=1:length(freqAxisPreFirstDimDwSamp)
                    PoolGCs.pValues.(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pvalues(ff)=...
                        signrank(PoolGCs.RF.(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcf(:,ff),...
                        PoolGCs.OUTrp.(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcf(:,ff));
                end
                pBarVec=PoolGCs.pValues.(['V' num2str(vv1) 'toV' num2str(vv2)]).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.(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcfSh(:,:,ff));
                        currLLsPwgcfRFShFreqWise(ff)=prctile(nanmean(currLLsFreqPwgcfRFSh,2),95);
                        currLLsFreqPwgcfOUTrpSh=squeeze(PoolGCs.OUTrp.(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcfSh(:,:,ff));
                        currLLsPwgcfOUTrpShFreqWise(ff)=prctile(nanmean(currLLsFreqPwgcfOUTrpSh,2),95);
                    end
                end
                
                subplot(2,2,(ceil(vv1/2)-1)*2+ceil(vv2/2))
                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.(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcf(:,freqIndices2Plot)),'r');
                plotmsem(freqAxisPreFirstDimDwSamp(freqIndices2Plot),(PoolGCs.OUTrp.(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcf(:,freqIndices2Plot)),'b');
                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(['V' num2str(vv1) ' to V' num2str(vv2)])
                if vv1==1 && vv2==4
                    hleg=legend({'RF','OUT','p \leq .05','95% sh. perc. RF','95% sh. perc. OUT'},'Location','northeast');
                    hleg.Position=hleg.Position.*[.98 .91 1 1]; hleg.EdgeColor=[1 1 1];
                end
                text(90, .9*(ylimplot(2)-ylimplot(1))+ylimplot(1),['tot n=' num2str(size(PoolGCs.RF.(['V' num2str(vv1) 'toV' num2str(vv2)]).AllLayers.pwgcf,1))],'color','k');

            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
    end
end