set(0, 'DefaultLineLineWidth', 3);
set(0,'defaultTextFontSize', 10);
set(0,'defaultAxesFontSize',10);
set(0,'DefaultFigureColormap',jet)
set(0,'defaultfigurecolor',[1 1 1])
set(0,'DefaultAxesTitleFontWeight', 'normal')
%%
clear all;clc

animalNames = {'ACh_OXT004_012';
    'ACh_OXT005_010';
    'ACh_OXT006_006';
    'ACh_OXT007_009'};

%
fs_fluo = 100;
for anID = 1:length(animalNames)
    Names_AchOXT = dir([animalNames{anID},'*Act.mat']);
    filenamesAchOXT= {Names_AchOXT.name};
    foldernameAchOXT = {Names_AchOXT.folder};
    load(fullfile(foldernameAchOXT{1},filenamesAchOXT{1}));

    Ndff = length(ACh_OXT_Act.Ach_dFF);
    tt = 0:1/fs_fluo:(Ndff/fs_fluo)-1/fs_fluo;
    FluoData.time = tt';
    FluoData.Ach_dFF= ACh_OXT_Act.Ach_dFF;
    FluoData.Oxt_dFF= ACh_OXT_Act.Oxt_dFF;

    Fluo_Categorized = Categorize_Responses(FluoData,ACh_OXT_Act.PulseData);
    % Test Plot
    fs_fluo = Fluo_Categorized.fs;
    Fluo = Fluo_Categorized.dFF;
    PulseData = ACh_OXT_Act.PulseData;
    plotoptions.handle  = figure();
    plotoptions.error = 'sem';
    plotoptions.line_width = 2;
    plotoptions.alpha      = 0.2;

    Stats = {};
    Ach_max = [];
    dFF_pulse_stats ={};
    indpulse=0;
    for i = 1:length(fieldnames(Fluo.Mixed))
        
        amp = mean(PulseData.pulse_amplitude(PulseData.pulse_id==i));
        dur = round(mean(PulseData.pulse_dur(PulseData.pulse_id==i)));
        eval(['xx = -10:1/fs_fluo:length(Fluo.Mixed.pulse' num2str(i) '.Ach(:,1))/fs_fluo-1/fs_fluo-10;']);
    
        plotoptions.x_axis = xx;
        eval(['y1 = Fluo.Mixed.pulse' num2str(i) '.Ach;']);
        eval(['y2 = Fluo.Mixed.pulse' num2str(i) '.Oxt;']);
        if size(y1,2)>10
            indpulse  = (round(amp/5) -1)*3 +round(dur/5) + 1;
            dFF_pulse_stats.amp(indpulse) = amp;
            dFF_pulse_stats.dur(indpulse) = dur;
            
            subplot(4,3,indpulse)
            y1mean = mean(y1,2);
            dFF_pulse_stats.Ach_max(indpulse) =  peak2peak(y1mean(9*fs_fluo:(10+dur)*fs_fluo));
            dFF_pulse_stats.Ach{indpulse} = y1mean;
    
     
            plotoptions.color_area = [255 0 0]./255;
            plotoptions.color_line = [ 255 0 0]./255;
            plot_areaerrorbar(y1',plotoptions)
            hold on
            plotoptions.color_area = [0 140 0]./255;
            plotoptions.color_line = [0 140 0]./255;
            
            y2mean = mean(y2,2);
            dFF_pulse_stats.oxt_min(indpulse) =  peak2peak(y2mean(9*fs_fluo:(10+dur+1)*fs_fluo));
            dFF_pulse_stats.oxt{indpulse} = y2mean;
            [y2_min,inds] = min(y2mean(1001:1200),[],1);
            dFF_pulse_stats.oxtmin(indpulse) = y2_min;
            dFF_pulse_stats.min_time(indpulse) = inds/fs_fluo;
            dFF_pulse_stats.ach_oxtmin(indpulse) = y1mean(1001+inds);
            plot_areaerrorbar(y2',plotoptions)
            title(['Amp:' num2str(amp), ' Dur(s):' num2str(dur)])
            
            
        end
    end
    save([animalNames{anID},'_stats.mat'], 'dFF_pulse_stats');


    unique_numbers = unique(PulseData.pulse_dur);
    
    for jj = 1:length(unique_numbers)
        Dur_numbers(jj,:) = unique(ACh_OXT_Act.PulseData.pulse_id(PulseData.pulse_dur == unique_numbers(jj)));
    
    end


end

%%
load('ACh_OXT004_012_230328_141439ACh_OXT_Act.mat')
Ndff = length(ACh_OXT_Act.Ach_dFF);

% Test Plot
PulseData  = ACh_OXT_Act.PulseData;
lfp_fbtime = ACh_OXT_Act.lfp_fbtime;
Npulses = length(PulseData.trial_id);

%%
durs  = [1,5,10];
figure
set(gcf,'Position',[50 50 250 900])
for j = 1:3
ids = PulseData.pulse_amplitude == 6 & round(PulseData.pulse_dur) ==durs(j);
locs_insec = PulseData.pulse_start_intan(ids);
[ wavAvg, lfpAvg ] = bz_eventWavelet (ACh_OXT_Act.LFP, locs_insec,'twin', [10,10+durs(j)],'plotWave',false);

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);
subplot(3,1,j)
plot(wavAvg.timestamps,wavAvg_smoothed,'color', [0,0,0,0.25],LineWidth=1)

hold on
ids = PulseData.pulse_amplitude == 11 & round(PulseData.pulse_dur) ==durs(j);
locs_insec = PulseData.pulse_start_intan(ids);
[ wavAvg, lfpAvg ] = bz_eventWavelet (ACh_OXT_Act.LFP, locs_insec,'twin', [10,10+durs(j)],'plotWave',false);
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,'color', [0,0,0,0.5],LineWidth=1)

