2024_Hasegawa_Huang_et_al_CODE/Fig4/Fig4.m

addpath(genpath('...\auxillary'));


load('...\network\DataAll.mat')
load('...\network\ZLowCa.mat')

SampleRate=30;
StimDur=2;%second
DelayDur=2;
ResponseWind=2;
NormBaseline=1;%second before stimOn 1s
AfterStimDur=2;%second after stimOff 2s
StageName={'Naive','Exp','RevNaive','RevExp'};
nTrainStage=length(StageName);% 'Naive','Exp','RevNaive','RevExp'
BatchName={'Aud1','Vis1'};
PeriodName={'Stim','Delay','Resp'};


if nTrainStage==4
    TaskCode{1,1}=[...
        1 1 12720;%Aud-Go   1
        1 1 180;  %Vis-Nogo 2
        1 2 12720;%Aud-Go   3
        1 2 180;  %Vis-Nogo 4
        2 3 12720;%Aud-Nogo 5
        2 3 180;  %Vis-Go   6
        2 4 12720;%Aud-Nogo 7
        2 4 180;  %Vis-Go   8
        ];
    % batch 2
    TaskCode{1,2}=[...
        2 1 12720;%Aud-NoGo 1
        2 1 180;  %Vis-Go   2
        2 2 12720;%Aud-NoGo 3
        2 2 180;  %Vis-Go   4
        1 3 12720;%Aud-Go   5
        1 3 180;  %Vis-NoGo 6
        1 4 12720;%Aud-Go   7
        1 4 180;  %Vis-NoGo 8
        ];
    
    GoCode{1,1}=TaskCode{1,1}([1 3 6 8],:);   %1
    GoCode{1,2}=TaskCode{1,2}([2 4 5 7],:);   %2
    NoGoCode{1,1}=TaskCode{1,1}([2 4 5 7],:); %3
    NoGoCode{1,2}=TaskCode{1,2}([1 3 6 8],:); %4
end



