2024_Hasegawa_Huang_et_al_CODE/Fig3/Fig3.m

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

load('...\PVC\DataAll.mat')
load('...\PVC\ZLowCa.mat')
set(0,'DefaultFigureColormap',colormap(redblue(256)))


SampleRate=30;
StimDur=2;%second
DelayDur=2;
ResponseWind=2;
NormBaseline=1;%second before stimOn 1s
AfterStimDur=2;%second after stimOff 2s
ButFiltBand=[5];
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*


%% plot all correlation matrices correlation matrix stim delay resp separate

for iPeriod=1:2
    i=1;
    figure
    for iBatch=1:2
        for iSub=1:size(DataAll.(BatchName{iBatch}),1)
            for iTrainStage=1:nTrainStage
                for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                    % Go-NoGo-blank
                    clear TrialTypeIdx
                    % Go
                    TrialTypeIdx(ismember(TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1},...
                        GoCode{1,iBatch}(iTrainStage,:),'rows'),1)=1;
                    % NoGo
                    TrialTypeIdx(ismember(TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1},...
                        NoGoCode{1,iBatch}(iTrainStage,:),'rows'),1)=2;
                    % Blank
                    TrialTypeIdx(TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1}(:,3)==720,1)=3;
                    
                    [~,TrialTypeInd]=sort(TrialTypeIdx);
                    
                    CorrBLInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=...
                        StimBaseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(TrialTypeInd,1);
                    CorrStimInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=...
                        StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(TrialTypeInd,1);
                    CorrDelayInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=...
                        DelayOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(TrialTypeInd,1);
                    CorrRespInd{1,iBatch}{iSub,iTrainStage}{iSession,1}=...
                        ResponseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(TrialTypeInd,1);
                    
                    CorrResponseWind=2;
                    
                    temp=ZLowCa{1,iBatch}{iSub,iTrainStage}{iSession,1};
                    
                    tempMat1=[];tempMat2=[];tempMat3=[];
                    for iTrial=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}){iSession}.StimStart)
                        f0=mean(temp(CorrBLInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                            CorrStimInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:),1);
                        %stim
                        temp1=mean(temp(CorrStimInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                            CorrStimInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial)+StimDur*SampleRate,:),1)-f0;
                        %delay
                        temp2=mean(temp(CorrDelayInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                            CorrDelayInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial)+DelayDur*SampleRate,:),1)-f0;
                        %resp
                        temp3=mean(temp(CorrRespInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                            CorrRespInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial)+CorrResponseWind*SampleRate,:),1)-f0;
                        
                        tempMat1=cat(1,tempMat1,temp1);
                        tempMat2=cat(1,tempMat2,temp2);
                        tempMat3=cat(1,tempMat3,temp3);
                    end
                    
                    tempMat1=zscore(tempMat1')';
                    tempMat2=zscore(tempMat2')';
                    tempMat3=zscore(tempMat3')';
                    CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,1}=...
                        tempMat1*tempMat1'/(size(tempMat1,2)-1);%stim
                    CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,2}=...
                        tempMat2*tempMat2'/(size(tempMat2,2)-1);%delay
                    CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,3}=...
                        tempMat3*tempMat3'/(size(tempMat3,2)-1);%resp
                    
                    subplot(8,15,(i-1)*15+(iTrainStage-1)*4+iSession)
%                     TimeAxis=[1:sum(TrialTypeIdx==1) 1:sum(TrialTypeIdx==2) 1:sum(TrialTypeIdx==3)];
                    
                    temp=sum(TrialTypeIdx==1)+sum(TrialTypeIdx==2);
                    imagesc(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(...
                        1:temp,1:temp))
                    
                    % be careful of 1 more index number
                    Lines(cumsum([sum(TrialTypeIdx==1)]),[],'k','-',2);
                    Lines([],cumsum([sum(TrialTypeIdx==1)]),'k','-',2);
%                     colorbar
                    caxis([-.3 .6])
                    xticks([1 temp])
                    yticks([1 temp])
                    colormap(bluewhitered(256))
                    axis square
                    %                 title([BatchName{iBatch},'-m',num2str(iSub),'-',StageName{iTrainStage},...
                    %                     '-s',num2str(iSession),'d=',num2str(dprime{1,iBatch}{iSub,iTrainStage}{iSession,1})])
                    clear temp*
                end
            end
            i=i+1;
        end
    end
