%% STATISTIC OXYTOCIN BEHAVIRO
dataname = ['G:\data\LAMP2\LAMP2_002_190405_154550']
            %'G:\data\LAMP3\LAMP3_006_190613_194519',]
dataname = [ %'G:\data\LAMP37\LAMP37_005_201111_204459'
            'G:\data\LAMP37\LAMP37_006_201201_192319']
dataname = [ %'F:\LIMP009\LIMP009_002_210714_091509',
            'F:\LIMP011\LIMP011_010_210705_120340',
            'F:\LIMP013\LIMP013_009_210630_101849',
            'F:\LIMP014\LIMP014_012_211005_112625',
            'F:\LIMP016\LIMP016_004_220101_144210',
            'F:\LIMP019\LIMP019_004_220303_115500']
ACh_REM_ALL = [];thetascore_all =[];REM_ALL = [];ach_all_rem = [];

for ii =1: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);
    if ii<4
        digitalin = read_Intan_digitalin_file('digitalin.dat',info.num_board_dig_in_channels);
        lfp_fbtime = find(digitalin> (min(digitalin)+max(digitalin))/2,1);lfp_fbtime(1,2) = find(digitalin> (min(digitalin)+max(digitalin))/2,1,'last');
        delete digitalin
    else 
        analogin = read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
        lfp_fbtime = find(analogin(2,:)> (min(analogin(1,:))+max(analogin(2,:)))/2,1);lfp_fbtime(1,2) = find(analogin(2,:)> (min(analogin(1,:))+max(analogin(2,:)))/2,1,'last');
        delete analogin
    end
    % sitmulation and signal
    if exist('dFoverF.mat')==0
        name = basename(1:11);
        ext = '.mat'
        if ii<4
            load([name,ext],[name,'_OpTrTTL']);
            ACh_sti = eval([name,'_OpTrTTL']);
            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'));
        else 
            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'));
        end
        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')
        figure; plot(dFoverF);xlim([0 inf]);
    end
    
    rippleCh = SleepScoreLFP.THchanID
    thetaCh = SleepScoreLFP.THchanID;
    THLFPs = bz_GetLFP(thetaCh,'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,5e5]);
    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])
    
    %% statistic
    ACh_REM = [];thetascore = [];
    REM_t = []; % need to manual input from excel
    peri = REM_t(:,2)-REM_t(:,1);
    for kk = 1:length(REM_t);
        if REM_t(kk,1)-100>0 & REM_t(kk,2)+ 200 < lfp_fbtime(2)/20000
            win = REM_t(kk,1) + [-peri(kk) 2*peri(kk)];;
            ACh_REM(kk,1) = mean(dFoverF(win(1)*100:(win(1)*100+peri(kk)*100)));
            ACh_REM(kk,2) = mean(dFoverF((win(1)*100+peri(kk)*100):(win(1)*100+2*peri(kk)*100)));
            ACh_REM(kk,3) = mean(dFoverF((win(1)*100+2*peri(kk)*100):(win(1)*100+3*peri(kk)*100)));
            win = REM_t(kk,1) + [-peri(kk) 2*peri(kk)]+round(lfp_fbtime(1)/20000);
            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
    end
    ACh_REM_ALL = [ACh_REM_ALL;ACh_REM];
    thetascore_all =[thetascore_all;thetascore];
    figure;plot(bands.ratios.hippocampus);xlim([0 inf])
    %% get rangking time
    ach_rem = [];
    for kk = 1:length(REM_t);
        if REM_t(kk,1)-100>0 & REM_t(kk,2)+ 200 < lfp_fbtime(2)/20000
            win = round(REM_t(kk,1)*100) + [-100 200]*100;;
            ach_rem(kk,:) = dFoverF(win(1):win(2));
            REM_ALL = [REM_ALL;REM_t(kk,:)];
        end
    end
    ach_all_rem = [ach_all_rem;ach_rem];
end
%% all amplitude
figure;
x = 1:3;
e = std(ACh_REM_ALL)/ sqrt(length(ACh_REM_ALL));
yyaxis left
h = errorbar(x,mean(ACh_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 = ach_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')