for iBatch=1:2
    for iSub=1:size(DataAll.(BatchName{iBatch}),1)
        for iTrainStage=1:nTrainStage%FWD REV
            tic
            for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                clear temp*
                temp=DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}){iSession};
                % pad nan before DLC vector to align rasberry pi cam
                tempLick=cat(1,nan(floor(abs(temp.CamStamp(1)-temp.RecStart)*SampleRate),1),temp.DLC_Lick);
                
                % find 2 consecutive frames with tongue deteected
                [~,LickIndDLC{1,iBatch}{iSub,iTrainStage}{iSession,1}]=...
                    findpeaks(tempLick,'MinPeakHeight',.5,'MinPeakWidth',2);
                
                % time vector is different in ScanStamp and CamStamp in rasberry pi cam
                t{1,iBatch}{iSub,iTrainStage}{iSession,1}=temp.ScanStamp;
                tCam{1,iBatch}{iSub,iTrainStage}{iSession,1}=temp.CamStamp;
                
                StimOnT{1,iBatch}{iSub,iTrainStage}{iSession,1}=temp.StimStart;
                RewardOnT{1,iBatch}{iSub,iTrainStage}{iSession,1}=temp.RewardOn;
                StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=dsearchn(t{1,iBatch}{iSub,iTrainStage}{iSession,1},StimOnT{1,iBatch}{iSub,iTrainStage}{iSession,1});
                StimOnIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1}=dsearchn(tCam{1,iBatch}{iSub,iTrainStage}{iSession,1},StimOnT{1,iBatch}{iSub,iTrainStage}{iSession,1});
                StimOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}+StimDur*SampleRate;
                StimOffIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1}=StimOnIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1}+StimDur*SampleRate;
                DelayOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=StimOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1};
                DelayOnIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1}=StimOffIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1};
                DelayOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=DelayOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}+DelayDur*SampleRate;
                DelayOffIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1}=DelayOnIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1}+DelayDur*SampleRate;
                ResponseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=DelayOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1};
                ResponseOnIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1}=DelayOffIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1};
                ResponseOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=ResponseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}+ResponseWind*SampleRate;
                ResponseOffIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1}=ResponseOnIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1}+ResponseWind*SampleRate;
                StimBaseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}-NormBaseline*SampleRate;
                RewardOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=dsearchn(t{1,iBatch}{iSub,iTrainStage}{iSession,1},RewardOnT{1,iBatch}{iSub,iTrainStage}{iSession,1});
                
                % electrical lick detector
                LickIndResp{1,iBatch}{iSub,iTrainStage}{iSession,1}=unique(dsearchn(t{1,iBatch}{iSub,iTrainStage}{iSession,1},temp.LickPulse));
                temp1=[ResponseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1} ResponseOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}];
                for iTrial=1:size(temp1,1)
                    if any(LickIndResp{1,iBatch}{iSub,iTrainStage}{iSession,1}>=temp1(iTrial,1)&...
                            LickIndResp{1,iBatch}{iSub,iTrainStage}{iSession,1}<=temp1(iTrial,2))
                        IfLickResp{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial,1)=true;
                    else
                        IfLickResp{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial,1)=false;
                    end
                end
                
                % DLC lick detector
                temp1=[DelayOnIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1} DelayOffIndCam{1,iBatch}{iSub,iTrainStage}{iSession,1}];
                for iTrial=1:size(temp1,1)
                    if any(LickIndDLC{1,iBatch}{iSub,iTrainStage}{iSession,1}>=temp1(iTrial,1)&...
                            LickIndDLC{1,iBatch}{iSub,iTrainStage}{iSession,1}<=temp1(iTrial,2))
                        IfLickDelayDLC{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial,1)=true;
                    else
                        IfLickDelayDLC{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial,1)=false;
                    end
                end
                
                
                tempIndHit=zeros(size(ZLowCa{1,iBatch}{iSub,iTrainStage}{iSession,1},1),1);
                for iTrial=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}){iSession}.StimStart)
                    tempIndHit(StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        ResponseOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),1)=1;
                end
                ITIindex{1,iBatch}{iSub,iTrainStage}{iSession,1}=find(tempIndHit==0);
                
                
                % If TaskType == 1, that means that Aud-Go and Vis-Nogo.
                % If TaskType == 2, that means that Vis-Go and Aud-Nogo.
                StimTypeTemp{1,iBatch}{iSub,iTrainStage}{iSession,1}=...
                    temp.Trial_Freq+temp.Trial_Grat;
                % AudTrial%12720 aud
                % VisTrial%180: vis
                % BlankTrial%720: vis&aud off blank
                for iTrial=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}){iSession}.StimStart)
                    TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial,1)=...
                        temp.TaskType;%TaskType
                    TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial,2)=...
                        iTrainStage;%TrainStage
                    TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial,3)=...
                        StimTypeTemp{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial);
                end
                
                % find Go trial in stim trial
                temp2=[ismember(TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1},...
                    GoCode{1,iBatch}(iTrainStage,:),'rows') ...
                    IfLickResp{1,iBatch}{iSub,iTrainStage}{iSession,1}];
                % Hit 1
                TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}...
                    (ismember(temp2,[1 1],'rows'),1)=1;
                % Miss 2
                TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}...
                    (ismember(temp2,[1 0],'rows'),1)=2;
                
                % find NoGo trial in stim trial
                temp2=[ismember(TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1},...
                    NoGoCode{1,iBatch}(iTrainStage,:),'rows') ...
                    IfLickResp{1,iBatch}{iSub,iTrainStage}{iSession,1}];
                
                % Correct Rejection in stim trial 3
                TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}...
                    (ismember(temp2,[1 0],'rows'),1)=3;
                % False Alarm in stim trial 4
                TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}...
                    (ismember(temp2,[1 1],'rows'),1)=4;
                
                % find Blank trial
                temp3=[TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1}(:,3)==720 ...
                    IfLickResp{1,iBatch}{iSub,iTrainStage}{iSession,1}];
                
                % CR in blank trial 5
                TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}...
                    (ismember(temp3,[1 0],'rows'),1)=5;
                % FA in blank trial 6
                TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}...
                    (ismember(temp3,[1 1],'rows'),1)=6;
                
                % HIT rate
                PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,1)=...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==1)/...
                    (sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==1)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==2));
                % MISS rate
                PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,2)=...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==2)/...
                    (sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==1)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==2));
                % CR rate (Stim &  Control trials)
                PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,3)=...
                    (sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==3)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==5))/...
                    (sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==3)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==4)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==5)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==6));
                % FA rate
                PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,4)=...
                    (sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==4)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==6))/...
                    (sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==3)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==4)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==5)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==6));
                % Blank CR rate
                PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,5)=...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==5)/...
                    (sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==5)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==6));
                % Blank FA rate
                PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,6)=...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==6)/...
                    (sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==5)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==6));
                % correcting dprime for HIT=1/0 or FA=1/0
                % 1/2N 
                temp=1/(2*(sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==1)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==2)));
                temp1=1/(2*(sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==3)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==4)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==5)+...
                    sum(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==6)));
                
                if PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,1)==1
                    PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,1)=...
                        PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,1)-temp;
                end
                if PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,1)==0
                    PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,1)=...
                        PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,1)+temp;
                end
                if PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,4)==1
                    PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,4)=...
                        PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,4)-temp1;
                end
                if PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,4)==0
                    PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,4)=...
                        PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,4)+temp1;
                end
                
                [dprime{1,iBatch}{iSub,iTrainStage}{iSession,1},~]=...
                    dprime_simple(PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,1),...
                    PerformanceAll{1,iBatch}{iSub,iTrainStage}{iSession,1}(1,4));
                clear temp*
            end
            toc
        end
    end