end
suplabel(PeriodName{iPeriod},'x',[.5 .99 0 0]);%top x


%% PV correlation Go/Nogo grouping

for iPeriod=1:2% stim delay resp
    for iBatch=1:2
        for iSub=1:size(DataAll.(BatchName{iBatch}),1)
            for iTrainStage=1:nTrainStage
                for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                    
                    clear TrialTypeIdx
                    % Go
                    TrialTypeIdx(ismember(TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1},...
                        GoCode{1,iBatch}(iTrainStage,:),'rows'),1)=1;
                    % NoGo
                    TrialTypeIdx(ismember(TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1},...
                        NoGoCode{1,iBatch}(iTrainStage,:),'rows'),1)=2;
                    % Blank
                    TrialTypeIdx(TrialType{1,iBatch}{iSub,iTrainStage}{iSession,1}(:,3)==720,1)=3;
                    
                    temp1=1:sum(TrialTypeIdx==1);
                    temp2=[1:sum(TrialTypeIdx==2)]+temp1(end);
                    temp3=[1:sum(TrialTypeIdx==3)]+temp2(end);
                    
                    CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}...
                        (logical(diag(ones(1,size(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod},1)),0)))=nan;
                    
                    % Go-Go corr
                    DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(1,1)=...
                        nanmean(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp1,temp1),'all');
                    DesignCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{1,1}=...
                        rmmissing(reshape(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp1,temp1),[],1));
%                     DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{1,1}(1,2:3)=...
%                         prctile(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp1,temp1),percentile,'all');
                    % NoGo-NoGo corr
                    DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(2,2)=...
                        nanmean(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp2,temp2),'all');
                    DesignCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{2,2}=...
                        rmmissing(reshape(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp2,temp2),[],1));
%                     DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{2,2}(1,2:3)=...
%                         prctile(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp2,temp2),percentile,'all');
                    % catch-catch corr
                    DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(3,3)=...
                        nanmean(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp3,temp3),'all');
                    DesignCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{3,3}=...
                        rmmissing(reshape(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp3,temp3),[],1));
%                     DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{3,3}(1,2:3)=...
%                         prctile(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp3,temp3),percentile,'all');
                    % Go-NoGo corr
                    DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(2,1)=...
                        nanmean(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp2,temp1),'all');
                    DesignCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{2,1}=...
                        rmmissing(reshape(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp2,temp1),[],1));
%                     DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{2,1}(1,2:3)=...
%                         prctile(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp2,temp1),percentile,'all');
                    % Go-catch corr
                    DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(3,1)=...
                        nanmean(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp3,temp1),'all');
                    DesignCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{3,1}=...
                        rmmissing(reshape(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp3,temp1),[],1));
%                     DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{3,1}(1,2:3)=...
%                         prctile(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp3,temp1),percentile,'all');
                    % NoGo-catch corr
                    DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(3,2)=...
                        nanmean(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp3,temp2),'all');
                    DesignCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{3,2}=...
                        rmmissing(reshape(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp3,temp2),[],1));
%                     DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}{3,2}(1,2:3)=...
%                         prctile(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp3,temp2),percentile,'all');
                    clear temp*
                end
            end
        end
    end
end


clear tmep*
for iPeriod=1:2%stim delay resp
    for iResult=1:2% Go-Go/NoGo-NoGo/catch-catch
        for iTrainStage=1:nTrainStage
            CatDesignCorr=[];CatDesignCorr1=[];
            for iBatch=1:2
                for iSub=1:size(DataAll.(BatchName{iBatch}),1)
                    for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                        temp(iSession,1)=...
                            DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(iResult,iResult);
                    end
                    CatDesignCorr=cat(1,CatDesignCorr,mean(temp));
                    CatDesignCorr1=cat(1,CatDesignCorr1,temp);
                    clear temp*
                end
            end
            DesignCorrPerMouse(:,iTrainStage,iPeriod,iResult)=CatDesignCorr;
            DesignCorrPerMouseSessions{iTrainStage,iPeriod,iResult}=CatDesignCorr1;
        end
    end
end

