addpath(genpath('./support_routines')); addpath(genpath('./support_tools')); clear all; close all; clc; sSubjectName='monkey 2'; sVisualArea='V4'; reComputeGCs=0; reComputeGCsShuffling=0; plotSinglePenGCs=0; plotPenAvgGCs=1; reportSignifThreshold=1; plotPenAvgGCsSingChRes=0; lamLayerTag={'Supra','Granr','Infra'}; gratCondTag={'RF','OUTrp','OUT1','OUT2'}; spectWinTag={'Full','Alpha','Beta','Gamma'}; %% cGC implementation settings numInformChs=-2; % numInformChs = [-Inf -2 -1 0 1 2 Inf]; % number of informative channels used as regressors in the cGC strategy % 0 = Unconditional GC % >0 = cGC to most informative channels % <0 = cGC to m. inf. outside laminar compartments of each X,Y channel pair % Inf = cGC to ALL channles (ignores most informative computation) compMethod='Mvgc'; % compMethod = 'Mvgc','Geweke','MatPart'; choose method to compute cGC % Mvgc = (Barnett & Seth, J. Neurosci. Methods, 2014) % Geweke = (Geweke, J. of the American Stat. Assoc., 1984) % MatPart = (Chen et al., J. Neurosci Methods 2006) %% DATA SAMPLING SETTINGS SAMPLING_FREQ=1017.375; NUM_TIME_SAMPLES_EXTRACTED=1024; timeAxis0=linspace(-NUM_TIME_SAMPLES_EXTRACTED/SAMPLING_FREQ,0,NUM_TIME_SAMPLES_EXTRACTED)*1e3; NUM_TIME_SAMPLES=512; timeAxisPreFirstDim=timeAxis0(NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end); freqAxisPreFirstDim=linspace(0,SAMPLING_FREQ/2,NUM_TIME_SAMPLES); DOWNSAMPLING_FACTOR=4; NUM_TIME_DWSAMPLES=NUM_TIME_SAMPLES/DOWNSAMPLING_FACTOR; DWSAMPLING_FREQ=SAMPLING_FREQ/DOWNSAMPLING_FACTOR; timeAxisPreFirstDimDwSamp=linspace(timeAxisPreFirstDim(1),timeAxisPreFirstDim(end),NUM_TIME_DWSAMPLES); freqAxisPreFirstDimDwSamp=linspace(0,DWSAMPLING_FREQ/2,NUM_TIME_DWSAMPLES+1); freqIndices2Plot=find(freqAxisPreFirstDimDwSamp>2); if reComputeGCs subjectInfos=getSubjectInfos(sSubjectName); subjectData=getSubjectData(sSubjectName,sVisualArea,'ALL','LFPs','ALL'); % REMOVE DATA IN TIME WINDOWS OTHER THAN PREFIRSTDIM TO SAVE MEMORY for pp=1:length(subjectData.LfpStruct) if ~isempty(subjectData.LfpStruct(pp).Sorted) rmfn=fieldnames(subjectData.LfpStruct(pp).Sorted.Full.Supra.RF); rmfn(cellfun(@(fnm) strcmpi(fnm,'PreFirstDimDataBiZSc'),rmfn))=[]; for ll=1:3 subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF=rmfield(subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).RF,rmfn); subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1=rmfield(subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1,rmfn); subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2=rmfield(subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2,rmfn); end end end %% RANDPERM FOR AVG SUBSAMPLING OF TRIALS AT OUTSIDE LOCATIONS rng(0); % will pick exactly the same trials at every run, % Note: varying this seed might give slightly different cGC magnitudes % (different trials are selected) but results are robust to permutation for pp=1:length(subjectData.LfpStruct) if ~isempty(subjectData.LfpStruct(pp).Sorted) currNtrOUT1=size(subjectData.LfpStruct(pp).Sorted.Full.Supra.OUT1.PreFirstDimDataBiZSc,2); currRandPermTrials=randperm(2*currNtrOUT1,currNtrOUT1); for ll=1:3 currOUTcat=cat(2,subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT1.PreFirstDimDataBiZSc,... subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUT2.PreFirstDimDataBiZSc); subjectData.LfpStruct(pp).Sorted.Full.(lamLayerTag{ll}).OUTrp.PreFirstDimDataBiZSc=currOUTcat(:,currRandPermTrials,:); end end end eval(clearoutmem('curr*')); %% Information Toolbox settings (support_tools) - DOI: 10.1186/1471-2202-10-81 numQuantBins=4; MImethod='dr'; MIbias='qe'; MIbtsp=0; bandWinsWidth=2; bandWinsStart=.01; % Starts at 1 Hz bandWinsEnd=100; % Stops at ~100 Hz bandWinsEnd=min(bandWinsEnd,SAMPLING_FREQ/2); NUM_BAND_WINS=floor((bandWinsEnd-bandWinsStart)/bandWinsWidth); bandWinsCutoffs=[linspace(bandWinsStart,bandWinsEnd-bandWinsWidth,NUM_BAND_WINS);... linspace(bandWinsStart,bandWinsEnd-bandWinsWidth,NUM_BAND_WINS)+bandWinsWidth]'; freqAxisBandWins=mean(bandWinsCutoffs,2); filterOrder=2; filterBCoeffs=nan(NUM_BAND_WINS,2*filterOrder+1); filterACoeffs=nan(NUM_BAND_WINS,2*filterOrder+1); filterFreqResp=nan(NUM_BAND_WINS,512); filterImpResp=nan(NUM_BAND_WINS,2048); filterFreqAxis=nan(1,512); filterTimeAxis=nan(1,512); for bw=1:NUM_BAND_WINS [filterBCoeffs(bw,:),filterACoeffs(bw,:)]=butter(filterOrder,bandWinsCutoffs(bw,:)./(SAMPLING_FREQ/2),'bandpass'); [filterFreqResp(bw,:),filterFreqAxis]=freqz(filterBCoeffs(bw,:),filterACoeffs(bw,:),512,SAMPLING_FREQ); [filterImpResp(bw,:),filterTimeAxis]=impz(filterBCoeffs(bw,:),filterACoeffs(bw,:),2048,SAMPLING_FREQ); end fresEstim=1/4*NUM_TIME_DWSAMPLES; lambdafreqAxis=linspace(0,freqAxisPreFirstDimDwSamp(end),fresEstim+1); %% MVGC Toolbox settings (support_tools) - DOI: 10.1016/j.jneumeth.2013.10.018 regmode='LWR'; morder=10; % set up by optimization of AIC and R² (see Suppl. Figure S7) acmaxlags=[]; % autocovariance max lags acdectol=[]; % autocovariance decay tolerance (sets memory of the process) NShuffles=100; % number of shuffles used to compute significance threshold parfor pp=1:length(subjectData.LfpStruct) if ~isempty(subjectData.LfpStruct(pp).Sorted) for cnd=1:2 disp(['Processing cGC for ' sSubjectName ' ' sVisualArea ' PEN ' num2str(subjectData.LfpStruct(pp).penNum) ', ' gratCondTag{cnd}]); currX0=[subjectData.LfpStruct(pp).Sorted.Full.Supra.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end);... subjectData.LfpStruct(pp).Sorted.Full.Granr.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end);... subjectData.LfpStruct(pp).Sorted.Full.Infra.(gratCondTag{cnd}).PreFirstDimDataBiZSc(:,:,NUM_TIME_SAMPLES_EXTRACTED-NUM_TIME_SAMPLES+1:end)]; currXp=permute(currX0,[1 3 2]); %eval(clearoutmem('currX0')); currSchs=1:size(subjectData.LfpStruct(pp).Sorted.Full.Supra.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1); currGchs=max([currSchs 0])+(1:size(subjectData.LfpStruct(pp).Sorted.Full.Granr.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1)); currIchs=max([currSchs currGchs 0])+(1:size(subjectData.LfpStruct(pp).Sorted.Full.Infra.(gratCondTag{cnd}).PreFirstDimDataBiZSc,1)); % very awful but required by parfor transparency constraints currSGIchs=nan(3,16); % set to 16, higher than every max_{j=I,G,S}{length(currjchs)}, allows to cat along 1st dim. currSGIchs(1,1:length(currSchs))=currSchs; currSGIchs(2,1:length(currGchs))=currGchs; currSGIchs(3,1:length(currIchs))=currIchs; currSGIalgns=nan(3,16); currSGIalgns(1,1:length(currSchs))=subjectData.LfpStruct(pp).Sorted.Labels.SupraChsLamAlignsBi; currSGIalgns(2,1:length(currGchs))=subjectData.LfpStruct(pp).Sorted.Labels.GranrChsLamAlignsBi; currSGIalgns(3,1:length(currIchs))=subjectData.LfpStruct(pp).Sorted.Labels.InfraChsLamAlignsBi; % De-trend along trial-dimension may help for stationarity aspects currXDetr=mvdetrend(currXp); % De-mean along time-dimension currXDemn=currXDetr-repmat(mean(currXDetr,2),[1 NUM_TIME_SAMPLES 1]); % Downsampling by movmean currXDemnp=permute(currXDemn,[2 1 3]); currXDownp=movmeandecimatrix(currXDemnp,DOWNSAMPLING_FACTOR); currXDwSamp=permute(currXDownp,[2 1 3]); % Downsampling by brute decimation %currXDwSamp=currXDemn(:,1:DOWNSAMPLING_FACTOR:end,:); [currNvars,currNobs,currNtrials]=size(currXDwSamp); currfpw=nan(currNvars,currNvars,NUM_TIME_DWSAMPLES+1); currfpwSh=nan(NShuffles,currNvars,currNvars,NUM_TIME_DWSAMPLES+1); currFFiltHilbPowDown=nan(currNvars,NUM_BAND_WINS,currNtrials,NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR); % BINNING AND INFORMATION EVALUATION paramsMI=[]; paramsMI.method=MImethod; paramsMI.bias=MIbias; paramsMI.btsp=MIbtsp; currMIMat=nan(NUM_BAND_WINS,NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR,currNvars,currNvars); for bw=1:NUM_BAND_WINS % LOOP OVER BAND WINDOWS for cch=1:currNvars % Filter for each band window currBwinFFilt=filtfilt(filterBCoeffs(bw,:),filterACoeffs(bw,:),squeeze(currXDwSamp(cch,:,:))); % Hilbert Transform and abs^2 currBwinFFiltHilbPow=(abs(hilbertpad(currBwinFFilt)).^2); % Downsampling by moving average [ch x NUM_BAND_WINS x trials x NUM_TIME_DWSAMAPLES] currFFiltHilbPowDown(cch,bw,:,:)=movmeandecimatrix(currBwinFFiltHilbPow,DOWNSAMPLING_FACTOR)'; end % [currNvars x NBWINS x currNtrials x NTIMEDWS] currBwXDwMat=reshape(vectorisen(currFFiltHilbPowDown(:,bw,:,:),[2 3]),[1 currNvars*currNtrials*NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR]); currChsLabels=reshape(repmat(1:currNvars,currNtrials*NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR,1)',[1 currNvars*currNtrials*NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR]); % [1 x (currNtrials*NTIMEDWS) x currNvars] [currRespMatAllChs, currnTrAllChs]= buildr(currChsLabels,currBwXDwMat); currRespMatBinAllChs = binr(currRespMatAllChs, currnTrAllChs, numQuantBins, 'eqpop'); currRespMatBinAllChsRes=reshape(currRespMatBinAllChs(1,:,:),[currNtrials NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR currNvars]); paramsMI.nt=currNtrials;%currnTrAllChs(1); for tt=1:(NUM_TIME_DWSAMPLES/DOWNSAMPLING_FACTOR) for ch1=1:currNvars for ch2=ch1:currNvars currMIMat(bw,tt,ch1,ch2)=information(reshape(currRespMatBinAllChsRes(:,tt,[ch1 ch2]),[1 currNtrials 2]), paramsMI, 'I'); end end end end currTimeCumulInformMat=squeeze(nansum(vectorisen(currMIMat,[1 2]),2));%squeeze(nansum(currMIMat,1)); currTimeCumulInformMat2=currTimeCumulInformMat+triu(currTimeCumulInformMat',-1); currSrtMIMat=nan(currNvars,currNvars); currSrtMIChs=nan(currNvars,currNvars); currSrtRangeNormMIMat=nan(currNvars,currNvars); for ch=1:currNvars [currSrtMIMat(ch,:),currSrtMIChs(ch,:)]=sort(currTimeCumulInformMat2(ch,:),'descend'); end for chY=1:currNvars for chX=1:currNvars if chX~=chY %required by parfor, ch2=ch1:nvars is not allowed chsZOut=[]; if ~(any(chX==currSchs) || any(chY==currSchs)) chsZOut=currSchs; end if ~(any(chX==currGchs) || any(chY==currGchs)) chsZOut=[chsZOut currGchs]; end if ~(any(chX==currIchs) || any(chY==currIchs)) chsZOut=[chsZOut currIchs]; end if numInformChs==0 chsZ=[]; elseif numInformChs>0 chsZMIwrtX=currSrtMIChs(chX,1:min(numInformChs,end)); chsZMIwrtY=currSrtMIChs(chY,1:min(numInformChs,end)); chsZ=setdiff(unique([chsZMIwrtX,chsZMIwrtY],'stable'),[chX chY]); else if numInformChs==-inf chsZ=chsZOut; else numInformChsAbs=abs(numInformChs); chsZMIwrtX=intersect(currSrtMIChs(chX,1:min(numInformChsAbs,end)),chsZOut,'stable'); chsZMIwrtY=intersect(currSrtMIChs(chY,1:min(numInformChsAbs,end)),chsZOut,'stable'); chsZ=unique([chsZMIwrtX,chsZMIwrtY],'stable'); end end numZvars=length(chsZ); disp([sSubjectName '-' sVisualArea '-PEN-' num2str(subjectData.LfpStruct(pp).penNum) '-' gratCondTag{cnd} '-' compMethod '-cGC_( ' num2str(chY) ' -> ' num2str(chX) ' \ ' num2str(chsZ) ')']); % FIT A VAR(p) MODEL FOR (X,Y,Z) - (Barnett & Seth, J. Neurosci. Methods, 2014) [BxyzMat,ExyzMat]=tsdata_to_var(currXDwSamp([chX chY chsZ],:,:),morder,regmode); assert(~isbad(BxyzMat),'VAR estimation failed'); if strcmpi(compMethod,'Mvgc') && (chX < chY) % Autocovariance calculation () [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_21 = autocov_to_smvgc(GammaMat,2,1,fresEstim); % from 1 to 2 f_multichs_fres_12 = autocov_to_smvgc(GammaMat,1,2,fresEstim); % from 2 to 1 % assert(~isbad(f_multichs_fres,false),'Mvgc spectral GC calculation failed'); % if fres=[] the freqAxis ends up being non-uniform % currch1ch2freqAxis=linspace(0,DWSAMPLING_FREQ/2,size(multichsfpw,3)); % so interp1 is used to set a uniform output fres f_y_to_x_given_z=interp1(lambdafreqAxis,squeeze(f_multichs_fres_12),freqAxisPreFirstDimDwSamp,'spline'); f_x_to_y_given_z=interp1(lambdafreqAxis,squeeze(f_multichs_fres_21),freqAxisPreFirstDimDwSamp,'spline'); currfpw(chX,chY,:)=abs(f_y_to_x_given_z); currfpw(chY,chX,:)=abs(f_x_to_y_given_z); if reComputeGCsShuffling %SHUFFLE TRIALS for nnsh=1:NShuffles currXDwSampShuf=nan(size(currXDwSamp([chX chY chsZ],:,:))); for chshuf=1:length([chX chY chsZ]) currXDwSampShuf(chshuf,:,:)=currXDwSamp(chshuf,:,randperm(currNtrials)); end [BxyzMatShuf,ExyzMatShuf]=tsdata_to_var(currXDwSampShuf,morder,regmode); assert(~isbad(BxyzMatShuf),'VAR estimation failed'); [GammaMatShuf,GammaInfoShuf] = var_to_autocov(BxyzMatShuf,ExyzMatShuf,acmaxlags,acdectol); %f_multichs_fres_Sh = autocov_to_spwcgc(GammaMatShuf,fresEstim); f_multichs_fres_Sh_21 = autocov_to_smvgc(GammaMatShuf,2,1,fresEstim); f_multichs_fres_Sh_12 = autocov_to_smvgc(GammaMatShuf,1,2,fresEstim); %assert(~isbad(f_multichs_fres_Sh,false),'Mvgc spectral GC calculation failed'); f_y_to_x_given_z_sh=interp1(lambdafreqAxis,squeeze(f_multichs_fres_Sh_12),freqAxisPreFirstDimDwSamp,'spline'); f_x_to_y_given_z_sh=interp1(lambdafreqAxis,squeeze(f_multichs_fres_Sh_21),freqAxisPreFirstDimDwSamp,'spline'); currfpwSh(nnsh,chX,chY,:)=abs(f_y_to_x_given_z_sh); currfpwSh(nnsh,chY,chX,:)=abs(f_x_to_y_given_z_sh); % max(currfpwSh(:,chX,chY,:),[],4) is to be used for prctile end end end if strcmpi(compMethod,'Geweke') % ORIGINAL GEWEKE SPECTRAL FACTORIZATION METHOD % (Geweke, J. of the American Stat. Assoc., 1984) % FIT A SECOND, REDUCED VAR(p) MODEL FOR (X,Z) [DxzMat,ExzMat]=tsdata_to_var(currXDwSamp([chX chsZ],:,:),morder,regmode); assert(~isbad(DxzMat),'VAR estimation failed'); % P MATRICES FOR COVARIANCE NORMALIZATION PxzMat=[1 zeros(1,numZvars); -ExzMat(1,2:end)'/(ExzMat(1,1)) eye(numZvars)]; PxyzMat=[1 0 zeros(1,numZvars); -ExyzMat(1,2)/(ExyzMat(1,1)) 1 zeros(1,numZvars); -ExyzMat(1,3:end)'/(ExyzMat(1,1)) zeros(numZvars,1) eye(numZvars)]*... [1 0 zeros(1,numZvars); 0 1 zeros(1,numZvars); zeros(numZvars,1) -(ExyzMat(3:end,2) - ExyzMat(3:end,1)/ExyzMat(1,1)*ExyzMat(1,2))/(ExyzMat(2,2) - ExyzMat(2,1)/ExyzMat(1,1)*ExyzMat(1,2)) eye(numZvars)]; % APPLYING NORMALIZATION BY P-MATRICES MULTIPLICATION IN TIME PExzMat=PxzMat*ExzMat*PxzMat'; PExyzMat=PxyzMat*ExyzMat*PxyzMat'; PDxzMat=nan(numZvars+1,numZvars+1,morder); PBxyzMat=nan(numZvars+2,numZvars+2,morder); for mm=1:morder PDxzMat(:,:,mm)=PxzMat*DxzMat(:,:,mm); PBxyzMat(:,:,mm)=PxyzMat*BxyzMat(:,:,mm); end [PSxzMat,PGxzMat]=var_to_cpsd(PDxzMat,PExzMat,fresEstim); [PSxyzMat,PHxyzMat]=var_to_cpsd(PBxyzMat,PExyzMat,fresEstim); % APPLYING NORMALIZATION BY P-MATRICES MULTIPLICATION IN FREQUENCY %[SxzMat,GxzMat]=var_to_cpsd(DxzMat,ExzMat,fresEstim); %[SxyzMat,HxyzMat]=var_to_cpsd(BxyzMat,ExyzMat,fresEstim); f_y_to_x_given_z_fres=nan(1,fresEstim+1); for ll=1:fresEstim+1 % NORMALIZED TRANSFER FUNCTIONS P*G(X,Z) and P2*P1*H(X,Y,Z) %PGxzMat=GxzMat(:,:,lambda)/(PxzMat); %G*inv(P) %PHxyzMat=HxyzMat(:,:,lambda)/(PxyzMat); %H*inv(P) % EXPANDED G MATRIX FOR SPECTRAL FACTORIZATION GxyzMat=[PGxzMat(1,1,ll) 0 PGxzMat(1,2:end,ll); 0 1 zeros(1,numZvars); PGxzMat(2:end,1,ll) zeros(numZvars,1) PGxzMat(2:end,2:end,ll)]; QxyzMat=GxyzMat\PHxyzMat(:,:,ll); %inv(G)*H %SIGxyz=abs(QxyzMat(1,1)*ExyzMat(1,1)*conj(QxyzMat(1,1)) + QxyzMat(1,2)*ExyzMat(2,2)*conj(QxyzMat(1,2)) + QxyzMat(1,3:end)*ExyzMat(3:end,3:end)*conj(QxyzMat(1,3:end))'); SIGxyz=abs(QxyzMat(1,1)*PExyzMat(1,1)*conj(QxyzMat(1,1)) + QxyzMat(1,2)*PExyzMat(2,2)*conj(QxyzMat(1,2)) + QxyzMat(1,3:end)*PExyzMat(3:end,3:end)*conj(QxyzMat(1,3:end))'); %SIGx=abs(QxyzMat(1,1)*ExyzMat(1,1)*conj(QxyzMat(1,1))); SIGx=abs(QxyzMat(1,1)*PExyzMat(1,1)*conj(QxyzMat(1,1))); if QxyzMat(1,1)*PExyzMat(1,1)*conj(QxyzMat(1,1)) <0 warning('Qxx product is negative'); end if QxyzMat(1,2)*PExyzMat(2,2)*conj(QxyzMat(1,2)) <0 warning('Qxy product is negative'); end if QxyzMat(1,3:end)*PExyzMat(3:end,3:end)*conj(QxyzMat(1,3:end)') <0 warning('Qxz product is negative'); end f_y_to_x_given_z_fres(ll)= log( SIGxyz ./ SIGx ); %f_y_to_x_given_z_fres(lambda)= log( abs(PExzMat(1,1)) ./ SIGx ); end f_y_to_x_given_z=interp1(lambdafreqAxis,f_y_to_x_given_z_fres,freqAxisPreFirstDimDwSamp,'spline'); currfpw(chX,chY,:)=f_y_to_x_given_z; end if strcmpi(compMethod,'MatPart') % MATRIX PARTITIONING METHOD (Chen et al., J. Neurosci Methods 2006) % P MATRIX FOR COVARIANCE NORMALIZATION PxyzMat=[1 0 zeros(1,numZvars); -ExyzMat(1,2)/(ExyzMat(1,1)) 1 zeros(1,numZvars); -ExyzMat(1,3:end)'/(ExyzMat(1,1)) zeros(numZvars,1) eye(numZvars)]*... [1 0 zeros(1,numZvars); 0 1 zeros(1,numZvars); zeros(numZvars,1) -(ExyzMat(3:end,2) - ExyzMat(3:end,1)/ExyzMat(1,1)*ExyzMat(1,2))/(ExyzMat(2,2) - ExyzMat(2,1)/ExyzMat(1,1)*ExyzMat(1,2)) eye(numZvars)]; % APPLYING NORMALIZATION BY P-MATRICES MULTIPLICATION IN TIME PExyzMat=PxyzMat*ExyzMat*PxyzMat'; PBxyzMat=nan(numZvars+2,numZvars+2,morder); for mm=1:morder PBxyzMat(:,:,mm)=PxyzMat*BxyzMat(:,:,mm); end [PSxyzMat,PHxyzMat]=var_to_cpsd(PBxyzMat,PExyzMat,fresEstim); f_y_to_x_given_z_fres=nan(1,fresEstim+1); for ll=1:fresEstim+1 %PHxyzMat=HxyzMat(:,:,lambda)/(PxyzMat); %H*inv(P) % MIXED TRANSFER FUNCTIONS HBARxy_zy=[HBARxy; HBARzy] HBAR_xy_zy=([PHxyzMat(1,1,ll) PHxyzMat(1,3:end,ll); PHxyzMat(3:end,1,ll) PHxyzMat(3:end,3:end,ll)]) \ [PHxyzMat(1,2,ll); PHxyzMat(3:end,2,ll)]; % REDUCED COVARIANCE MATRIX SIGBAR_xz=[SIGBARxx SIGBARxz; SIGBARzx SIGBARzz] EBAR_xz=[PExyzMat(1,1) PExyzMat(1,3:end); PExyzMat(3:end,1) PExyzMat(3:end,3:end)]+HBAR_xy_zy*[PExyzMat(1,2)' PExyzMat(3:end,2)']+[PExyzMat(1,2)' PExyzMat(3:end,2)']'*conj(HBAR_xy_zy')+PExyzMat(2,2)*HBAR_xy_zy*conj(HBAR_xy_zy'); % REDUCED P MATRIX FOR COVARIANCE NORMALIZATION PBARMat_xz=[1 zeros(1,numZvars); -EBAR_xz(2:end,1)/EBAR_xz(1,1) eye(numZvars)]; GMat_xz=[PHxyzMat(1,1,ll) PHxyzMat(1,3:end,ll); PHxyzMat(3:end,1,ll) PHxyzMat(3:end,3:end,ll)]/(PBARMat_xz); % G/P = G*inv(P) % EXPANDED G MATRIX FOR SPECTRAL FACTORIZATION GMat_xyz=[GMat_xz(1,1) 0 GMat_xz(1,2:end); 0 1 zeros(1,numZvars); GMat_xz(2:end,1) zeros(numZvars,1) GMat_xz(2:end,2:end)]; Q_xyz=(GMat_xyz)\PHxyzMat(:,:,ll); % G\H = inv(G)*H %SIGxyz=abs( Q_xyz(1,1) * ExyzMat(1,1) * conj(Q_xyz(1,1)) + Q_xyz(1,2) * ExyzMat(2,2) * conj(Q_xyz(1,2)) + Q_xyz(1,3:end) * ExyzMat(3:end,3:end) * conj(Q_xyz(1,3:end))'); SIGxyz=abs( Q_xyz(1,1) * PExyzMat(1,1) * conj(Q_xyz(1,1)) + Q_xyz(1,2) * PExyzMat(2,2) * conj(Q_xyz(1,2)) + Q_xyz(1,3:end) * PExyzMat(3:end,3:end) * conj(Q_xyz(1,3:end))'); %SIGx=abs( Q_xyz(1,1) * ExyzMat(1,1) * conj(Q_xyz(1,1)) ); SIGx=abs( Q_xyz(1,1) * PExyzMat(1,1) * conj(Q_xyz(1,1)) ); f_y_to_x_given_z_fres(ll)=log( SIGxyz ./ SIGx ); %f_y_to_x_given_z_fres(lambda)=log( abs(EBAR_xz(1,1)) ./ SIGx ); end f_y_to_x_given_z=interp1(lambdafreqAxis,f_y_to_x_given_z_fres,freqAxisPreFirstDimDwSamp,'spline'); currfpw(chX,chY,:)=f_y_to_x_given_z; end end end end currGCsStruct(pp).(gratCondTag{cnd}).pwgcf=currfpw; currGCsStruct(pp).(gratCondTag{cnd}).pwgcfSh=currfpwSh; end currGCsStruct(pp).penNum=subjectData.LfpStruct(pp).penNum; currGCsStruct(pp).currSGIchs=currSGIchs; currGCsStruct(pp).currSGIalgns=currSGIalgns; end end % SORTING BY VISUAL AREA / LAMINAR LAYER / GRATC CONDITION for pp=1:length(currGCsStruct) if ~isempty(currGCsStruct(pp).currSGIchs) currSGIchs=currGCsStruct(pp).currSGIchs; currSGIalgns=currGCsStruct(pp).currSGIalgns; for cnd=1:2 for ll1=1:3 for ll2=1:3 % eval is awful but it saves from extra indexing currll1chs=currSGIchs(ll1,~isnan(currSGIchs(ll1,:))); currll2chs=currSGIchs(ll2,~isnan(currSGIchs(ll2,:))); % NOTE: reverse to get fi->j|[ij] instead of fj->i|[ij] currllsfpw=currGCsStruct(pp).(gratCondTag{cnd}).pwgcf(currll2chs,currll1chs,:); currllsfpwSh=currGCsStruct(pp).(gratCondTag{cnd}).pwgcfSh(:,currll2chs,currll1chs,:); if ~isempty(currllsfpw) % pool chsxchs within ll1xll2 currllsfpwRes=reshape(currllsfpw,size(currllsfpw,1)*size(currllsfpw,2),size(currllsfpw,3)); currllsfpwResSh=reshape(currllsfpwSh,NShuffles,size(currllsfpwSh,2)*size(currllsfpwSh,3),size(currllsfpwSh,4)); % get rid of all-NaNs rowsacdpt[] currllsfpwRes((sum(isnan(currllsfpwRes),2)==size(currllsfpwRes,2)),:)=[]; currllsfpwResSh(:,sum(isnan(currllsfpwResSh(1,:,:)),3)==size(currllsfpwResSh,3),:)=[]; sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcf=currllsfpwRes; sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcfSh=currllsfpwResSh; %eval(clearoutmem('currlls*')); else sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcf=[]; sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcfSh=[]; end currll1algns=currSGIalgns(ll1,~isnan(currSGIalgns(ll1,:))); currll2algns=currSGIalgns(ll2,~isnan(currSGIalgns(ll2,:))); if length(currll1chs)~=length(currll1algns) || length(currll2chs)~=length(currll2algns); error('Layer labels do not match alignments.'); end for algn1=1:length(currll1algns) for algn2=1:length(currll2algns) currllalgnsfpw=currGCsStruct(pp).(gratCondTag{cnd}).pwgcf(currll2chs(algn2),currll1chs(algn1),:); currllalgnsfpwRes=reshape(currllalgnsfpw,1,size(currllalgnsfpw,3)); currllalgnsfpwSh=currGCsStruct(pp).(gratCondTag{cnd}).pwgcfSh(:,currll2chs(algn2),currll1chs(algn1),:); currllalgnsfpwResSh=reshape(currllalgnsfpwSh,NShuffles,1,size(currllalgnsfpwSh,4)); if currll1algns(algn1)<0 strcurrll1algn1=['m' num2str(abs(currll1algns(algn1)))]; else strcurrll1algn1=['p' num2str(currll1algns(algn1))]; end if currll2algns(algn2)<0 strcurrll2algn2=['m' num2str(abs(currll2algns(algn2)))]; else strcurrll2algn2=['p' num2str(currll2algns(algn2))]; end if sum(isnan(currllalgnsfpw),2)~=size(currllalgnsfpw,2) sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcf'])=currllalgnsfpwRes; sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcfSh'])=currllalgnsfpwResSh; else sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcf'])=[]; sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([strcurrll1algn1 'to' strcurrll2algn2 'pwgcfSh'])=[]; end end end end if cnd==2 sGCsStruct(pp).penNum=currGCsStruct(pp).penNum; sGCsStruct(pp).Labels.(['num' lamLayerTag{ll1}(1) 'chs'])=length(currll1chs); sGCsStruct(pp).Labels.([lamLayerTag{ll1} 'ChsAlgns'])=currll1algns; end end end end end save(['./' upper(sSubjectName(1)) sSubjectName(end) sVisualArea '_cGCs_' num2str(numInformChs) '_MostInfChs_' compMethod '_PreFirstDim.mat'],'sGCsStruct','-v7.3'); else load(['./' upper(sSubjectName(1)) sSubjectName(end) sVisualArea '_cGCs_' num2str(numInformChs) '_MostInfChs_' compMethod '_PreFirstDim.mat']); end %% POOL ACROSS PENs GCs PoolGCs=[]; for cnd=1:2 PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcf=[]; PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcfSh=[]; for ll1=1:3 for ll2=1:3 currll1nchs=[]; PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}])=[]; for pp=1:length(sGCsStruct) if ~isempty(sGCsStruct(pp).RF) currSubFieldNames=fieldnames(sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}])); currSubFieldNames=currSubFieldNames(~contains(currSubFieldNames,'Sh')); for sbfs=1:length(currSubFieldNames) if ~isfield(PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]),currSubFieldNames{sbfs}) PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs})=[]; PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh'])=[]; end PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs})=... [PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs});... sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).(currSubFieldNames{sbfs})]; PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh'])=... cat(2,PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh']),... sGCsStruct(pp).(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).([currSubFieldNames{sbfs} 'Sh'])); end currll1nchs=cat(1,currll1nchs,sGCsStruct(pp).Labels.(['num' lamLayerTag{ll1}(1) 'chs'])); end end if cnd==2 PoolGCs.Labels.(['num' lamLayerTag{ll1}(1) 'chs'])=currll1nchs; end PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcf=[PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcf;... PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcf]; PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcfSh=[PoolGCs.(gratCondTag{cnd}).([lamLayerTag{ll1} 'to' lamLayerTag{ll2}]).pwgcfSh;... PoolGCs.(gratCondTag{cnd}).AllLayers.pwgcfSh]; end end end %% PLOT PEN AVG GCs Layers x Layers (3x3 subplot) if plotPenAvgGCs ylimplot=[0 0]; for lltx=1:3 for llrx=1:3 if ~isempty(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf) ylimplot(1)=min([nanmean(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)-... nansem(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)... nanmean(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)-... nansem(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1) ylimplot(1)]); ylimplot(2)=max([nanmean(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)+... nansem(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)... nanmean(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1)+... nansem(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot),1) ylimplot(2)]); end end end ylimplot=1*ylimplot; ylimplot(1)=-1e-2; hfig=figure('units','normalized','outerposition',[0 0 1 1]); for lltx=1:3 for llrx=1:3 if ~isempty(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf) for ff=1:length(freqAxisPreFirstDimDwSamp) PoolGCs.pValues.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pvalues(ff)=... signrank(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,ff),... PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,ff)); end pBarVec=PoolGCs.pValues.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pvalues(freqIndices2Plot); [~,~,~,pBarVecFDR]=fdr_bh(pBarVec); currLLsPwgcfRFShFreqWise=nan(1,NUM_TIME_DWSAMPLES+1); currLLsPwgcfOUTrpShFreqWise=nan(1,NUM_TIME_DWSAMPLES+1); if reportSignifThreshold for ff=1:(NUM_TIME_DWSAMPLES+1) currLLsFreqPwgcfRFSh=squeeze(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcfSh(:,:,ff)); currLLsPwgcfRFShFreqWise(ff)=prctile(nanmean(currLLsFreqPwgcfRFSh,2),95); currLLsFreqPwgcfOUTrpSh=squeeze(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcfSh(:,:,ff)); currLLsPwgcfOUTrpShFreqWise(ff)=prctile(nanmean(currLLsFreqPwgcfOUTrpSh,2),95); end end subplot(3,3,(lltx-1)*3+llrx) plotcmapdots([1 0 0; 0 0 1;.2 .6 .2]); plot([-inf -inf],[-inf -inf],'r:'); plot([-inf -inf],[-inf -inf],'b:'); plotmsem(freqAxisPreFirstDimDwSamp(freqIndices2Plot),(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot)),'r');%[.73 0 0]); plotmsem(freqAxisPreFirstDimDwSamp(freqIndices2Plot),(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf(:,freqIndices2Plot)),'b');%[0 .73 .73]); plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),currLLsPwgcfRFShFreqWise(freqIndices2Plot),'r:'); plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),currLLsPwgcfOUTrpShFreqWise(freqIndices2Plot),'b:'); plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),double(pBarVecFDR<=0.05)-1+ylimplot(1),'color',[.2 .6 .2],'linewidth',3) xlim([0 100]); ylim(ylimplot); title([lamLayerTag{lltx} ' to ' lamLayerTag{llrx}]) if llrx<=lltx text(100, .9*(ylimplot(2)-ylimplot(1))+ylimplot(1),['tot n=' num2str(size(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcf,1))],'color','k'); end if llrx==3 && lltx==1 legend({'RF','OUT','p \leq 0.05','95% sh. perc. RF','95% sh. perc. OUT'},'Location','northeast'); end if llrx==1 text(100, .8*(ylimplot(2)-ylimplot(1))+ylimplot(1),['avg n=' sprintf('%.2f',mean(PoolGCs.Labels.(['num' lamLayerTag{lltx}(1) 'chs'])))],'color','k'); end end end end if numInformChs==0 supertitle([sSubjectName ' ' sVisualArea ' Unconditional ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']); elseif numInformChs==Inf supertitle([sSubjectName ' ' sVisualArea ' Conditional (to all Complementary Channels) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']); elseif numInformChs>0 supertitle([sSubjectName ' ' sVisualArea ' Conditional (to n=' num2str(numInformChs) ' Most Inform. Complem. Chs) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']); elseif numInformChs==-Inf supertitle([sSubjectName ' ' sVisualArea ' Conditional (to Complementary Chs outside current Layers) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']); else supertitle([sSubjectName ' ' sVisualArea ' Conditional (to n=' num2str(abs(numInformChs)) ' Most Inf. Compl. Chs outside current Layers) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']); end end %% PLOT PEN AVG GCs Layers x Layers (3x3 subplot) if plotPenAvgGCs ylimplot=[-1e-2 0.12]; hfig=figure('units','normalized','outerposition',[0 0 1 1]); for ff=1:length(freqAxisPreFirstDimDwSamp) PoolGCs.pValues.AllLayers.pvalues(ff)=signrank(PoolGCs.RF.AllLayers.pwgcf(:,ff),PoolGCs.OUTrp.AllLayers.pwgcf(:,ff)); end pBarVec=PoolGCs.pValues.AllLayers.pvalues(freqIndices2Plot); [~,~,~,pBarVecFDR]=fdr_bh(pBarVec); currLLsPwgcfRFShFreqWise=nan(1,NUM_TIME_DWSAMPLES+1); currLLsPwgcfOUTrpShFreqWise=nan(1,NUM_TIME_DWSAMPLES+1); if reportSignifThreshold for ff=1:(NUM_TIME_DWSAMPLES+1) currLLsFreqPwgcfRFSh=squeeze(PoolGCs.RF.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcfSh(:,:,ff)); currLLsPwgcfRFShFreqWise(ff)=prctile(nanmean(currLLsFreqPwgcfRFSh,2),95); currLLsFreqPwgcfOUTrpSh=squeeze(PoolGCs.OUTrp.([lamLayerTag{lltx} 'to' lamLayerTag{llrx}]).pwgcfSh(:,:,ff)); currLLsPwgcfOUTrpShFreqWise(ff)=prctile(nanmean(currLLsFreqPwgcfOUTrpSh,2),95); end end plotcmapdots([1 0 0; 0 0 1;.2 .6 .2]); plot([-inf -inf],[-inf -inf],'r:'); plot([-inf -inf],[-inf -inf],'b:'); plotmsem(freqAxisPreFirstDimDwSamp(freqIndices2Plot),(PoolGCs.RF.AllLayers.pwgcf(:,freqIndices2Plot)),'r');%[.73 0 0]); plotmsem(freqAxisPreFirstDimDwSamp(freqIndices2Plot),(PoolGCs.OUTrp.AllLayers.pwgcf(:,freqIndices2Plot)),'b');%[0 .73 .73]); plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),currLLsPwgcfRFShFreqWise(freqIndices2Plot),'r:'); plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),currLLsPwgcfOUTrpShFreqWise(freqIndices2Plot),'b:'); plot(freqAxisPreFirstDimDwSamp(freqIndices2Plot),double(pBarVecFDR<=0.05)-1+ylimplot(1),'color',[.2 .6 .2],'linewidth',3) xlim([0 100]); ylim(ylimplot); text(90, .9*(ylimplot(2)-ylimplot(1))+ylimplot(1),['tot n=' num2str(size(PoolGCs.RF.AllLayers.pwgcf,1))],'color','k'); if llrx==3 && lltx==1 legend({'RF','OUT','p \leq 0.05','95% sh. perc. RF','95% sh. perc. OUT'},'Location','northeast'); end if numInformChs==0 supertitle([sSubjectName ' ' sVisualArea ' Unconditional ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']); elseif numInformChs==Inf supertitle([sSubjectName ' ' sVisualArea ' Conditional (to all Complementary Channels) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']); elseif numInformChs>0 supertitle([sSubjectName ' ' sVisualArea ' Conditional (to n=' num2str(numInformChs) ' Most Inform. Complem. Chs) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']); elseif numInformChs==-Inf supertitle([sSubjectName ' ' sVisualArea ' Conditional (to Complementary Chs outside current Layers) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']); else supertitle([sSubjectName ' ' sVisualArea ' Conditional (to n=' num2str(abs(numInformChs)) ' Most Inf. Compl. Chs outside current Layers) ' compMethod ' GC - Pre First Dim LFPs [-500,0 ms] - PEN AVG']); end end