end
clear temp*


%% adjacency matrix (concatenated trials)
clear AdjMatAll temp*
for iBatch=1:2
    for iSub=1:size(DataAll.(BatchName{iBatch}),1)
        for iTrainStage=1:nTrainStage% INI EXP REV
            for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                % Go
                TrialTypeIdxTemp(ismember(TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1},...
                    GoCode{1,iBatch}(iTrainStage,:),'rows'),1)=1;
                % NoGo
                TrialTypeIdxTemp(ismember(TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1},...
                    NoGoCode{1,iBatch}(iTrainStage,:),'rows'),1)=2;
                % Blank
                TrialTypeIdxTemp(TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1}(:,3)==720,1)=3;
                [~,TrialTypeIndTemp]=sort(TrialTypeIdxTemp);
                
%                 temp=zscore(RawSpk{1,iBatch}{iSub,iTrainStage}{iSession,1});
                temp=ZLowCa{1,iBatch}{iSub,iTrainStage}{iSession,1};
                
                for iTrial=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}){iSession}.StimStart)
                    tempBL=mean(temp(StimBaseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:),1);
                    tempStim(:,:,iTrial)=temp(StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        StimOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:)-tempBL;
                    tempDelay(:,:,iTrial)=temp(DelayOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        DelayOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:)-tempBL;
                    tempResp(:,:,iTrial)=temp(ResponseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        ResponseOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:)-tempBL;
                    tempStimDelay(:,:,iTrial)=temp(StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        DelayOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:)-tempBL;
                end
                
                TrialTypeInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=1:sum(TrialTypeIdxTemp==1);
                TrialTypeInd{1,iBatch}{iSub,iTrainStage}{iSession,2}=[1:sum(TrialTypeIdxTemp==2)]+...
                    TrialTypeInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(end);
                TrialTypeInd{1,iBatch}{iSub,iTrainStage}{iSession,3}=[1:sum(TrialTypeIdxTemp==3)]+...
                    TrialTypeInd{1,iBatch}{iSub,iTrainStage}{iSession,2}(end);
                
                % concatenate all trials time series 
                for iType=1:3%Go/NoGo/catch
                    tempStim1=tempStim(:,:,TrialTypeIndTemp(...
                        TrialTypeInd{1,iBatch}{iSub,iTrainStage}{iSession,iType}));
                    tempDelay1=tempDelay(:,:,TrialTypeIndTemp(...
                        TrialTypeInd{1,iBatch}{iSub,iTrainStage}{iSession,iType}));
                    tempResp1=tempResp(:,:,TrialTypeIndTemp(...
                        TrialTypeInd{1,iBatch}{iSub,iTrainStage}{iSession,iType}));
                    tempStimDelay1=tempStimDelay(:,:,TrialTypeIndTemp(...
                        TrialTypeInd{1,iBatch}{iSub,iTrainStage}{iSession,iType}));
                    
                    tempStim1=zscore(reshape(permute(tempStim1,[1 3 2]),size(tempStim1,1)*size(tempStim1,3),size(tempStim1,2)));
                    tempDelay1=zscore(reshape(permute(tempDelay1,[1 3 2]),size(tempDelay1,1)*size(tempDelay1,3),size(tempDelay1,2)));
                    tempResp1=zscore(reshape(permute(tempResp1,[1 3 2]),size(tempResp1,1)*size(tempResp1,3),size(tempResp1,2)));
                    tempStimDelay1=zscore(reshape(permute(tempStimDelay1,[1 3 2]),size(tempStimDelay1,1)*size(tempStimDelay1,3),size(tempStimDelay1,2)));
                    
                    tempN=size(tempStim1,1)-1;
                    AdjMatStim(:,:,iType)=tempStim1'*tempStim1/tempN;
                    AdjMatDelay(:,:,iType)=tempDelay1'*tempDelay1/tempN;
                    AdjMatResp(:,:,iType)=tempResp1'*tempResp1/tempN;
                    AdjMatStimDelay(:,:,iType)=tempStimDelay1'*tempStimDelay1/tempN;
                end
                
                AdjMatAll{1,iBatch}{iSub,iTrainStage}{iSession,1}=AdjMatStim;
                AdjMatAll{1,iBatch}{iSub,iTrainStage}{iSession,2}=AdjMatDelay;
                AdjMatAll{1,iBatch}{iSub,iTrainStage}{iSession,3}=AdjMatResp;
                AdjMatAll{1,iBatch}{iSub,iTrainStage}{iSession,4}=AdjMatStimDelay;
                
                clear temp* AdjMatStim AdjMatDelay AdjMatResp AdjMatStimDelay TrialTypeIdxTemp
            end
        end
    end