% change iTrainStage to 1:2 for comparing Initial learning stage or 3:4 for
% reversal learning stage
% plot individual session in each training phase
clear Cat* label
ResultName={'Go','NoGo'};
label=[];CatDesignCorrPerMouseSession=[];
for iPeriod=1:2%stim delay resp
    for iResult=1:2% Go-Go/NoGo-NoGo/catch-catch
        for iTrainStage=3:4% change here
            CatDesignCorrPerMouseSession=cat(1,CatDesignCorrPerMouseSession,...
                DesignCorrPerMouseSessions{iTrainStage,iPeriod,iResult});
            label=cat(1,label,string(repmat([StageName{iTrainStage},'-',PeriodName{iPeriod},'-',ResultName{iResult}],...
                length(DesignCorrPerMouseSessions{iTrainStage,iPeriod,iResult}),1)));
        end
    end
end

figure
[~,ViolinOutput]=violinplot(CatDesignCorrPerMouseSession,label,'ViolinColor',...
    repmat([hex2rgb('#D35400');hex2rgb('#2980B9')],4,1));
xtickangle(45)

% plot Aud or Vis Go first group

close all
BatchGroup=[12 12 11 11;13 13 12 12];% Aud Vis train sequence and bumber of session
figure
hold on
bar(cellfun(@mean,ViolinOutput(:,1)),'w','BarWidth',.5,'LineWidth',2);
errorbar([1:8],cellfun(@mean,ViolinOutput(:,1)),sem(ViolinOutput(:,1),0),'Color',[0 0 0],'LineStyle','none','LineWidth',2,'CapSize',10)
for i=1:8
    switch ViolinOutput{i,2}(1:4)
        case 'Naiv'
            scatter(ones(length(ViolinOutput{i,1}(1:BatchGroup(1,1))),1)*i,...
                ViolinOutput{i,1}(1:BatchGroup(1,1)),'r','filled','jitter','on','jitterAmount',0.2);
            scatter(ones(length(ViolinOutput{i,1}(BatchGroup(2,1):end)),1)*i,...
                ViolinOutput{i,1}(BatchGroup(2,1):end),'g','filled','jitter','on','jitterAmount',0.2);
        case 'Exp-'
            scatter(ones(length(ViolinOutput{i,1}(1:BatchGroup(1,2))),1)*i,...
                ViolinOutput{i,1}(1:BatchGroup(1,2)),'r','filled','jitter','on','jitterAmount',0.2);
            scatter(ones(length(ViolinOutput{i,1}(BatchGroup(2,2):end)),1)*i,...
                ViolinOutput{i,1}(BatchGroup(2,2):end),'g','filled','jitter','on','jitterAmount',0.2);
        case 'RevN'
            scatter(ones(length(ViolinOutput{i,1}(1:BatchGroup(1,3))),1)*i,...
                ViolinOutput{i,1}(1:BatchGroup(1,3)),'r','filled','jitter','on','jitterAmount',0.2);
            scatter(ones(length(ViolinOutput{i,1}(BatchGroup(2,3):end)),1)*i,...
                ViolinOutput{i,1}(BatchGroup(2,3):end),'g','filled','jitter','on','jitterAmount',0.2);
        case 'RevE'
            scatter(ones(length(ViolinOutput{i,1}(1:BatchGroup(1,4))),1)*i,...
                ViolinOutput{i,1}(1:BatchGroup(1,4)),'r','filled','jitter','on','jitterAmount',0.2);
            scatter(ones(length(ViolinOutput{i,1}(BatchGroup(2,4):end)),1)*i,...
                ViolinOutput{i,1}(BatchGroup(2,4):end),'g','filled','jitter','on','jitterAmount',0.2);
    end
end

temp=gcf;
temp.Parent.Children.CurrentAxes.XTick=[1:8];
temp.Parent.Children.CurrentAxes.XTickLabel=ViolinOutput(:,2);
ylim([0 .5])
temp.Parent.Children.CurrentAxes.YTick=0:.1:.5;
xtickangle(45)


%% statistics

clear TableDesignMeanCorrSession
TableDesignMeanCorrSession=table('Size',[0 5],'VariableTypes',{'double','string','string','string','string'},...
    'VariableNames',{'R','stage','period','design','ID'});
