
set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultTextFontSize', 10);
set(0,'defaultAxesFontSize',10);
set(0,'DefaultFigureColormap',jet)
set(0,'defaultfigurecolor',[1 1 1])
set(0,'DefaultAxesTitleFontWeight', 'normal')
set(0, 'defaulttextinterpreter', 'none') % or 'tex'

%% Load the data
clear all;clc
filename = dir('*_OXT_NoArti_ripples.mat');

for ii = 4 %:length(filename)-1
    %% load data
    basename = filename(ii).name;
    load(basename);
    %%
    fs = 1e4;
    Vscale = 1/(2.0^15)*10;
    OXT_EM = double(Sp2.OXT_Em)*Vscale;
 
    %%
    p = [];
    p = polyfit(Sp2.times,OXT_EM,2);
    
    Ybaseline = polyval(p,Sp2.times);
    figure
    plot(downsample(Sp2.times,10),downsample(OXT_EM - Ybaseline,10)); hold on;
    dFoverF = smoothdata((OXT_EM - Ybaseline)./Ybaseline,'movmean',10);
    hold on;
    plot(downsample(Sp2.times,10), downsample(dFoverF,10),'r');xlim([0 inf]);
    
    Ns = (floor(length(dFoverF)/100))*100;
    dFoverF_matrix = reshape(dFoverF(1:Ns), 100, []);
    dFoverF_ds = mean(dFoverF_matrix, 1);
    %% aritifical ripples dFoverF
    Npp = 1;
    NREM_flag = 3;
    WAKE_flag = 1;
    REM_flag = 5;
    pre = 5; % in seconds
    post = 5;  % in seconds
    smoothing_sigma = 15; % smoothing sigma in samples
    lfp_time = Sp2.lfp_time;
    
    fs_down = 100;
    windowSize = fs_down;
    %% Spontaneous ripples dFoverF

    ripple_times = Sp2.ripples.timestamps(:,1:2); 
    lfp_time = Sp2.lfp_time;
    ripples_synched= round(0.5*(ripple_times(:,1)+ripple_times(:,2)));
    ripples_synched = ripples_synched(ripples_synched>0);
    SleepState = Sp2.SleepState;
    ripples_sleepstate = SleepState.idx.states(ripples_synched);
    
    ripples_locs = round((Sp2.ripples.peaks-lfp_time(1))*100);
    stim_NREM  = ripples_locs(ripples_sleepstate==NREM_flag & circshift(ripples_sleepstate,5)==NREM_flag & circshift(ripples_sleepstate,-5)==NREM_flag);
    dFoverF_NREM = zeros(length(stim_NREM), pre*windowSize+post*windowSize+1);
    stim_NREM = stim_NREM((stim_NREM-pre*windowSize)>1);

    for i = 1:Npp:length(stim_NREM)-1
        eventTime = stim_NREM(i); % Get the current event time
        % Extract dFoverF values within the window around the event time
        windowStart = max(1, eventTime - pre*windowSize);
        windowEnd = min(length(dFoverF_ds), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds)
            dFoverF_NREM(i,:) = dFoverF_ds(windowStart:windowEnd) - mean(dFoverF_ds(windowStart+fs_down*(pre-3):windowStart+fs_down*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    figure
    set(gcf,'Position',[100 100 700 400])
    time_window = (-pre*windowSize:post*windowSize)/100;
    plot(time_window,zscore(smoothdata(mean(dFoverF_NREM),"gaussian",smoothing_sigma)),'b');
    hold
    
    stim_WAKE = ripples_locs(ripples_sleepstate==WAKE_flag & circshift(ripples_sleepstate,5)==WAKE_flag & circshift(ripples_sleepstate,-5)==WAKE_flag);
    dFoverF_WAKE = [];
    stim_WAKE = stim_WAKE((stim_WAKE-pre*windowSize)>1);
    for j = 1:Npp:length(stim_WAKE)
        eventTime = stim_WAKE(j); % Get the current event time
        % Extract dFoverF values within the window around the event time
        windowStart = max(1, eventTime - pre*windowSize);
        windowEnd = min(length(dFoverF_ds), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds)
            dFoverF_WAKE(j,:) = dFoverF_ds(windowStart:windowEnd) -mean(dFoverF_ds(windowStart+fs_down*(pre-3):windowStart+fs_down*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    time_window = (-pre*windowSize:post*windowSize)/100;
    plot(time_window,zscore(smoothdata(mean(dFoverF_WAKE),"gaussian",smoothing_sigma)),'k');
    
    stim_REM = ripples_locs(ripples_sleepstate==REM_flag & circshift(ripples_sleepstate,5)==REM_flag & circshift(ripples_sleepstate,-5)==REM_flag);
    dFoverF_REM = [];
    stim_REM = stim_REM((stim_REM-pre*windowSize)>1);
    for k = 1:Npp:length(stim_REM)
        eventTime = stim_REM(k); % Get the current event time
        % Extract dFoverF values within the window around the event time
        windowStart = max(1, eventTime - pre*windowSize);
        windowEnd = min(length(dFoverF_ds), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds)
            dFoverF_REM(k,:) = dFoverF_ds(windowStart:windowEnd) - mean(dFoverF_ds(windowStart+fs_down*(pre-3):windowStart+fs_down*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    time_window = (-pre*windowSize:post*windowSize)/100;
    plot(time_window,zscore(smoothdata(mean(dFoverF_REM),"gaussian",smoothing_sigma)),'m');
    ylabel('dF/F (zscored)');xlabel('Time (s)');xlim([-pre post]);ylim([-3 3]);
    box off
    legend('NREM','WAKE','REM',Location='best')
    legend('boxoff') 
    yline(0,'r--')
    xline(0,'r--')
    %exportgraphics(gcf,[basename(1:45),'dFoverF_Ripple_sleepstate.pdf'],'BackgroundColor','none','ContentType','vector');
    
    dFoverF_SleepState.NREM_Ripple = dFoverF_NREM;
    dFoverF_SleepState.WAKE_Ripple = dFoverF_WAKE;
    dFoverF_SleepState.REM_Ripple = dFoverF_REM;

    %%Save the data
    save([basename(1:45),'_dFoverF_Ripple_sleepstate.mat'],'dFoverF_SleepState'); 
     
    %% save fig
% Set figure size and position
figure('Position', [100, 100, 1200, 600]);
movmean_size = 30;
% Set plot options
options.color_area = [0 150 0]./255;    % Green theme
options.color_line = [0 150 0]./255;
options.alpha      = 0.1;
options.line_width = 2;
options.error      = 'sem';

% Define common x-axis
options.x_axis = [1:size(dFoverF_SleepState.NREM_Ripple,2)]./100-5;

% Plot NREM_ArtRipple_OXT
subplot(2,3,1);
plot_areaerrorbar(smoothdata(dFoverF_SleepState.NREM_Ripple,2,'movmean',movmean_size), options);
xlim([-5 5]); box off;
title('NREM_ArtRipple_OXT');
xline(0,'r--');
yline(0,'r--');

% Plot WAKE_ArtRipple_OXT
if isempty(dFoverF_SleepState.WAKE_Ripple)~=1
    subplot(2,3,2);
    plot_areaerrorbar(smoothdata(dFoverF_SleepState.WAKE_Ripple,2,'movmean',movmean_size), options);
    xlim([-5 5]); box off;
    title('WAKE-ArtRipple-OXT');
    xline(0,'r--');
end

% Plot REM_ArtRipple_OXT




% Export figure
exportgraphics(gcf,[basename(1:45),'_Control_Ripple_OXT.pdf'],'BackgroundColor','none','ContentType','vector');

end

