
set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultTextFontSize', 12);
set(0,'defaultAxesFontSize',12);
set(0,'DefaultFigureColormap',jet)
set(0,'defaultfigurecolor',[1 1 1])
set(0,'DefaultAxesTitleFontWeight', 'normal')
set(0, 'DefaultAxesFontName', 'Arial');

%%
clear all
close all
clc
fname = 'OXYT014_001_220429_155806';
load([fname,'Oxt_popRate_noRip.mat']');
load([fname, '.SleepState.states.mat']);
%
N_dFF = length(dFoverF);

fs = 20e3;
fs_dff = 100;
Tstart = floor(lfp_fbtime(1)/fs);
Nstart = floor(lfp_fbtime(1)/fs*fs_dff);
Tend = floor(lfp_fbtime(2)/fs);
Nend = floor(lfp_fbtime(2)/fs*fs_dff);
popRateAligned = double(popRate(Nstart:Nend));

%
sleep_ints = bz_INTtoIDX(SleepState.ints);
sleeps = sleep_ints.states;
sleeps(sleeps ==3) =5;
sleeps(sleeps ==2) =3;
sleepAligned = round(interp(sleeps(Tstart:Tend),fs_dff)*2)/2;
N_dff = length(dFoverF);
if N_dff< length(sleepAligned)
    sleepAligned = sleepAligned(1:N_dff);
end
%


histfit(popRateAligned,100,"rayleigh")
set(gca,'YScale','log')
box off
hold on
xline(mean(popRateAligned),'--k',"LineWidth",2)
colrs=colorGradient([0.75 0.75 0.75],[0,0,0],4);
xline(mean(popRateAligned)+0.5*std(popRateAligned),"Color",colrs(1,:),"LineWidth",2)
xline(mean(popRateAligned)+std(popRateAligned),"Color",colrs(2,:),"LineWidth",2)
xline(mean(popRateAligned)+2*std(popRateAligned),"Color",colrs(3,:),"LineWidth",2)
xline(mean(popRateAligned)+3*std(popRateAligned),"Color",colrs(4,:),"LineWidth",2)
xlabel('MUA (firing rate)')

%%

windowSize = 1;
thresholds =[1,2];
time_window = [15,20];
[dFF_eventNREM,dFF_eventWAKE,dFF_event_mean] = getFiringRateEvents(popRateAligned,(dFoverF),sleepAligned,thresholds,time_window);
figure
t_plot = -1*time_window(1):1/fs_dff:time_window(2);
plot(t_plot,(mean(dFF_event_mean,2)),'k')
hold on
plot(t_plot,mean(dFF_eventNREM,2),'b')
plot(t_plot,mean(dFF_eventWAKE,2),'m')

legend('ALL', 'NREM' ,'WAKE')


%%
clear all
clc
fnames = {'OXYT014_001_220429_155806';
         'OXYT011_005_220224_160927';
         'OXYT010_004_220126_112539'};
    %
    %
fs = 20e3;
fs_dff = 100;
time_window = [10,10];

fig1 = figure();
set(fig1,'Position',[100 100 400 1000])

ths =[0.5,1;1,2;2,3;3,10];
alphas_ = [0.25,0.5,0.75,1.0];
for kk = 1:size(ths,1)
dFF_NREM_all = zeros([sum(time_window*fs_dff)+1,length(fnames) ]);
dFF_REM_all = zeros([sum(time_window*fs_dff)+1,length(fnames) ]);
dFF_WAKE_all = zeros([sum(time_window*fs_dff)+1,length(fnames) ]);
for i = 1:length(fnames)
    load([fnames{i},'Oxt_popRate_noRip.mat']');
    load([fnames{i}, '.SleepState.states.mat']);
    N_dFF = length(dFoverF);
    
    fs = 20e3;
    fs_dff = 100;
    Tstart = floor(lfp_fbtime(1)/fs);
    Nstart = floor(lfp_fbtime(1)/fs*fs_dff);
    Tend = floor(lfp_fbtime(2)/fs);
    Nend = floor(lfp_fbtime(2)/fs*fs_dff);
    popRateAligned = double(popRate(Nstart:Nend));
    
    %
    sleep_ints = bz_INTtoIDX(SleepState.ints,'statename',{'WAKE','','NREM','','REM'});
    sleeps = sleep_ints.states;
    sleepAligned = round(interp(sleeps(Tstart:Tend),fs_dff)*2)/2;
    N_dff = length(dFoverF);
    if N_dff< length(sleepAligned)
        sleepAligned = sleepAligned(1:N_dff);
    end

    windowSize = 1;
    thresholds = ths(kk,:);
    [dFF_eventNREM,dFF_eventREM,dFF_eventWAKE,dFF_event_mean] = getFiringRateEvents(popRateAligned,(dFoverF),sleepAligned,thresholds,time_window);
    
    dFF_WAKE_all (:,i) = mean(dFF_eventWAKE,2);
    dFF_NREM_all (:,i) = mean(dFF_eventNREM,2);
    if ~isempty(dFF_eventREM)
        dFF_REM_all (:,i)  =  mean(dFF_eventREM,2);
    end
        i
end
t_plot = -1*time_window(1):1/fs_dff:time_window(2);

tbaseline = 5;
subplot(3,1,1)
dFF_event_all = dFF_NREM_all;
meanBaseLines = squeeze(mean(dFF_event_all(1:tbaseline*fs_dff,:),1));
dFF_plot_Oxt= dFF_event_all-meanBaseLines;
plot(t_plot,mean(dFF_plot_Oxt,2),'Color',[0,0.4, 0, alphas_(kk)]);ylim([-1.2 0.2])
hold on
title('NREM')
box off
ave_NREM(:,kk) = mean(dFF_plot_Oxt,2);

subplot(3,1,2)
dFF_event_all = dFF_WAKE_all;
meanBaseLines = squeeze(mean(dFF_event_all(1:tbaseline*fs_dff,:),1));
dFF_plot_Oxt= dFF_event_all-meanBaseLines;
plot(t_plot,mean(dFF_plot_Oxt,2),'Color',[0,0.4, 0, alphas_(kk)]);ylim([-1.2 0.2])
hold on
title('WAKE')
box off
ave_WAKE(:,kk) = mean(dFF_plot_Oxt,2);


subplot(3,1,3)
dFF_event_all = dFF_REM_all ;
meanBaseLines = squeeze(mean(dFF_event_all(1:tbaseline*fs_dff,:),1));
dFF_plot_Oxt= dFF_event_all-meanBaseLines;
plot(t_plot,mean(dFF_plot_Oxt,2),'Color',[1.0,0, 0, alphas_(kk)]);
hold on
title('REM')
box off
ave_NREM(:,kk) = mean(dFF_plot_Oxt,2);


end
legend('0.5-1', '1-2', '2-3','3-max',Location='best')
exportgraphics(gcf,'All_OXT_MUAs.pdf','BackgroundColor','none','ContentType','vector')
%% check the statistic
data = [data_group1(:); data_group2(:); data_group3(:)];
groups = [ones(size(data_group1(:))); 2*ones(size(data_group2(:))); 3*ones(size(data_group3(:)))];

[p,tbl,stats] = anova1(data, groups);

fvalue = tbl{2,5};
pvalue = tbl{2,6};

comparison = multcompare(stats);

fprintf('F-value: %f, P-value: %f\n', fvalue, pvalue);
%%
    figure
    thres = [3 2.5 2 1.5 1];
    alpha_val = linspace(1.0,0.2, length(thres));
    for ii = 1:length(thres)
        durations = 10;
        if ii>1 && ii<6
            [pks locs_a width] = findpeaks(popRate,'minpeakheight',mean(popRate)+thres(ii)*std(popRate),...
                'WidthReference','halfprom','MinPeakDistance',durations(1)*2);        
            [pks locs_b width] = findpeaks(popRate,'minpeakheight',mean(popRate)+thres(ii-1)*std(popRate),...
                'WidthReference','halfprom','MinPeakDistance',durations(1)*2);
            locs = setdiff(locs_a,locs_b);            
        else
            [pks locs width] = findpeaks(popRate,'minpeakheight',mean(popRate)+thres(ii)*std(popRate),...
            'WidthReference','halfprom','MinPeakDistance',durations(1)*2);
        end
        
        win = [-25 30];
        events = round(locs+[win(1) win(2)]*100-Tstart*100);
        dF_pop = [];
        events = events(events(:,1)>0,:);
        events(:,2) = events(:,1) + diff(win)*100;
        events = events(events(:,2)<length(dFoverF),:); 
        for jj = 1:length(events)
           dF_pop(:,jj) = dFoverF(events(jj,1):events(jj,2));
        end
        miu = mean(dF_pop(1:1000,:));
        sigma = std(dF_pop(1:1000,:));
        dF_pop_norm = (dF_pop - miu)./sigma;
        hold on;
        cc(ii) = plot(mean(dF_pop_norm,2), 'k');
        cc(ii).Color = [cc(ii).Color alpha_val(ii)];
        title('PYR-ACh');xlim([1 inf]);
    end