% TableDesignCorrTrial=table('Size',[0 5],'VariableTypes',{'double','string','string','string','string'},...
%     'VariableNames',{'R','stage','period','design','ID'});
tempMouseID{1,1}=1:4;tempMouseID{1,2}=5:8;
ResultName={'Go','NoGo'};

% change iPeriod iResult iTrainStage for subgroup statistical analysis
i=1;
for iPeriod=2%stim delay resp
    for iResult=1% Go-Go/NoGo-NoGo/catch-catch
        for iTrainStage=1:2
            for iBatch=1:2
                for iSub=1:size(DataAll.(BatchName{iBatch}),1)
                    for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                        % within session mean (2-3 sessions) as observation
                        TableDesignMeanCorrSession(i,:)={...
                            DesignMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(iResult,iResult)...
                            ,StageName{iTrainStage},PeriodName{iPeriod},ResultName{iResult},...
                            string(tempMouseID{1,iBatch}(iSub))};
                        i=i+1;
                    end
                end
            end
        end
    end
end

LMM=fitlme(TableDesignMeanCorrSession(:,[1 2 5]),'R~stage+(1|ID)')


%% sort trial result Hit/FA

%% sort trial result Hit/FA
for iPeriod=1:2
    i=1;
    figure
    for iBatch=1:2
        for iSub=1:size(DataAll.(BatchName{iBatch}),1)
            for iTrainStage=1:nTrainStage%FWD REV
                for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                    temp=ZLowCa{1,iBatch}{iSub,iTrainStage}{iSession,1};
                    tempMat1=[];tempMat2=[];tempMat3=[];
                    for iTrial=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}){iSession}.StimStart)
                        f0=mean(temp(StimBaseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                            StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:),1);
                        %stim
                        temp1=mean(temp(StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                            StimOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:))-f0;
                        %delay
                        temp2=mean(temp(DelayOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                            DelayOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:))-f0;
                        %resp
                        temp3=mean(temp(ResponseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                            ResponseOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:))-f0;
                        
                        tempMat1=cat(1,tempMat1,temp1);
                        tempMat2=cat(1,tempMat2,temp2);
                        tempMat3=cat(1,tempMat3,temp3);
                    end
                    
                    
                    % [Hit; FA]
                    TrialResultIdx{1,1}=find(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==1);
                    TrialResultIdx{1,2}=find(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==4);
                    tempMat1=zscore([tempMat1(TrialResultIdx{1,1},:);tempMat1(TrialResultIdx{1,2},:)]')';
                    tempMat2=zscore([tempMat2(TrialResultIdx{1,1},:);tempMat2(TrialResultIdx{1,2},:)]')';
                    tempMat3=zscore([tempMat3(TrialResultIdx{1,1},:);tempMat3(TrialResultIdx{1,2},:)]')';
                    
                    CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,1}=...
                        tempMat1*tempMat1'/(size(tempMat1,2)-1);%stim
                    CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,2}=...
                        tempMat2*tempMat2'/(size(tempMat2,2)-1);%delay
                    CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,3}=...
                        tempMat3*tempMat3'/(size(tempMat3,2)-1);%resp
                    
                    subplot(8,15,(i-1)*15+(iTrainStage-1)*4+iSession)
                    imagesc(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod})
                    % be careful of 1 more index number
                    
                    Lines(length(TrialResultIdx{1,1}),[],'k','',1.5);
                    Lines([],length(TrialResultIdx{1,1}),'k','',1.5);
                    axis square
                    caxis([-.3 .6])
                    xticks([1 sum(cellfun(@length,TrialResultIdx))])
                    yticks([1 sum(cellfun(@length,TrialResultIdx))])
                    colormap(bluewhitered(256))
                end
            end
            i=i+1;
        end
    end
    colorbar
end