hold on
ids = PulseData.pulse_amplitude == 16 & round(PulseData.pulse_dur) ==durs(j);
locs_insec = PulseData.pulse_start_intan(ids);
[ wavAvg, lfpAvg ] = bz_eventWavelet (ACh_OXT_Act.LFP, locs_insec,'twin', [10,10+durs(j)],'plotWave',false);
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,'color', [0,0,0,0.7],LineWidth=1)

hold on
ids = PulseData.pulse_amplitude == 20 & round(PulseData.pulse_dur) ==durs(j);
locs_insec = PulseData.pulse_start_intan(ids);
[ wavAvg, lfpAvg ] = bz_eventWavelet (ACh_OXT_Act.LFP, locs_insec,'twin', [10,10+durs(j)],'plotWave',false);
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,'color', [0,0,0,1],LineWidth=1)
xlim([-3, 10+durs(j)])
end
exportgraphics(gcf ,'AchOXT005_RipplePower_Pulse.pdf','ContentType','vector')

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([40, 200])
%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(99)-wavAvg.freqs(87);
wavAvg_smoothed = smoothdata(mean(wavAvg2_band(:,87:99),2),'sgolay',1e3);
plot(wavAvg.timestamps,wavAvg_smoothed,'k')
%%  Group Data for Oxtmin
m1 = load('ACh_OXT004_012_stats.mat','dFF_pulse_stats');
m2 = load('ACh_OXT005_010_stats.mat','dFF_pulse_stats');
m4 = load('ACh_OXT007_009_stats.mat','dFF_pulse_stats');

OXTmin = [m1.dFF_pulse_stats.oxtmin/min(m1.dFF_pulse_stats.oxtmin); m2.dFF_pulse_stats.oxtmin/min(m2.dFF_pulse_stats.oxtmin);...
    m4.dFF_pulse_stats.oxtmin/min(m4.dFF_pulse_stats.oxtmin)];
figure
set(gcf,'Position',[50 50 500 400])

bar(-mean(OXTmin,1),'w')

hold on
scatter(1:12, -OXTmin,50,'k','o','filled')
errorbar(1:12, -mean(OXTmin,1), std(OXTmin,1), std(OXTmin,1), 'k', 'LineStyle','none','LineWidth',2)


Achmax = [m1.dFF_pulse_stats.ach_oxtmin/max(m1.dFF_pulse_stats.ach_oxtmin); m2.dFF_pulse_stats.ach_oxtmin/max(m2.dFF_pulse_stats.ach_oxtmin);...
    m4.dFF_pulse_stats.ach_oxtmin/max(m4.dFF_pulse_stats.ach_oxtmin)];

bar(mean(Achmax,1),'w')

hold on
scatter(1:12, Achmax,50,'k','o','filled')
errorbar(1:12, mean(Achmax,1), std(Achmax,1), std(Achmax,1), 'k', 'LineStyle','none','LineWidth',2)
ylim([-1.1,1.1])
box off
exportgraphics(gcf ,'Allmice_OxtminAch.pdf','ContentType','vector')


figure
set(gcf,'Position',[50 50 400 300])

Tmin = [m1.dFF_pulse_stats.min_time; m2.dFF_pulse_stats.min_time;m4.dFF_pulse_stats.min_time];
bar(mean(Tmin,1),'w')
hold on
scatter(1:12, Tmin,50,'k','o','filled')
errorbar(1:12, mean(Tmin,1), std(Tmin,1), std(Tmin,1), 'k', 'LineStyle','none','LineWidth',2)
box off
exportgraphics(gcf ,'Allmice_TminAch.pdf','ContentType','vector')

%% get data based on same duration with diffrecnt amplitude
data_sta = OXTmin;
Group_dur1 = data_sta(:,[1,4,7,10]);
Group_dur2 = data_sta(:,[2,5,8,11]);
Group_dur3 = data_sta(:,[3,6,9,12]);


