
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

animalNames = {'OXYT010_004';
    'OXYT011_005';
    'OXYT014_001';
    'OXYT012_001'};

%
fs1 = 100;
animalID = [];
sessionID =[];
filenamesID =[];
sessionCount = 0;
rippleDurs = {};
for anID = 1:length(animalNames)
    Names_LIMP = dir([animalNames{anID},'*Oxt_LFP.mat']);
    filenamesLIMP= {Names_LIMP.name};
    foldernameLIMP = {Names_LIMP.folder};
    load(fullfile(foldernameLIMP{1},filenamesLIMP{1}));
    Names_LIMPSleep = dir([animalNames{anID},'*SleepState.states.mat']);
    filenamesLIMPSleep = {Names_LIMPSleep.name};
    load(fullfile(foldernameLIMP{1},filenamesLIMPSleep{1}));
    %
    fs_ripple = RipLFPs.samplingRate;
    
    N1 = ceil(lfp_fbtime(1)/2e4*fs_ripple);
    N2 = floor(lfp_fbtime(2)/2e4*fs_ripple);
    ripple_ = RipLFPs.data(N1:N2);
    Theta_ = THLFPs.data(N1:N2);

    % Find Throughs
    fs1 = 100;
    tstart = ceil(lfp_fbtime(1)/2e4);
    tend = floor(lfp_fbtime(2)/2e4);
    sleep_ints = bz_INTtoIDX(SleepState.ints,statenames ={'WAKE','','NREM','','REM'});
    sleeps = sleep_ints.states;
    sleep_ = round(round(interp(sleeps(tstart:tend),fs1)*5)/5);
    N_dff = length(dFoverF);
    if N_dff< length(sleep_)
        sleep_ = sleep_(1:N_dff);
    end
    dFoverF_ = zeros(size(dFoverF));
    dFoverF_(sleep_ == 1) = dFoverF(sleep_ == 1);%
    peak_mean_diff = 20;
    threshold_ = 0.03;
    [pks, locs, width] = findpeaks(dFoverF_,'MinPeakProminence',threshold_,...
        'WidthReference','halfprom','MinPeakDistance',peak_mean_diff*fs1);
  
    dF_event_mean = [];
    dFoverF_ = zeros(size(dFoverF));
    dFoverF_(sleep_ == 1) = dFoverF(sleep_ == 1);
    dFoverF_(dFoverF_ ==0) = NaN;
    locs_min =[];
    count_ = 0;
    for i = 1:length(locs)
        if (locs(i)-peak_mean_diff*fs1 >0) && (locs(i)+peak_mean_diff*fs1<N_dff)
            a = dFoverF_(locs(i)-peak_mean_diff*fs1:locs(i)+peak_mean_diff*fs1);
            if ~anynan(a)
                count_= count_ + 1;
                dF_event_mean = [dF_event_mean a];
                locs_min = [locs_min locs(i)];
            end
        end
    end
    %
    options.handle     =  figure();
    set(gcf,'Position',[100 100 500 250])
    options.color_area = [255 0 0]./255;
    options.color_line = [255 0 0]./255;
    options.alpha      = 0.5;
    options.line_width = 2;
    options.error      = 'sem';
    options.x_axis = -peak_mean_diff:1/fs1:peak_mean_diff;
    plot_areaerrorbar(dF_event_mean',options)
    box off
    hold on


    % locs_insec = (locs_min/fs1 + lfp_fbtime(1)/2e4)';
    % size(locs_insec)
    % [ wavAvg, lfpAvg ] = bz_eventWavelet (THLFPs, locs_insec,'twin', [20,20],'plotWave',true);
    % 
    % fnames = split(filenamesLIMPSleep,'.');
    % save([fnames{1},'WAKE_Peak_wavelet_lfpChan.mat'],'wavAvg','lfpAvg')

end

%% WAKE Gamma
figure;
set(gcf,'Position',[100 100 500 300])
h = pcolor(downsample(wavAvg.timestamps,50),wavAvg.freqs,downsample(wavAvg.data,50)');hold on;
set(h, 'EdgeColor', 'none');
shading interp;
set(gca,'YDir','normal')
set(gca,'YScale','log');ylim([30 250]);
% ylim([wavespec.freqs(1) wavespec.freqs(end)]);
colormap('jet')
xlabel('time (ms)'); ylabel('frequency (Hz)');
colorbar
clim([70, 270])
%xline(0,'m','LineWidth',2)
yyaxis right
df = diff(wavAvg.freqs);
wavAvg2 = 0.5*(wavAvg.data(:,2:end) + wavAvg.data(:,1:end-1));
wavAvg2_band = wavAvg2.*df;
Fband = wavAvg.freqs(76)-wavAvg.freqs(87);
wavAvg_smoothed = smoothdata(mean(wavAvg2_band(:,76:87),2),'sgolay',2e3);
plot(wavAvg.timestamps,wavAvg_smoothed,'k')
fs_ripple = 1250;
a = (-mean(wavAvg_smoothed(1:5*fs_ripple)) +  mean(wavAvg_smoothed(15*fs_ripple:20*fs_ripple)))/mean(wavAvg_smoothed(1:5*fs_ripple))
%set(gcf,'renderer','Painters')


%%

animalNames = {'OXYT010_004';
    'OXYT011_005';
    'OXYT014_001';
    'OXYT012_001'};

for anID = 1:length(animalNames)
    Names_LIMP = dir([animalNames{anID},'*NREM_Peak*.mat']);
    filenamesLIMP= {Names_LIMP.name};
    foldernameLIMP = {Names_LIMP.folder};
    load(fullfile(foldernameLIMP{1},filenamesLIMP{1}));
    figure;
    set(gcf,'Position',[100 100 500 300])
    h = pcolor(downsample(wavAvg.timestamps,50),wavAvg.freqs,downsample(wavAvg.data,50)');hold on;
    set(h, 'EdgeColor', 'none');
    shading interp;
    set(gca,'YDir','normal')
    set(gca,'YScale','log');
    colormap('jet')
    xlabel('time (ms)'); ylabel('frequency (Hz)');
    colorbar
    clim([40, 120])
    yyaxis right
    df = diff(wavAvg.freqs);
    wavAvg2 = 0.5*(wavAvg.data(:,2:end) + wavAvg.data(:,1:end-1));
    wavAvg2_band = wavAvg2.*df;
    Fband = wavAvg.freqs(99)-wavAvg.freqs(87);
    wavAvg_smoothed = smoothdata(mean(wavAvg2_band(:,87:99),2),'sgolay',1e3);
    plot(wavAvg.timestamps,wavAvg_smoothed,'k')
    
    fs_ripple = 1250;
    PNREM_Peak(anID) = 100*(-mean(wavAvg_smoothed(10*fs_ripple:15*fs_ripple)) +  mean(wavAvg_smoothed(17*fs_ripple:21*fs_ripple)))/mean(wavAvg_smoothed(1:5*fs_ripple))


end
%
for anID = 1:length(animalNames)
    Names_LIMP = dir([animalNames{anID},'*NREM_Trough*.mat']);
    filenamesLIMP= {Names_LIMP.name};
    foldernameLIMP = {Names_LIMP.folder};
    load(fullfile(foldernameLIMP{1},filenamesLIMP{1}));
    figure;
    set(gcf,'Position',[100 100 500 300])
    h = pcolor(downsample(wavAvg.timestamps,50),wavAvg.freqs,downsample(wavAvg.data,50)');hold on;
    set(h, 'EdgeColor', 'none');
    shading interp;
    set(gca,'YDir','normal')
    set(gca,'YScale','log');
    colormap('jet')
    xlabel('time (ms)'); ylabel('frequency (Hz)');
    colorbar
    clim([40, 120])
    yyaxis right
    df = diff(wavAvg.freqs);
    wavAvg2 = 0.5*(wavAvg.data(:,2:end) + wavAvg.data(:,1:end-1));
    wavAvg2_band = wavAvg2.*df;
    Fband = wavAvg.freqs(99)-wavAvg.freqs(87);
    wavAvg_smoothed = smoothdata(mean(wavAvg2_band(:,87:99),2),'sgolay',2e3);
    plot(wavAvg.timestamps,wavAvg_smoothed,'k')
    
    fs_ripple = 1250;
    PNREM_Trough(anID) = 100*(-mean(wavAvg_smoothed(1*fs_ripple:5*fs_ripple)) +  mean(wavAvg_smoothed(17*fs_ripple:21*fs_ripple)))/mean(wavAvg_smoothed(1:5*fs_ripple))
end

%
for anID = 1:length(animalNames)
    Names_LIMP = dir([animalNames{anID},'*WAKE_Trough*.mat']);
    filenamesLIMP= {Names_LIMP.name};
    foldernameLIMP = {Names_LIMP.folder};
    load(fullfile(foldernameLIMP{1},filenamesLIMP{1}));
    figure;
    set(gcf,'Position',[100 100 500 300])
    h = pcolor(downsample(wavAvg.timestamps,50),wavAvg.freqs,downsample(wavAvg.data,50)');hold on;
    set(h, 'EdgeColor', 'none');
    shading interp;
    set(gca,'YDir','normal')
    set(gca,'YScale','log');
    colormap('jet')
    xlabel('time (ms)'); ylabel('frequency (Hz)');
    colorbar
    clim([40, 120])
    yyaxis right
    df = diff(wavAvg.freqs);
    wavAvg2 = 0.5*(wavAvg.data(:,2:end) + wavAvg.data(:,1:end-1));
    wavAvg2_band = wavAvg2.*df;
    Fband = wavAvg.freqs(76)-wavAvg.freqs(87);
    wavAvg_smoothed = smoothdata(mean(wavAvg2_band(:,76:87),2),'sgolay',2e3);
    plot(wavAvg.timestamps,wavAvg_smoothed,'k')
    
    fs_ripple = 1250;
    PWake_Trough(anID) = 100*(-mean(wavAvg_smoothed(10*fs_ripple:15*fs_ripple)) +  mean(wavAvg_smoothed(17*fs_ripple:21*fs_ripple)))/mean(wavAvg_smoothed(1:5*fs_ripple))
end


%
for anID = 1:length(animalNames)
    Names_LIMP = dir([animalNames{anID},'*WAKE_Peak*.mat']);
    filenamesLIMP= {Names_LIMP.name};
    foldernameLIMP = {Names_LIMP.folder};
    load(fullfile(foldernameLIMP{1},filenamesLIMP{1}));
    figure;
    set(gcf,'Position',[100 100 500 300])
    h = pcolor(downsample(wavAvg.timestamps,50),wavAvg.freqs,downsample(wavAvg.data,50)');hold on;
    set(h, 'EdgeColor', 'none');
    shading interp;
    set(gca,'YDir','normal')
    set(gca,'YScale','log');
    colormap('jet')
    xlabel('time (ms)'); ylabel('frequency (Hz)');
    colorbar
    yyaxis right
    df = diff(wavAvg.freqs);
    wavAvg2 = 0.5*(wavAvg.data(:,2:end) + wavAvg.data(:,1:end-1));
    wavAvg2_band = wavAvg2.*df;
    Fband = wavAvg.freqs(76)-wavAvg.freqs(86);
    wavAvg_smoothed = smoothdata(mean(wavAvg2_band(:,76:87),2),'sgolay',2e3);
    plot(wavAvg.timestamps,wavAvg_smoothed,'k')
    
    fs_ripple = 1250;
    PWake_Peak(anID) = 100*(-mean(wavAvg_smoothed(10*fs_ripple:15*fs_ripple)) +  mean(wavAvg_smoothed(17*fs_ripple:21*fs_ripple)))/mean(wavAvg_smoothed(1:5*fs_ripple))
end



%% 
AnimalNames = animalNames;
Wake_Peak = PWake_Peak';
Wake_Trough = PWake_Trough';

NREM_Peak = PNREM_Peak';
NREM_Trough = PNREM_Trough';

Oxt_Results = table(NREM_Trough,NREM_Peak,Wake_Trough,Wake_Peak);

%%
colr = [0,0.5,0];
scatter([0.5,1.5,2.5,3.5],Oxt_Results{3,:},70,colr,'o',LineWidth=2)
hold on
scatter([0.5,1.5,2.5,3.5],Oxt_Results{4,:},70,colr,'d',LineWidth=2)
scatter([0.5,1.5,2.5,3.5],Oxt_Results{2,:},70,colr,'s',LineWidth=2)
scatter([0.5,1.5,2.5,3.5],Oxt_Results{1,:},70,colr,'p',LineWidth=2)

xlim([-0.5,4])
yline([0],'--')