% summary data
TrialNumThr=6;
for iBatch=1:2
    for iSub=1:size(DataAll.(BatchName{iBatch}),1)
        for iTrainStage=1:nTrainStage%FWD REV
            for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                clear temp*
                temp=ZLowCa{1,iBatch}{iSub,iTrainStage}{iSession,1};
                tempMat1=[];tempMat2=[];tempMat3=[];
                for iTrial=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}){iSession}.StimStart)
                    f0=mean(temp(StimBaseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:),1);
                    %stim
                    temp1=mean(temp(StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        StimOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:))-f0;
                    %delay
                    temp2=mean(temp(DelayOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        DelayOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:))-f0;
                    %resp
                    temp3=mean(temp(ResponseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        ResponseOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:))-f0;
                    
                    tempMat1=cat(1,tempMat1,temp1);
                    tempMat2=cat(1,tempMat2,temp2);
                    tempMat3=cat(1,tempMat3,temp3);
                end
                
                % [Hit;FA]%CR
                TrialResultIdx{1,1}=find(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==1);
                TrialResultIdx{1,2}=find(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==4);
                % TrialResultIdx{1,3}=find(TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1}==3);

                tempMat1=zscore([tempMat1(TrialResultIdx{1,1},:);tempMat1(TrialResultIdx{1,2},:)]')';%;tempMat1(TrialResultIdx{1,3},:)
                tempMat2=zscore([tempMat2(TrialResultIdx{1,1},:);tempMat2(TrialResultIdx{1,2},:)]')';%;tempMat2(TrialResultIdx{1,3},:)
                tempMat3=zscore([tempMat3(TrialResultIdx{1,1},:);tempMat3(TrialResultIdx{1,2},:)]')';%;tempMat3(TrialResultIdx{1,3},:)
                
                CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,1}=...
                    tempMat1*tempMat1'/(size(tempMat1,2)-1);%stim
                CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,2}=...
                    tempMat2*tempMat2'/(size(tempMat2,2)-1);%delay
                CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,3}=...
                    tempMat3*tempMat3'/(size(tempMat3,2)-1);%resp
                
                ResultInd=cumsum([length(TrialResultIdx{1,1}) length(TrialResultIdx{1,2})]);% length(TrialResultIdx{1,2}) 
                temp1=[1:ResultInd(1)];
                temp2=[ResultInd(1):ResultInd(2)];
                for iPeriod=1:3  %stim delay resp
                    % Hit corr
                    ResultMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(1,1)=...
                        mean(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp1,temp1),'all');
                    % FA corr
                    ResultMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(2,2)=...
                        mean(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp2,temp2),'all');
                    % CR corr
                    % ResultMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(3,3)=...
                    % mean(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp3,temp3),'all');
                    
                    % Avoid bias of small num of trial (<TrialNumThr will be discard)
                    % If there is no trial being in Hit/FA/CR, diagonal "1" will be selected
                    
                    if length(temp1)<TrialNumThr
                        ResultMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(1,1)=nan;
                    end
                    if length(temp2)<TrialNumThr
                        ResultMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(2,2)=nan;
                    end
%                     if length(temp3)<TrialNumThr
%                         ResultMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(3,3)=1;
%                     end
                    
                end
            end
        end
    end
end


clear Cat* temp* ResultCorrPerMouse ResultCorrPerMouseSessions
for iPeriod=1:2%stim delay resp
    for iResult=1:2% Hit/FA 
        for iTrainStage=1:nTrainStage
            CatResultCorr=[];CatResultCorr1=[];
            for iBatch=1:2
                for iSub=1:size(DataAll.(BatchName{iBatch}),1)
                    i=1;
                    for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                        if ~isnan(ResultMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(iResult,iResult))
                            temp(i,1)=ResultMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(iResult,iResult);
                            i=i+1;
                        end
                    end
                    CatResultCorr=cat(1,CatResultCorr,nanmean(temp));
                    CatResultCorr1=cat(1,CatResultCorr1,temp);
                    clear temp
                end
            end
            ResultCorrPerMouse(:,iTrainStage,iPeriod,iResult)=CatResultCorr;
            ResultCorrPerMouseSessions{iTrainStage,iPeriod,iResult}=CatResultCorr1;
        end
    end
end



% statistics for individual session data in suppl.
% change iTrainStage to plot Initial and Reversal training stage
clear Cat* label
ResultName={'Hit','FA'};
label=[];CatResultCorrPerMouseSession=[];
for iPeriod=1:2%stim delay resp
    for iResult=1:2% Hit-Hit/FA-FA
        for iTrainStage=1:2
            CatResultCorrPerMouseSession=cat(1,CatResultCorrPerMouseSession,...
                ResultCorrPerMouseSessions{iTrainStage,iPeriod,iResult});
            label=cat(1,label,string(repmat([StageName{iTrainStage},'-',PeriodName{iPeriod},'-',ResultName{iResult}],...
                length(ResultCorrPerMouseSessions{iTrainStage,iPeriod,iResult}),1)));
        end
    end