end




%% clustering coef, Geodesic path length,co-firing strength
% contatenate all trials to one time series

SignFactor=[1 -1];% split corr mat to positive R and negative R
for iBatch=1:2
    for iSub=1:size(DataAll.(BatchName{iBatch}),1)
        tic
        for iTrainStage=1:nTrainStage% INI EXP REV
            for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                for iPeriod=1:4 %BL Stim Delay Resp StimDelay(concaternated time series from stim to delay)
                    for iType=1:3%Go/NoGo/catch
                        for iSign=1:2
                            % not fixed with signed weight CluCoef
                            % clustering coef
                            temp=SignFactor(iSign)*AdjMatAll{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,:,iType);
                            temp(logical(diag(ones(1,size(temp,1)),0)))=0;
                            temp(temp<0)=0;
                            tempMax=max(temp,[],'all');% clu coef normalisation
                            % exclude negative weight
                            AdjMatAllNorm=temp/tempMax;
                            CluCoef{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,iType,iSign)=...
                                clustering_coef_wu(AdjMatAllNorm);
                            
                            % co-firing strength (Hubness)
                            temp=SignFactor(iSign)*AdjMatAll{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,:,iType);
                            temp(find(temp<0))=0;
                            Weight2All{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,iType,iSign)=...
                                mean(temp,2);
                            
                            % Geodesic path length
                            temp=SignFactor(iSign)*AdjMatAll{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,:,iType);
                            temp1=1./temp;
                            temp1(temp1<0)=inf;
                            temp1(logical(diag(ones(1,size(temp,1)),0)))=0;
                            temp2=distance_wei(temp1);
                            temp3=temp2(find(triu(temp2)~=0));
                            PathLength{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,iType,iSign)=...
                                temp3;
                            clear temp*
                        end
                    end
                end
            end
        end
        toc
    end
end

clear temp* AdjMatAllNorm CluCoefPerMouse Weight2AllPerMouse PathLengthPerMouse identifiers identifiers1
TypeNameTemp={'G','N','C'};
PeriodNameTemp={'S','D','R','SD'};
StageNameTemp={'1','2','3','4'};
RValueName={'p','n'};