Group_amp1 = data_sta(:,[1,5,9]);
Group_amp12 = data_sta(:,[2,6,10]);
Group_amp13 = data_sta(:,[3,7,11]);
Group_amp14 = data_sta(:,[4,8,12]);



%% RIPPLE
Fluo_Categorized_Ripple = Categorize_Ripple_Responses(FluoData,PulseData);
%%
pulse_types = Fluo_Categorized_Ripple.pulse_type(:,1:2);
Fluo = Fluo_Categorized_Ripple.dFF;
ripple_dt = Fluo_Categorized_Ripple.ripple_dt;
options.handle  = figure();
set(gcf,'Position',[50 50 250 900])
options.alpha      = 0.1;
options.line_width = 2;
options.error = 'sem';
criteria_dur = [1,5,10]; 

trial_bound = 1:30;
alphas= linspace(0.2,0.8,4);
for jj = 1:3
    subplot(3,1,jj)
    pulse_ids = find(pulse_types(:,2) ==criteria_dur(jj) );
    [~,II] = sort(pulse_types(pulse_ids,1));
    for i = 1:length(pulse_ids)
        ids = pulse_ids(II(i));
        eval(['xx = -10:ripple_dt:length(Fluo.Mixed.pulse' num2str(ids) '.Ripple(:,1))*ripple_dt-ripple_dt-10;']);
        options.x_axis = xx;
        eval(['y1 = Fluo.Mixed.pulse' num2str(ids) '.Ripple;']);
        options.color_area = [1,1,1];   
        options.color_line = [0,0,0,alphas(i)];
        plot_areaerrorbar(y1',options)
        ylim([0 0.8])
        xlim([-3,20])
        hold on

    end
end
exportgraphics(options.handle ,'Ripple_Pulse.pdf','ContentType','vector')
%%
pulse_id = [1,2,12];
options.handle  = figure();
options.alpha      = 0.1;
options.line_width = 2;
options.error = 'sem';
subplot(1,4,1)
for i = 1:3
    eval(['xx = -10:1/fs_fluo:length(Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Ach(:,1))/fs_fluo-1/fs_fluo-10;']);
    options.x_axis = xx;
    eval(['y1 = Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Ach;']);
    options.color_area = [255 0 0]./255;
    options.color_line = [ 255 0 0]./255;

    plot_areaerrorbar(y1',options)
    hold on
    options.color_area = [0 150 0]./255;
    options.color_line = [0 150 0]./255;
    eval(['y2 = Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Oxt;']);
    plot_areaerrorbar(y2',options)
    hold on
end
ylim([-0.015 0.04])

pulse_id = [4,7,9];
subplot(1,4,2)
for i = 1:3
    eval(['xx = -10:1/fs_fluo:length(Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Ach(:,1))/fs_fluo-1/fs_fluo-10;']);
    options.x_axis = xx;
    eval(['y1 = Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Ach;']);
    options.color_area = [255 0 0]./255;
    options.color_line = [ 255 0 0]./255;
    plot_areaerrorbar(y1',options)
    hold on
    options.color_area = [0 150 0]./255;
    options.color_line = [0 150 0]./255;
    eval(['y2 = Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Oxt;']);
    plot_areaerrorbar(y2',options)
    hold on
end
set(gca,'YTickLabel',[]);
ylim([-0.015 0.04])

pulse_id = [6,10,11];
subplot(1,4,3)
for i = 1:3
    eval(['xx = -10:1/fs_fluo:length(Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Ach(:,1))/fs_fluo-1/fs_fluo-10;']);
    options.x_axis = xx;
    eval(['y1 = Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Ach;']);
    options.color_area = [255 0 0]./255;
    options.color_line = [ 255 0 0]./255;
    plot_areaerrorbar(y1',options)
    hold on
    options.color_area = [0 150 0]./255;
    options.color_line = [0 150 0]./255;
    eval(['y2 = Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Oxt;']);
    plot_areaerrorbar(y2',options)
    hold on
end
set(gca,'YTickLabel',[]);
ylim([-0.015 0.04])

pulse_id = [5,8,13];
subplot(1,4,4)
for i = 1:3
    eval(['xx = -10:1/fs_fluo:length(Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Ach(:,1))/fs_fluo-1/fs_fluo-10;']);
    options.x_axis = xx;
    eval(['y1 = Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Ach;']);
    options.color_area = [255 0 0]./255;
    options.color_line = [ 255 0 0]./255;
    plot_areaerrorbar(y1',options)
    hold on
    options.color_area = [0 150 0]./255;
    options.color_line = [0 150 0]./255;
    eval(['y2 = Fluo.Mixed.pulse' num2str(pulse_id(i)) '.Oxt;']);
    plot_areaerrorbar(y2',options)
    hold on
end
set(gca,'YTickLabel',[]);
ylim([-0.015 0.04])