end

close all

figure
hold on
bar(cellfun(@mean,ViolinOutput(:,1)),'w','BarWidth',.5,'LineWidth',2);
errorbar([1:8],cellfun(@mean,ViolinOutput(:,1)),sem(ViolinOutput(:,1),0),'Color',[0 0 0],'LineStyle','none','LineWidth',2,'CapSize',0)
for i=1:8
    scatter(ones(length(ViolinOutput{i,1}),1)*i,...
        ViolinOutput{i,1},'k','filled','jitter','on','jitterAmount',0.2);
end

temp=gcf;
temp.Parent.Children.CurrentAxes.XTick=[1:8];
temp.Parent.Children.CurrentAxes.XTickLabel=ViolinOutput(:,2);
xtickangle(45)
ylim([0 .5])
temp.Parent.Children.CurrentAxes.YTick=0:.1:.5;


clear TableResultMeanCorrSession
TableResultMeanCorrSession=table('Size',[0 5],'VariableTypes',{'double','string','string','string','string'},...
    'VariableNames',{'R','stage','period','design','ID'});
tempMouseID{1,1}=1:4;tempMouseID{1,2}=5:8;
ResultName={'Hit','FA'};

% change iPeriod iResult iTrainStage for subgroup statistical analysis
i=1;
for iPeriod=2%stim delay resp
    for iResult=2% Hit-Hit/FA-FA
        for iTrainStage=1:2
            for iBatch=1:2
                for iSub=1:size(DataAll.(BatchName{iBatch}),1)
                    for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                        % within session mean (2-3 sessions) as observation
                        temp=ResultMeanCorr{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(iResult,iResult);
                        if ~isnan(temp)
                            TableResultMeanCorrSession(i,:)={...
                                temp,StageName{iTrainStage},PeriodName{iPeriod},ResultName{iResult},...
                                string(tempMouseID{1,iBatch}(iSub))};
                            i=i+1;
                        end
                    end
                end
            end
        end
    end
end

LMM=fitlme(TableResultMeanCorrSession(:,[1 2 5]),'R~stage+(1|ID)')




%% pool correlation for Hit trial for Lick +/-
clear label* Cat* CorrMat
HitLickMeanCorr=[];HitNoLickMeanCorr=[];labelLick=[];labelNoLick=[];
TrialNumThr=6;
figure
i=1;
for iBatch=1:2 
    for iSub=1:size(DataAll.(BatchName{iBatch}),1)
        for iTrainStage=[1 2]%Exp Rev-Exp
            for iSession=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}))
                clear temp*
                temp=ZLowCa{1,iBatch}{iSub,iTrainStage}{iSession,1};
                tempMat1=[];tempMat2=[];tempMat3=[];
                for iTrial=1:length(DataAll.(BatchName{iBatch})(iSub).BHV.(StageName{iTrainStage}){iSession}.StimStart)
                    f0=mean(temp(StimBaseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:),1);
                    %stim
                    temp1=mean(temp(StimOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        StimOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:))-f0;
                    %delay
                    temp2=mean(temp(DelayOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        DelayOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:))-f0;
                    %resp
                    temp3=mean(temp(ResponseOnInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial):...
                        ResponseOffInd{1,iBatch}{iSub,iTrainStage}{iSession,1}(iTrial),:))-f0;
                    
                    tempMat1=cat(1,tempMat1,temp1);
                    tempMat2=cat(1,tempMat2,temp2);
                    tempMat3=cat(1,tempMat3,temp3);
                end
                
                % Hit
                temp1=[TrialResult{1,iBatch}{iSub,iTrainStage}{iSession,1} ...
                    IfLickDelayDLC{1,iBatch}{iSub,iTrainStage}{iSession,1}];
                tempLickIndHit=find(ismember(temp1,[1 1],'rows'));
                tempNoLickIndHit=find(ismember(temp1,[1 0],'rows'));
                % include session with Lick and Nolick >6 trials
                if (length(tempLickIndHit)>TrialNumThr)&(length(tempNoLickIndHit)>TrialNumThr)