for iBatch=1:2
    for iSub=1:size(DataAll.(BatchName{iBatch}),1)
        for iTrainStage=1:4% INI EXP REV-naive REV-expert
            tempClu=[];tempWeight=[];tempPath=[];tempLabel=[];tempLabel1=[];
            for iPeriod=1 %Stim Delay Resp
                for iType=[2 1]%:2%Go/NoGo/catch
                    for iSign=1% +R/-R
                        for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                            
                            label=[StageNameTemp{iTrainStage},'-',PeriodNameTemp{iPeriod},'-',...
                                TypeNameTemp{iType},'-',RValueName{iSign}];
                            
                            temp=CluCoef{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,iType,iSign);
                            tempClu=cat(1,tempClu,temp);
                            temp1=Weight2All{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,iType,iSign);
                            tempWeight=cat(1,tempWeight,temp1);
                            tempLabel=cat(1,tempLabel,string(repmat(label,length(temp),1)));
                            temp2=PathLength{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,iType,iSign);
                            tempPath=cat(1,tempPath,temp2);
                            tempLabel1=cat(1,tempLabel1,string(repmat(label,length(temp2),1)));
                        end
                    end
                end
            end
            CluCoefPerMouse{1,iBatch}{iSub,iTrainStage}=tempClu;
            Weight2AllPerMouse{1,iBatch}{iSub,iTrainStage}=tempWeight;
            PathLengthPerMouse{1,iBatch}{iSub,iTrainStage}=tempPath;
            identifiers{1,iBatch}{iSub,iTrainStage}=tempLabel;
            identifiers1{1,iBatch}{iSub,iTrainStage}=tempLabel1;
            clear temp*
        end
    end
end


clear temp* Group* label*
GroupClu=[];GroupWeight=[];GroupPath=[];GroupLabel=[];GroupLabel1=[];GroupLabel2=[];
for iBatch=1:2
    for iSub=1:size(DataAll.(BatchName{iBatch}),1)
        
        % pool mice
        tempDim=numel(CluCoefPerMouse{1,iBatch});
        data=cell2mat(reshape(CluCoefPerMouse{1,iBatch},tempDim,1));
        data1=cell2mat(reshape(Weight2AllPerMouse{1,iBatch},tempDim,1));
        data2=cell2mat(reshape(PathLengthPerMouse{1,iBatch},tempDim,1));
        temp=reshape(identifiers{1,iBatch},tempDim,1);
        temp1=reshape(identifiers1{1,iBatch},tempDim,1);
        label=cellstr(cat(1,temp{:}));
        label1=cellstr(cat(1,temp{:}));
        label2=cellstr(cat(1,temp1{:}));
        
        GroupClu=cat(1,GroupClu,data);
        GroupWeight=cat(1,GroupWeight,data1);
        GroupPath=cat(1,GroupPath,data2);
        GroupLabel=cat(1,GroupLabel,label);
        GroupLabel1=cat(1,GroupLabel1,label1);
        GroupLabel2=cat(1,GroupLabel2,label2);
        
    end
end

% labels

%1-2-3-4:INI EXP REV-naive REV-expert
%S-D:Stimulus period, Delay period
%G-N:Go-Nogo trial type
%p-n:positive-R, negative-R

%example
%1-S-N-p: Naive-Stimulus-Nogo-positiveR
%4-D-G-n: REV-expert-Delay-Go-negativeR

% only check fig1
% plot path length
close all
figure;
FscatJit2(GroupLabel2,GroupPath,[],'N',1);
suplabel('Path','x',[.5 1 0 0]);%top x

% plot Hubness
close all
figure
FscatJit2(GroupLabel1,GroupWeight,[],'N',5);
suplabel('Hubness','x',[.5 .99 0 0]);%top x

close all
figure
FscatJit2(GroupLabel,GroupClu,[],'N',5);
suplabel('CluCoef','x',[.5 .99 0 0]);%top x

xtickangle(45)

% circular plot
tempSession=[1 3 1 2];
iBatch=1;iSub=3;iPeriod=2;

for iType=1:2 % Go,Nogo,Catch
    figure
    i=1;
    for iTrainStage=1:nTrainStage
        iSession=tempSession(iTrainStage);
        temp=AdjMatAll{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,:,iType);
        subplot(2,4,i)
        imagesc(temp);
        title(['C',num2str(round(mean(CluCoef{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,iType)),2))...
            ' L',num2str(round(mean(PathLength{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,iType)),2))...
            ' H',num2str(round(mean(Weight2All{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(:,iType)),3))])
        axis square
        caxis([0 .3])
        
        subplot(2,4,i+4)
        temp(find(temp<0.3|temp>.99))=0;
        G=graph(temp,'omitselfloops','upper');
        plot(G,'Layout','circle','EdgeColor','r','NodeLabel',{},...
            'LineWidth',5*G.Edges.Weight);%,'WeightEffect','inverse'
        title(["d'=",num2str(dprime{1,iBatch}{iSub,iTrainStage}{iSession,1})])
        axis square
        i=i+1;
        
    end
end