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');
set(0, 'DefaultLineLineWidth', 1.5);

%%
clear all
close all
clc
fname = 'OXYT025_003_230101_131022';
load([fname,'OXT_CNO.mat']);
%%
fs = 20e3;
Tstart = OXT_CNO.lfp_fbtime(1)/fs;

fs_dff = 100;
Oxt_ripples_dFF = OXT_CNO.Oxt_dFF;
N_dFF = length(Oxt_ripples_dFF);
taxis = (0:1:N_dFF-1)/fs_dff;
ripple_times = OXT_CNO.ripples.peaks -Tstart;

%%
figure();
set(gcf,'Position',[100 100 600 250])
plot(taxis,100*Oxt_ripples_dFF,'r');
hold on
xline(ripple_times,'-.','Color',[0 0.4470 0.7410])
xlim([7400 7900])
box off
%%
dFoverF = zscore(Oxt_ripples_dFF);
dFF_event_mean = [];
count_ = 0;
time_window = 20;
for i = 1:(length(ripple_times))
    idx = round(ripple_times(i)*fs_dff);
    if (idx-time_window*fs_dff >0) && (idx+time_window*fs_dff<N_dFF)
        a = dFoverF(idx-time_window*fs_dff:idx+time_window*fs_dff);
        dFF_event_mean = [dFF_event_mean a];
    end
end

