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 () [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 () %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