
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_Arti_ripples.mat');

for ii = 10  %2:length(filename)
    %% load data
    basename = filename(ii).name;
    load(basename);
    %%
    fs = Sp2.FluoData.sampling_rate;
    LaserShutter = double(Sp2.FluoData.laser_times); %%camsignal is the laser signal, digital_2 is cam signal
    
    %%
    OXT_dff = double(Sp2.FluoData.dFF_OXT);
    figure;plot(Sp2.FluoData.time_seconds,OXT_dff+0.1);hold on
    peakLocation = round(0.5*(LaserShutter(:,1)+LaserShutter(:,2))*fs);
    samplesToRemove = 16;
    for i = 1:length(peakLocation)
        % Calculate the range of samples to remove
        rangeStart = max(1, peakLocation(i) - samplesToRemove);
        rangeEnd = min(length(OXT_dff), peakLocation(i) + round(1.5*samplesToRemove));
        
        % Remove the samples within the defined range
        OXT_dff(rangeStart:rangeEnd) = NaN;
    end
    OXT_dff = fillmissing(OXT_dff, 'nearest');
    
    plot(Sp2.FluoData.time_seconds,OXT_dff);
    
    %%
    
    pre = 5;
    post = 5;
    eventAveragedDF = []; % Initialize array for event-averaged dFoverF
    Ns = length(OXT_dff);
    windowSize = fs; % Define the window size around each event time
    dFoverFWindow = [];
    stim_locs = round(0.5*(LaserShutter(:,1)+LaserShutter(:,2))*130);
    for i = 1:length(stim_locs)
        eventTime = stim_locs(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(OXT_dff), eventTime + post*windowSize);
        if windowEnd<length(OXT_dff)
            dFoverFWindow(i,:) = OXT_dff(windowStart:windowEnd) - mean(OXT_dff(windowStart+fs*3:windowStart+fs*5));
        end
        % Calculate the average of the extracted dFoverF values
    end
    
    % plot the event-averaged dFoverF
    figure
    time_window = (-pre*windowSize:post*windowSize)/fs;
    plot(time_window,smoothdata(mean(dFoverFWindow),"gaussian",20),'r');
    
%% alighed dff with cam signal
    cam_signal = double(Sp2.FluoData.cam_signal); 
    % Define your threshold
    threshold = 0.5;                                                        
    % Calculate the difference between consecutive elements
    diff_signal = diff(cam_signal > threshold);
    % Find the start times (where diff_signal is 1)
    start_times = find(diff_signal == 1,1) + 1;
    % Find the end times (where diff_signal is -1)
    end_times = find(diff_signal == -1,1,'last');
    dFoverF_ds_aligned = OXT_dff(start_times:end_times);
    LaserShutter_aligned = LaserShutter - Sp2.FluoData.time_seconds(start_times);
    stim_locs_aligned = stim_locs - start_times;
%% 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;
    pulsetimes_synched= round(0.5*(LaserShutter_aligned(:,1)+LaserShutter_aligned(:,2)) + lfp_time(1));
    pulsetimes_synched = pulsetimes_synched(pulsetimes_synched>0);
    SleepState = Sp2.SleepState;
    pulse_sleepstate = SleepState.idx.states(pulsetimes_synched);
    
    stim_NREM  = stim_locs_aligned(pulse_sleepstate==NREM_flag & circshift(pulse_sleepstate,5)==NREM_flag & circshift(pulse_sleepstate,-5)==NREM_flag);
    dFoverF_NREM = zeros(length(stim_NREM), pre*windowSize+post*windowSize+1);
    %stim_locs = round(0.5*(pulse_times_aligned(:,1)+pulse_times_aligned(:,2))*130);
    
    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_aligned), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds_aligned)
            dFoverF_NREM(i,:) = dFoverF_ds_aligned(windowStart:windowEnd) - mean(dFoverF_ds_aligned(windowStart+fs*(pre-3):windowStart+fs*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)/130;
    plot(time_window,zscore(smoothdata(mean(dFoverF_NREM),"gaussian",smoothing_sigma)),'b');
    hold
    
    stim_WAKE = stim_locs_aligned(pulse_sleepstate==WAKE_flag & circshift(pulse_sleepstate,5)==WAKE_flag & circshift(pulse_sleepstate,-5)==WAKE_flag);
    dFoverF_WAKE = [];
    %stim_locs = round(0.5*(pulse_times(:,1)+pulse_times(:,2))*100);
    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_aligned), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds_aligned)
            dFoverF_WAKE(j,:) = dFoverF_ds_aligned(windowStart:windowEnd) -mean(dFoverF_ds_aligned(windowStart+fs*(pre-3):windowStart+fs*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    time_window = (-pre*windowSize:post*windowSize)/130;
    plot(time_window,zscore(smoothdata(mean(dFoverF_WAKE),"gaussian",smoothing_sigma)),'k');
    
    stim_REM = stim_locs_aligned(pulse_sleepstate==REM_flag & circshift(pulse_sleepstate,5)==REM_flag & circshift(pulse_sleepstate,-5)==REM_flag);
    dFoverF_REM = [];
    %stim_locs = round(0.5*(pulse_times(:,1)+pulse_times(:,2))*100);
    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_aligned), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds_aligned)
            dFoverF_REM(k,:) = dFoverF_ds_aligned(windowStart:windowEnd) - mean(dFoverF_ds_aligned(windowStart+fs*(pre-3):windowStart+fs*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    time_window = (-pre*windowSize:post*windowSize)/130;
    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_ArtRipple_sleepstate.pdf'],'BackgroundColor','none','ContentType','vector');

    dFoverF_SleepState.NREM_ArtRipple = dFoverF_NREM;
    dFoverF_SleepState.WAKE_ArtRipple = dFoverF_WAKE;
    dFoverF_SleepState.REM_ArtRipple = dFoverF_REM;
    %% 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))*fs);
    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_aligned), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds_aligned)
            dFoverF_NREM(i,:) = dFoverF_ds_aligned(windowStart:windowEnd) - mean(dFoverF_ds_aligned(windowStart+fs*(pre-3):windowStart+fs*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)/fs;
    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_aligned), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds_aligned)
            dFoverF_WAKE(j,:) = dFoverF_ds_aligned(windowStart:windowEnd) -mean(dFoverF_ds_aligned(windowStart+fs*(pre-3):windowStart+fs*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    time_window = (-pre*windowSize:post*windowSize)/fs;
    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_aligned), eventTime + post*windowSize);
        if windowEnd<length(dFoverF_ds_aligned)
            dFoverF_REM(k,:) = dFoverF_ds_aligned(windowStart:windowEnd) - mean(dFoverF_ds_aligned(windowStart+fs*(pre-3):windowStart+fs*pre));
        end
        % Calculate the average of the extracted dFoverF values
    end
    time_window = (-pre*windowSize:post*windowSize)/fs;
    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 = 15;
% Set plot options
options.handle = figure;
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_ArtRipple,2)]./fs-5;

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

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

% Plot REM_ArtRipple_OXT
subplot(2,3,3);
plot_areaerrorbar(smoothdata(dFoverF_SleepState.REM_ArtRipple,2,'movmean',movmean_size), options);
xlim([-5 5]); box off;
title('REM-ArtRipple-OXT');
xline(0,'r--');
yline(0,'r--');

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

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

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

end



