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')

%% Add Util Files in to your folder
FolderName = cd;
AboveFolder = fullfile(FolderName, '..');
addpath(fullfile(AboveFolder,'Ach_Oxt_Utils'));

addpath('G:\My Drive\Papers_Buzsaki\Paper_ACh_OXT\Ach_Oxt_Utils')

%% Get FluoData
FluoData = Get_OxtAch_Fluos();
%% Get StimPulse data
PulseData = CreatePulseTable();

%% Get pulse and state categorized dFFs
Fluo_Categorized = Categorize_Responses(FluoData,PulseData);

%%
save('CIRP001_004_processed.mat', 'Fluo_Categorized', 'PulseData', 'FluoData');
%% Test Plot
fs_fluo = Fluo_Categorized.fs;
Fluo = Fluo_Categorized.dFF;
options.handle  = figure();
options.error = 'sem';
options.line_width = 2;
options.color_area = [255 0 0]./255;   
options.color_line = [ 255 0 0]./255;
options.alpha      = 0.5;

for i = 1:15
    subplot(5,3,i)
     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;']);
    
    options.x_axis = xx;
    eval(['y1 = Fluo.Mixed.pulse' num2str(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 140 0]./255;    
    options.color_line = [0 140 0]./255;
    eval(['y2 = Fluo.Mixed.pulse' num2str(i) '.Oxt;']);
    plot_areaerrorbar(y2',options)
    title(['Amp:' num2str(amp), ' Dur(s):' num2str(dur)])
end


%%
pulse_id = [1,2,12];
options.handle  = figure();
options.alpha      = 0.1;
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])

%%
pulse_types = Fluo_Categorized.pulse_type(:,1:2);
options.handle  = figure();
set(gcf,'Position',[50 50 250 900])
options.alpha      = 0.1;
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:1/fs_fluo:length(Fluo.Mixed.pulse' num2str(ids) '.Ach(:,1))/fs_fluo-1/fs_fluo-10;']);
        options.x_axis = xx;
        eval(['y1 = Fluo.Mixed.pulse' num2str(ids) '.Ach;']);
        options.color_area = [1,0,0];   
        options.color_line = [1,0,0,alphas(i)];
        plot_areaerrorbar(y1(:,trial_bound)',options)

        hold on
        options.color_area = [0,0.8,0];    
        options.color_line = [0,0.6,0,alphas(i)];
        eval(['y2 = Fluo.Mixed.pulse' num2str(ids) '.Oxt;']);
        plot_areaerrorbar(y2(:,trial_bound)',options)
        hold on
        ylim([-0.02 0.06])
        xlim([-3,20])

        %
        idend = round(criteria_dur(jj)*fs_fluo);
        x1 = xx(1001:1101);
        ym = mean(y1,2);
        ym= ym(1001:1101)';
        [fitted] = RiseFit(x1,ym);
        Trise_fit(jj,i) = fitted.riserate.*fitted.amplitude;

        %
        paramsout =FindPeak_Trough(y1(:,trial_bound),y2(:,trial_bound));
        Ach_max_mean(jj,i) = paramsout.Ach_max_mean;
        Ach_max_sem(jj,i) = paramsout.Ach_max_sem;
        Oxt_min_mean(jj,i) = paramsout.Oxt_min_mean;
        Oxt_min_sem(jj,i) = paramsout.Oxt_min_sem;
        Tmin(jj,i) = paramsout.min_time;

        %
        y1_mean = mean(y1(idend+980:idend+1000,:),1);
        y2_mean = mean(y2(idend+980:idend+1000,:),1);
        Ach_end_mean(jj,i) = mean(y1_mean(trial_bound));
        Ach_end_sem(jj,i) = std(y1_mean(trial_bound))/sqrt(length(trial_bound));
        Oxt_end_mean(jj,i) = mean(y2_mean(trial_bound));
        Oxt_end_sem(jj,i) = std(y2_mean(trial_bound))/sqrt(length(trial_bound));
%         figure(2)
%         subplot(3,1,jj)
%         plot(x1,ym, 'ko',x1,fitted.fx','b');
%         hold on

    end
end
%exportgraphics(options.handle ,'PulseAmplitude_first30trials.pdf','ContentType','vector')
%%
color_ = ['k','b','m','r'];
for kk = 1:4
    scatter(abs(Oxt_min_mean(:,kk)),Ach_max_mean(:,kk));
    hold on;
    errorbar(abs(Oxt_min_mean(:,kk)),Ach_max_mean(:,kk),Ach_max_sem(:,kk),Ach_max_sem(:,kk), 'linestyle', 'none', 'color', color_(kk), 'linewidth', 2);
    errorbar(abs(Oxt_min_mean(:,kk)),Ach_max_mean(:,kk),Oxt_min_sem(:,kk),Oxt_min_sem(:,kk),'horizontal', 'linestyle', 'none', 'color', color_(kk), 'linewidth', 2);
    hold on
end

%% BAR PLOT
 
PPdata = Oxt_min_mean';
PPerr = Oxt_min_sem';
h = figure
set(gcf,'Position',[50 50 500 900])
subplot(2,1,2)
hBar = bar(PPdata);
hBar(1).FaceColor = [.2,.8,0];
hBar(2).FaceColor = [.2 .6 .0];
hBar(3).FaceColor = [.2 .4 .00];

ctr = [];
ydt = [];
for k1 = 1:size(PPdata,2)
    ctr(k1,:) = bsxfun(@plus, hBar(k1).XData, hBar(k1).XOffset');   % Note: ‘XOffset’ Is An Undocumented Feature, This Selects The ‘bar’ Centres
    ydt(k1,:) = hBar(k1).YData;                                     % Individual Bar Heights
end
hold on
errorbar(ctr, ydt, PPerr.', "s","MarkerSize",8,...
    "MarkerEdgeColor","blue","MarkerFaceColor",'black','color', 'k', 'linewidth', 2)
PPx=categorical({'\it P0/18','\it P0/4','\it P0/2','\it P0'});
set(gca, 'XTickLabel', PPx)
ylim([-0.02 0])
title('Oxytocin minimum')
%
subplot(2,1,1)
PPdata = Ach_max_mean';
PPerr = Ach_max_sem';
hBar = bar(PPdata);
hBar(3).FaceColor = [0.5 0 0];
hBar(2).FaceColor = [0.75 0 0];
hBar(1).FaceColor = [1 0 0];
ctr = [];
ydt = [];
for k1 = 1:size(PPdata,2)
    ctr(k1,:) = bsxfun(@plus, hBar(k1).XData, hBar(k1).XOffset');   % Note: ‘XOffset’ Is An Undocumented Feature, This Selects The ‘bar’ Centres
    ydt(k1,:) = hBar(k1).YData;                                     % Individual Bar Heights
end
hold on
errorbar(ctr, ydt, PPerr.', "s","MarkerSize",8,...
    "MarkerEdgeColor","blue","MarkerFaceColor",'black','color', 'k', 'linewidth', 2)
PPx=categorical({'\it P0/18','\it P0/4','\it P0/2','\it P0'});
set(gca, 'XTickLabel', PPx)
title('Ach maximum')
%legend({'1 s','5 s','10 s'},'Location','southoutside')
ylim([0,0.035])

%exportgraphics(h ,'Ach_Oxt_minmax.pdf','ContentType','vector')

%% Tmin PLot

h = figure()
hBar = bar(Tmin');
hBar(1).FaceColor = [.2,.8,0];
hBar(2).FaceColor = [.2 .6 .0];
hBar(3).FaceColor = [.2 .4 .00];
exportgraphics(h ,'Oxt_mintime.pdf','ContentType','vector')

%% BAR PLOT for End of Stimulation
 
PPdata = Oxt_end_mean';
PPerr = Oxt_end_sem';
h = figure
set(gcf,'Position',[50 50 500 900])
subplot(2,1,2)
hBar = bar(PPdata);
hBar(1).FaceColor = [.2,.8,0];
hBar(2).FaceColor = [.2 .6 .0];
hBar(3).FaceColor = [.2 .4 .00];

ctr = [];
ydt = [];
for k1 = 1:size(PPdata,2)
    ctr(k1,:) = bsxfun(@plus, hBar(k1).XData, hBar(k1).XOffset');   % Note: ‘XOffset’ Is An Undocumented Feature, This Selects The ‘bar’ Centres
    ydt(k1,:) = hBar(k1).YData;                                     % Individual Bar Heights
end
hold on
errorbar(ctr, ydt, PPerr.', "s","MarkerSize",8,...
    "MarkerEdgeColor","blue","MarkerFaceColor",'black','color', 'k', 'linewidth', 2)
PPx=categorical({'\it P0/18','\it P0/4','\it P0/2','\it P0'});
set(gca, 'XTickLabel', PPx)
ylim([-0.01 0.01])
title('Oxytocin ')
%
subplot(2,1,1)
PPdata = Ach_end_mean';
PPerr = Ach_end_sem';
hBar = bar(PPdata);
hBar(3).FaceColor = [0.5 0 0];
hBar(2).FaceColor = [0.75 0 0];
hBar(1).FaceColor = [1 0 0];
ctr = [];
ydt = [];
for k1 = 1:size(PPdata,2)
    ctr(k1,:) = bsxfun(@plus, hBar(k1).XData, hBar(k1).XOffset');   % Note: ‘XOffset’ Is An Undocumented Feature, This Selects The ‘bar’ Centres
    ydt(k1,:) = hBar(k1).YData;                                     % Individual Bar Heights
end
hold on
errorbar(ctr, ydt, PPerr.', "s","MarkerSize",8,...
    "MarkerEdgeColor","blue","MarkerFaceColor",'black','color', 'k', 'linewidth', 2)
PPx=categorical({'\it P0/18','\it P0/4','\it P0/2','\it P0'});
set(gca, 'XTickLabel', PPx)
title('Ach')
%legend({'1 s','5 s','10 s'},'Location','southoutside')
ylim([0,0.06])

%exportgraphics(h ,'Ach_Oxt_EndofStimulation.pdf','ContentType','vector')

%%
figure;
h = figure()
set(gcf,'Position',[50 50 500 900])
subplot(2,1,1)
hBar = bar(Trise_fit');
hBar(3).FaceColor = [0.5 0 0];
hBar(2).FaceColor = [0.75 0 0];
hBar(1).FaceColor = [1 0 0];

PPdata = Oxt_min_mean';
PPerr = Oxt_min_sem';
hold on
subplot(2,1,2)
hBar = bar(PPdata);
hBar(1).FaceColor = [.2,.8,0];
hBar(2).FaceColor = [.2 .6 .0];
hBar(3).FaceColor = [.2 .4 .00];


%% 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;
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')