t_plot = -time_window:1/fs_dff:time_window;
figure();
subplot(2,1,1)
set(gcf,'Position',[100 100 400 700])
h = pcolor(t_plot,1:size(dFF_event_mean,2),dFF_event_mean');
set(h, 'EdgeColor', 'none');
shading interp;
set(gca,'YDir','reverse')
clim([-3,3])
colormap('parula')
%colormap(brewermap([],'-RdBu'));
colorbar
subplot(2,1,2)
plot(t_plot,(mean(dFF_event_mean,2)),'Color','r')
box off


%%
clear all
clc

fnames = {'OXYT025_004_230105_135258';
         'OXYT026_005_230113_165251';
         'OXYT026_008_230117_172559';
         'OXYT027_003_230104_141857'};
    %
fs = 20e3;
fs_dff = 100;
time_window = 20;
dFF_event_all = zeros([2*time_window*fs_dff+1,3,length(fnames) ]);
Ripple_event_all = zeros([81,3,length(fnames)]);
CNOTimes= [2 3]*3600;
for i = 1:length(fnames)
    load([fnames{i},'OXT_CNO.mat']);
    Tstart = OXT_CNO.lfp_fbtime(1)/fs;
    
    fs_dff = 100;
    Oxt_ripples_dFF = OXT_CNO.Oxt_dFF;
    N_dFF = length(Oxt_ripples_dFF);
    taxis = (0:1:N_dFF-1)/fs_dff;
    ripple_times = OXT_CNO.ripples.peaks -Tstart;

    
    LFPs = OXT_CNO.RipLFP.data(OXT_CNO.RipLFP.timestamps>Tstart);
    OXT_CNO.RipLFP.data = LFPs;
    times_ = OXT_CNO.RipLFP.timestamps( OXT_CNO.RipLFP.timestamps>Tstart)-Tstart;
    OXT_CNO.RipLFP.timestamps = times_;
    burtsWindow = 1;
    [dFF_event_mean1,ripplePower1] = GetRate(ripple_times,zscore(Oxt_ripples_dFF), OXT_CNO.RipLFP.data,burtsWindow,1,20,CNOTimes);
    [dFF_event_mean2,ripplePower2] = GetRate(ripple_times,zscore(Oxt_ripples_dFF), OXT_CNO.RipLFP.data,burtsWindow,2,20,CNOTimes);
    [dFF_event_mean3,ripplePower3] = GetRate(ripple_times,zscore(Oxt_ripples_dFF), OXT_CNO.RipLFP.data,burtsWindow,3,20,CNOTimes);
    
    dFF_event_all (:,1,i) = mean(dFF_event_mean1,2);
    dFF_event_all (:,2,i) = mean(dFF_event_mean2,2);
    dFF_event_all (:,3,i) = mean(dFF_event_mean3,2);
    Ripple_event_all (:,1,i) = ripplePower1;
    Ripple_event_all (:,2,i) = ripplePower2;
    Ripple_event_all (:,3,i) = ripplePower3;
    i
end

%%
t_plot = -time_window:1/fs_dff:time_window;
figure();
set(gcf,'Position',[100 100 400 600])
[Nt,Nc,Na] = size(dFF_event_all);
meanBaseLines = squeeze(mean(dFF_event_all(1:5*fs_dff,1,:),1));
dFF_plot = squeeze(dFF_event_all(:,1,:))-meanBaseLines';
options.handle     = figure(1);
subplot(2,1,1)
options.color_area = [0, 0.4,  0];   
options.color_line = [0, 0.4,  0, 0.3];
options.alpha      = 0.1;
options.line_width = 2;
options.error      = 'sem';
options.x_axis = t_plot;
plot_areaerrorbar(dFF_plot', options)
%plot(t_plot,mean(dFF_plot,2),'Color',[1, 0, 0, 0.4]);
hold on
meanBaseLines = squeeze(mean(dFF_event_all(1:5*fs_dff,2,:),1));
dFF_plot = squeeze(dFF_event_all(:,2,:))-meanBaseLines';
options.color_area = [0, 0.4,  0];   
options.color_line = [0, 0.4,  0, 0.6];
options.alpha      = 0.2;
plot_areaerrorbar(dFF_plot', options)
%plot(t_plot,mean(dFF_plot,2),'Color',[1, 0, 0, 0.7]);

meanBaseLines = squeeze(mean(dFF_event_all(1:5*fs_dff,3,:),1));
dFF_plot = squeeze(dFF_event_all(:,3,:))-meanBaseLines';
options.color_area = [0, 0.4,  0];   
options.color_line = [0, 0.4,  0, 1.0];
options.alpha      = 0.3;
plot_areaerrorbar(dFF_plot', options)
%plot(t_plot,mean(dFF_plot,2),'Color',[1, 0, 0, 1.0]);
box off
legend('Low Rate','Mid','High Rate')

tplot2 = linspace(-1*time_window,time_window, size(Ripple_event_all,1));
options.handle     = figure(1);
subplot(2,1,2)
options.color_area = [0 0 0];   
options.color_line = [0 0.0 0.0,0.3];
options.alpha      = 0.1;
options.line_width = 2;
options.error      = 'sem';
options.x_axis = tplot2;
plot_areaerrorbar(squeeze(Ripple_event_all(:,1,:))', options)
%plot(tplot2,squeeze(mean(Ripple_event_all(:,1,:),3)),'Color',[0 0.4470 0.7410,0.4]);
hold on
options.color_line = [0 0.0 0.0,0.6];
options.alpha      = 0.2;
plot_areaerrorbar(squeeze(Ripple_event_all(:,2,:))', options)
options.color_line = [0 0.0 0.0,1];
options.alpha      = 0.3;
plot_areaerrorbar(squeeze(Ripple_event_all(:,3,:))', options)
%plot(tplot2,squeeze(mean(Ripple_event_all(:,2,:),3)),'Color',[0 0.4470 0.7410,0.7]);
%plot(tplot2,squeeze(mean(Ripple_event_all(:,3,:),3)),'Color',[0 0.4470 0.7410,1.0]);
box off
legend('Low Rate','Mid','High Rate')
xlabel('Time (s)')
exportgraphics(gcf,'AllCNO_Oxt_ripples.pdf','BackgroundColor','none','ContentType','vector')

%%
%%
clear all
clc

fnames = {'OXYT025_004_230105_135258';
         'OXYT026_005_230113_165251';
         'OXYT026_008_230117_172559';
         'OXYT027_003_230104_141857'};
    %
fs = 20e3;
fs_dff = 100;

CNOTimes= [2 3]*3600;
RipPowerCNO = zeros([7200,length(fnames)]);
thetaPowerCNO = zeros([7200,length(fnames)]);

for i = 1:length(fnames)
    load([fnames{i},'OXT_CNO.mat']);
    Tstart = OXT_CNO.lfp_fbtime(1)/fs;
    ripple_times = OXT_CNO.ripples.peaks -Tstart;

    
    LFPs = OXT_CNO.RipLFP.data(OXT_CNO.RipLFP.timestamps>Tstart);
    [spectrogrammt,tm,fm] = MTSpectrogram(LFPs,'window',1);
    bands = SpectrogramBands(spectrogrammt,fm,'theta', [5,10]);
    tCNOIds = find((tm>CNOTimes(1)).*(tm<CNOTimes(2)));
    Ripp_power = (bands.ripples);
    RipPowerCNO(:,i) = Ripp_power(tCNOIds);
    thetaPower = (bands.theta);
    thetaPowerCNO(:,i) = thetaPower(tCNOIds);
    i
end

%%
figure
for i =1:4
    subplot(4,1,i)
histogram((RipPowerCNO(:,i)),200,'Normalization','pdf')
end