%                     if true
                    tempMat1=zscore([tempMat1(tempLickIndHit,:);tempMat1(tempNoLickIndHit,:)]')';
                    tempMat2=zscore([tempMat2(tempLickIndHit,:);tempMat2(tempNoLickIndHit,:)]')';
                    tempMat3=zscore([tempMat3(tempLickIndHit,:);tempMat3(tempNoLickIndHit,:)]')';
                    
                    CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,1}=...
                        tempMat1*tempMat1'/(size(tempMat1,2)-1);%stim
                    CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,2}=...
                        tempMat2*tempMat2'/(size(tempMat2,2)-1);%delay
                    CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,3}=...
                        tempMat3*tempMat3'/(size(tempMat3,2)-1);%resp
                    
                    temp1=[1:length(tempLickIndHit)];
                    temp2=[length(tempLickIndHit)+1:(length(tempNoLickIndHit)+length(tempLickIndHit))];
                    
                    for iPeriod=2  %stim delay resp
                        % HitLick corr
                        HitLickMeanCorr=cat(1,HitLickMeanCorr,...
                            mean(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp1,temp1),'all'));
                        % HitNoLick corr
                        HitNoLickMeanCorr=cat(1,HitNoLickMeanCorr,...
                            mean(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod}(temp2,temp2),'all'));
                        labelLick=cat(1,labelLick,cellstr([StageName{iTrainStage}]));
                        labelNoLick=cat(1,labelNoLick,cellstr([StageName{iTrainStage}]));
                    end
                end
%                 subplot(8,15,(i-1)*15+(iTrainStage-1)*4+iSession)
%                 imagesc(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod})
%                 Lines(length(tempLickIndHit)+.5,[],'k','',1.5);
%                 Lines([],length(tempLickIndHit)+.5,'k','',1.5);
%                 axis square
%                 caxis([-.3 .6])
%                 xticks([1 size(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod},1)])
%                 yticks([1 size(CorrMat{1,iBatch}{iSub,iTrainStage}{iSession,iPeriod},1)])
%                 colormap(bluewhitered(256))
            end
        end
        i=i+1;
    end
end


g1=[labelLick;labelNoLick];
g2=[repmat(cellstr('Lick'),length(labelLick),1);repmat(cellstr('NoLick'),length(labelNoLick),1)];
anovan([HitLickMeanCorr;HitNoLickMeanCorr],{g1,g2})
figure
violinplot([HitLickMeanCorr;HitNoLickMeanCorr],append(g1,'-',g2));

tempCorr=[HitLickMeanCorr;HitNoLickMeanCorr];
[~,~,UniIdx]=unique(append(g1,'-',g2));
[ranksum(tempCorr(UniIdx==2),tempCorr(UniIdx==1)),ranksum(tempCorr(UniIdx==4),tempCorr(UniIdx==3))]

% Expert 2-nolick 1-lick
[~,~,stats]=ranksum(tempCorr(UniIdx==2),tempCorr(UniIdx==1))
% Naive 4-nolick 3-lick
[~,~,stats]=ranksum(tempCorr(UniIdx==4),tempCorr(UniIdx==3))

PlotValue={tempCorr(UniIdx==2),tempCorr(UniIdx==1),tempCorr(UniIdx==4),tempCorr(UniIdx==3)};

PlotLabel={'Exp-NoLick','Exp-Lick','Naive-NoLick','Naive-Lick'};

close all
figure
hold on
bar(cellfun(@mean,PlotValue),'w','BarWidth',.5,'LineWidth',2);
errorbar([1:4],cellfun(@mean,PlotValue),sem(PlotValue,0),'LineStyle','none','LineWidth',2,'CapSize',10)
for i=[1 3]
    scatter(ones(length(PlotValue{1,i}),1)*i,PlotValue{1,i},'k','filled','jitter','on','jitterAmount',0.2);
end
for i=[2 4]
    scatter(ones(length(PlotValue{1,i}),1)*i,PlotValue{1,i},'k','filled','LineWidth',2,'jitter','on','jitterAmount',0.2);
end
temp=gcf;
temp.Parent.Children.CurrentAxes.XTick=[1:4];
temp.Parent.Children.CurrentAxes.XTickLabel=PlotLabel;
ylim([0 .5])
xtickangle(45)