%% STATISTIC OXYTOCIN BEHAVIRO
dataname = ['F:\oxytocin\OXYT010\OXYT010_004_220126_112539',
            'F:\oxytocin\OXYT011\OXYT011_005_220224_160927',
            'F:\oxytocin\OXYT012\OXYT012_001_220323_130702',
            'F:\oxytocin\OXYT013\OXYT013_002_220412_124441',
            'F:\oxytocin\OXYT014\OXYT014_001_220429_155806',]
anaCh = [23;72;45;15;16]; OXT_REM_ALL = [];thetascore_all =[];REM_ALL = [];oxt_all_rem = [];
for ii = 1:5 %size(dataname,1)
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
%     load([basename,'.EMGFromLFP.LFP.mat']);
    load([basename,'.SleepScoreLFP.LFP.mat']);
    load([basename,'.SleepState.states.mat']);
    
    info = read_Intan_RHD2000_file(basepath);
    analogin = read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
    lfp_fbtime = find(analogin(1,:)>0.6,1);lfp_fbtime(1,2) = find(analogin(1,:)>0.6,1,'last');
    delete analogin
    % sitmulation and signal
    if exist('dFoverF.mat')==0
        name = basename(1:11);
        ext = '.mat'
        load([name,ext],[name,'_AChSti']);
        ACh_sti = eval([name,'_AChSti']);
        fbtime = ACh_sti.times(find(ACh_sti.values(:,1)>0.5,1));fbtime(1,2) = ACh_sti.times(find(ACh_sti.values(:,1)>0.5,1,'last'));
        fbIDX = round(fbtime*20000);
        load([name,ext],[name,'_AChSig1']);
        ACh_signal = eval([name,'_AChSig1']);
        ACh_signal.data = downsample(ACh_signal.values((fbIDX(1,1)):(fbIDX(1,2))),200);
        ACh_signal.time = downsample(ACh_signal.times((fbIDX(1,1)):(fbIDX(1,2))),200)- fbtime(1,1);
        p = [];
        p = polyfit(ACh_signal.time,ACh_signal.data,5);
        Ybaseline = polyval(p,ACh_signal.time);
        figure
        plot(ACh_signal.time, ACh_signal.data); hold on;
        plot(ACh_signal.time,Ybaseline);
        dFoverF = smoothdata((ACh_signal.data - Ybaseline)./Ybaseline,'movmean',10);
        hold on;
        plot(ACh_signal.time, dFoverF);xlim([0 inf]);
    else load('dFoverF.mat')
    end
    
    rippleCh = anaCh(ii)
    thetaCh = SleepScoreLFP.THchanID;
    THLFPs = bz_GetLFP(rippleCh,'basepath',basepath,'noPrompts',true);
    thetarange = linspace(0.5,250,1000);
    lfp = double(THLFPs.data);
    lfp(abs(lfp)>8000) = nan;
    lfp = fillmissing(lfp,'linear');
    [ss,ft,tt] = spectrogram(lfp,2*THLFPs.samplingRate,0.6*THLFPs.samplingRate,thetarange,THLFPs.samplingRate);
    [spectrogrammt,tm,fm] = MTSpectrogram(lfp,'window',2);
    bands = SpectrogramBands(spectrogrammt,fm);

    figure;
    colormap('parula')
    imagesc(tt,thetarange,abs(ss));clim([0,5e4]);
    yticks([0.5,20,50,100,200]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
    set(gca,'YDir','normal');
    hold on;plot([1:length(dFoverF)]/100+lfp_fbtime(1)/20000, dFoverF*500+100,'w');xlim([0 inf])
    
    %%% figure for REM
    if ii == 1
        figure;
        subplot(3,1,3)
        colormap('parula')
        imagesc(tt,thetarange,abs(ss));ylim([0 15]);clim([0 3e5]);xlim([9150 9500])
        yticks([0,5,10,20]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
        set(gca,'YDir','normal');
        subplot(3,1,2);plot(tm,bands.ratios.hippocampus );xlim([9150 9500])
        subplot(3,1,1);plot([1:length(dFoverF)]/100+lfp_fbtime(1)/20000, dFoverF);xlim([9150 9500])
    else
    end
    
    %% statistic
    OXT_REM = [];thetascore = [];
    REM = SleepState.ints.REMstate;peri = REM(:,2)-REM(:,1);
    for kk = 1:length(REM);
        win = REM(kk,1) + [-peri(kk) 2*peri(kk)]-round(lfp_fbtime(1)/20000);;
        OXT_REM(kk,1) = mean(dFoverF(win(1)*100:(win(1)*100+peri(kk)*100)));
        OXT_REM(kk,2) = mean(dFoverF((win(1)*100+peri(kk)*100):(win(1)*100+2*peri(kk)*100)));
        OXT_REM(kk,3) = mean(dFoverF((win(1)*100+2*peri(kk)*100):(win(1)*100+3*peri(kk)*100)));
        win = REM(kk,1) + [-peri(kk) 2*peri(kk)];
        thetascore(kk,1) = mean(bands.ratios.hippocampus(win(1):(win(1)+peri(kk))));
        thetascore(kk,2) = mean(bands.ratios.hippocampus((win(1)+peri(kk)):(win(1)+2*peri(kk))));
        thetascore(kk,3) = mean(bands.ratios.hippocampus((win(1)+2*peri(kk)):(win(1)+3*peri(kk))));
    end
    OXT_REM_ALL = [OXT_REM_ALL;OXT_REM];
    thetascore_all =[thetascore_all;thetascore];
    %% get rangking time
    oxt_rem = [];
    for kk = 1:length(REM);
        if REM(kk,2)+200<lfp_fbtime(2)/20000
            win = round(REM(kk,1)*100) + [-100 200]*100-round(lfp_fbtime(1)*100/20000);;
            oxt_rem(kk,:) = dFoverF(win(1):win(2));
        else
            REM(kk,:) = [];
        end
    end
    REM_ALL = [REM_ALL;REM];
    oxt_all_rem = [oxt_all_rem;oxt_rem];
end
%% all amplitude
figure;
x = 1:3;
e = std(OXT_REM_ALL)/ sqrt(length(OXT_REM_ALL));
yyaxis left
h = errorbar(x,mean(OXT_REM_ALL),e, 'o-');
set(h, 'MarkerSize', 10, 'MarkerFaceColor', [0.9290 0.6940 0.1250], ...
    'MarkerEdgeColor', [0.9290 0.6940 0.1250]);
xlim([0.5 3.5]);

x = 1:3;
e = std(thetascore_all)/sqrt(length(thetascore_all));
yyaxis right;k = errorbar(x,mean(thetascore_all),e, 'o-');
set(k, 'MarkerSize', 10, 'MarkerFaceColor', [0.47,0.67,0.19], ...
    'MarkerEdgeColor', [0.47,0.67,0.19]);

xlim([0.5 3.5]);ylim([0 5]);
%%
REM_duration = REM_ALL(:,2)-REM_ALL(:,1);
[a_sort, a_order] = sort(REM_duration);
oxt_rem_all = oxt_all_rem(a_order,:);
colormap('jet');
ts = 0;te = (a_sort+100)*100;
REM_ALL(:,3) =  a_order;

Ntrial = size(oxt_rem_all,1);
ave_oxt = mean(oxt_rem_all(:,1:10000),2);
miu = std(oxt_rem_all(:,1:10000),[],2);
oxt_rem_all_norm = (oxt_rem_all - ave_oxt)./miu;

figure;
imagesc(oxt_rem_all_norm)
zlim = [-8 3];
caxis(zlim)
colormap(NegativeEnhancingColormap(256, zlim, ...
        [0 0 1], [1 0 0],[1 1 1],1.25));
ylabel('Cell Number');colorbar
hold on
plot(te,1:Ntrial,'-.k','LineWidth',1);
hold on;
xline(10000,'-.k')
