% %% File ID
% 
fid = '2013_12_21_Smiley-24729';
% 
% %% Load data
load([fid '_COMBO.SS'],'-mat')
load([fid '_OC'])
SS.RemoveUnit(0);
% 
% %% Get burst shapes/times
% 
% Timing parameters
dt = 0.01; % ASDR time bin
pre_sec = 1;% seconds
post_sec = 3;% seconds
min_sep = 5; % seconds
minIBI = 10; % seconds

% Raw ASDR burst threshold
asdr_thresh = 1000; % 0.01;

% Burst snippets
pre_samp = ceil(pre_sec/dt);
post_samp = ceil(post_sec/dt);
min_sep_samp = min_sep/dt;
n = -pre_samp:post_samp;

% Stim snippets
pre_samp_stim = ceil(pre_sec*oc.meta.loop_rate_hz);
post_samp_stim = ceil(post_sec*oc.meta.loop_rate_hz);
n_stim = -pre_samp_stim:post_samp_stim;

% % Storage
% bursts.dt = dt;
% bursts.pre_sec = pre_sec;
% bursts.post_sec = post_sec;
% bursts.dt_stim = 1/oc.meta.loop_rate_hz;
% bursts.pre_drug.times = [];
% bursts.pre_drug.waveforms = [];
% bursts.drug.times = [];
% bursts.drug.waveforms = [];
% bursts.drug_noctl.times = [];
% bursts.drug_noctl.waveforms = [];
% bursts.drug.stim = [];
% bursts.drug.stim_i = [];
% bursts.post_drug.times = [];
% bursts.post_drug.waveforms = [];
% 
% % Pre-drug
% steps = 0:3600:1*3600;
% 
% for j = 1:numel(steps)-1
%     
%     %     close all
%     %     figure
%     %     hold on
%     SS.ASDR(dt,[steps(j) steps(j+1)],[],[],[],0);
%     
%     x = find(SS.asdr.asdr > asdr_thresh); %find threshold crossings
%     if ~isempty(x)
%         b_i = x(1); %index threshold crossings
%         c = 1;
%         for i = 1:length(x)
%             if SS.asdr.bin(x(i))-minIBI > SS.asdr.bin(b_i(c))
%                 c = c+1;
%                 b_i(c) = x(i); %ensure that there is a minimum interburst interval between bursts (otherwise 2nd burst is removed)
%             end
%         end;
%         b_i = b_i(b_i > pre_samp & b_i < numel(SS.asdr.bin) - post_samp);
%         bursts.pre_drug.times = [bursts.pre_drug.times; SS.asdr.bin(b_i)];
%         %     plot(SS.asdr.bin(b_i),SS.asdr.asdr(b_i),'r*')
%         
%         
%         B = zeros(numel(b_i),numel(n));
%         
%         for i = 1:length(b_i)
%             
%             idx = b_i(i) + n;
%             B(i,:) = SS.asdr.asdr(idx);
%             
%         end
%         
%         bursts.pre_drug.waveforms = [bursts.pre_drug.waveforms; B];
%     end
% end
% 
% 
% bursts.pre_drug.waveforms =  bursts.pre_drug.waveforms/(numel(unique(SS.unit))-1);
% 
% % Drug - control
% steps = [3600 1.166*3600];
% 
% for j = 1:numel(steps)-1
%     
%     %     close all
%     %     figure
%     %     hold on
%     SS.ASDR(dt,[steps(j) steps(j+1)],[],[],[],0);
%     
%     x = find(SS.asdr.asdr > asdr_thresh); %find threshold crossings
%     if ~isempty(x)
%         b_i = x(1); %index threshold crossings
%         c = 1;
%         for i = 1:length(x)
%             if SS.asdr.bin(x(i))-minIBI > SS.asdr.bin(b_i(c))
%                 c = c+1;
%                 b_i(c) = x(i); %ensure that there is a minimum interburst interval between bursts (otherwise 2nd burst is removed)
%             end
%         end;
%         b_i = b_i(b_i > pre_samp & b_i < numel(SS.asdr.bin) - post_samp);
%         bursts.drug_noctl.times = [bursts.drug_noctl.times; SS.asdr.bin(b_i)];
%         %     plot(SS.asdr.bin(b_i),SS.asdr.asdr(b_i),'r*')
%         
%         
%         B = zeros(numel(b_i),numel(n));
%         
%         for i = 1:length(b_i)
%             
%             idx = b_i(i) + n;
%             B(i,:) = SS.asdr.asdr(idx);
%             
%         end
%         
%         bursts.drug_noctl.waveforms = [bursts.drug_noctl.waveforms; B];
%     end
% end
% 
% 
% bursts.drug_noctl.waveforms =  bursts.drug_noctl.waveforms/(numel(unique(SS.unit))-1);
% 
% % Drug + control
% steps = (1.166*3600):3600:(4.166*3600);
% 
% for j = 1:numel(steps)-1
%     
%     close all
%     figure
%     hold on
%     SS.ASDR(dt,[steps(j) steps(j+1)],[],[],[],1);
%     
%     x = find(SS.asdr.asdr > asdr_thresh); %find threshold crossings
%     if ~isempty(x)
%         b_i = x(1); %index threshold crossings
%         c = 1;
%         for i = 1:length(x)
%             if SS.asdr.bin(x(i))-minIBI > SS.asdr.bin(b_i(c))
%                 c = c+1;
%                 b_i(c) = x(i); %ensure that there is a minimum interburst interval between bursts (otherwise 2nd burst is removed)
%             end
%         end;
%         b_i = b_i(b_i > pre_samp & b_i < numel(SS.asdr.bin) - post_samp);
%         
%         s_i = zeros(size(b_i));
%         for i = 1:numel(b_i)
%             [~, s_i(i)] = min(abs(oc.time + 3600 - SS.asdr.bin(b_i(i))));
%         end
%         
%         bursts.drug.times = [bursts.drug.times; SS.asdr.bin(b_i)];
%         plot(SS.asdr.bin(b_i),SS.asdr.asdr(b_i),'r*')
%         plot(oc.time(s_i),1000*ones(size(s_i)),'g*')
%         
%         B = zeros(numel(b_i),numel(n));
%         S = zeros(numel(b_i),numel(n_stim));
%         I = zeros(numel(b_i),numel(n_stim));
%         
%         for i = 1:length(b_i)
%             
%             idx = b_i(i) + n;
%             B(i,:) = SS.asdr.asdr(idx);
%             
%             idx = s_i(i) + n_stim;
%             I(i,:) = s_i(i);
%             S(i,:) = oc.out2(idx);
%             
%             
%         end
%         
%         bursts.drug.waveforms = [bursts.drug.waveforms; B];
%         bursts.drug.stim = [bursts.drug.stim; S];
%         bursts.drug.stim_i = [bursts.drug.stim_i; I];
%     end
% end
% 
% bursts.drug.waveforms = bursts.drug.waveforms/(numel(unique(SS.unit))-1);
% 
% % Post-drug
% steps = 4*3600:3600:10*3600;
% 
% for j = 1:numel(steps)-1
%     
%     %     close all
%     %     figure
%     %     hold on
%     SS.ASDR(dt,[steps(j) steps(j+1)],[],[],[],0);
%     
%     x = find(SS.asdr.asdr > asdr_thresh); %find threshold crossings
%     if ~isempty(x)
%         b_i = x(1); %index threshold crossings
%         c = 1;
%         for i = 1:length(x)
%             if SS.asdr.bin(x(i))-minIBI > SS.asdr.bin(b_i(c))
%                 c = c+1;
%                 b_i(c) = x(i); %ensure that there is a minimum interburst interval between bursts (otherwise 2nd burst is removed)
%             end
%         end;
%         b_i = b_i(b_i > pre_samp & b_i < numel(SS.asdr.bin) - post_samp);
%         bursts.post_drug.times = [bursts.post_drug.times; SS.asdr.bin(b_i)];
%         %     plot(SS.asdr.bin(b_i),SS.asdr.asdr(b_i),'r*')
%         
%         
%         B = zeros(numel(b_i),numel(n));
%         for i = 1:length(b_i)
%             
%             idx = b_i(i) + n;
%             B(i,:) = SS.asdr.asdr(idx);
%             
%         end
%         
%         bursts.post_drug.waveforms = [bursts.post_drug.waveforms; B];
%     end
%     
% end
% 
% bursts.post_drug.waveforms = bursts.post_drug.waveforms/(numel(unique(SS.unit))-1);
% save([fid '_BURSTS'],'bursts')

%% Plots
load([fid '_BURSTS'])

% 53.534027778 mm total height

dt_shift = 0.020; % s (there is no true burst 'start')

close all
figure('units','centimeter','position',[5 5 20 20]);
xlimit = [-0.5 1.0]; % s
% xlimit = [-1 3]; % s


% Rasters
ylimit = [1 max(SS.unit)];
axes('units','centimeter','position',[5 13+0.60 2.5 0.5]);
idx = 4;
b = SS.time >= bursts.pre_drug.times(idx) - pre_sec & SS.time <= bursts.pre_drug.times(idx) + post_sec; 
plot((SS.time(b) - bursts.pre_drug.times(idx)) + dt_shift,SS.unit(b),'.','Color',[57 83 164]/255,'MarkerSize',2);
set(gca,'Xlim',xlimit,'Ylim',ylimit,'XTick',[],'Fontsize',7,'FontName','Arial','TickDir','out','box','off')

ylimit = [1 max(SS.unit)];
axes('units','centimeter','position',[5 13 2.5 0.5]);
idx = 4;
b = SS.time >= bursts.drug_noctl.times(idx) - pre_sec & SS.time <= bursts.drug_noctl.times(idx) + post_sec; 
plot((SS.time(b) - bursts.drug_noctl.times(idx)) + dt_shift,SS.unit(b),'.','Color',[217 30 38]/255,'MarkerSize',2);
set(gca,'Xlim',xlimit,'Ylim',ylimit,'XTick',[],'Fontsize',7,'FontName','Arial','TickDir','out','box','off')

axes('units','centimeter','position',[5 13-0.60 2.5 0.5]);
hold on
idx = 120;
b = SS.time >= bursts.drug.times(idx) - pre_sec & SS.time <= bursts.drug.times(idx) + post_sec; 
plot((SS.time(b) - bursts.drug.times(idx)) + dt_shift,SS.unit(b),'.','Color',[16 169 74]/255,'MarkerSize',2);

t = (1:numel(bursts.drug.stim(idx,:)))/oc.meta.loop_rate_hz;
stairs(t - pre_sec + dt_shift,ylimit(2)*bursts.drug.stim(idx,:))
set(gca,'Xlim',xlimit,'Ylim',ylimit,'XTick',[],'Fontsize',7,'FontName','Arial','TickDir','out','box','off')

axes('units','centimeter','position',[5 13-1.2 2.5 0.5]);
idx = 143;
b = SS.time >= bursts.post_drug.times(idx) - pre_sec & SS.time <= bursts.post_drug.times(idx) + post_sec; 
plot((SS.time(b) - bursts.post_drug.times(idx)) + dt_shift,SS.unit(b),'.','Color',[255 204 4]/255,'MarkerSize',2);
set(gca,'Xlim',xlimit,'Ylim',ylimit,'XTick',[],'Fontsize',7,'FontName','Arial','TickDir','out','box','off')

% Burst profile
t = n*bursts.dt;
ylimit = [0 200];
axes('units','centimeter','position',[5 9.9 2.5 1.75]);
hold on

stairs(t+dt_shift,mean(bursts.pre_drug.waveforms),'Color',[57 83 164]/255);
stairs(t+dt_shift,mean(bursts.drug_noctl.waveforms),'Color',[217 30 38]/255);
stairs(t+dt_shift,mean(bursts.drug.waveforms),'Color',[16 169 74]/255);
stairs(t+dt_shift,mean(bursts.post_drug.waveforms),'Color',[255 204 4]/255);

% s = std(bursts.pre_drug.waveforms);
% boundedline(t+dt_shift,mean(bursts.pre_drug.waveforms),s,'b','alpha');
% s = std(bursts.drug_noctl.waveforms);
% boundedline(t+dt_shift,mean(bursts.drug_noctl.waveforms),s,'r','alpha');
% s = std(bursts.drug.waveforms);
% boundedline(t+dt_shift,mean(bursts.drug.waveforms),s,'g','alpha');
% s = std(bursts.pre_drug.waveforms);
% boundedline(t+dt_shift,mean(bursts.post_drug.waveforms),s,'y','alpha');

set(gca,'Xlim',xlimit,'Ylim',ylimit,'Fontsize',7,'FontName','Arial','TickDir','out','box','off','XTickLabel',{}, 'YScale','linear')


% Burst-triggered control signal average
t = n_stim*bursts.dt_stim;
u = 10.5*mean(bursts.drug.stim);
s = 10.5*std(bursts.drug.stim);
ylimit = [4 11];

axes('units','centimeter','position',[5 8.75 2.5 1]);
plot(t+dt_shift,u,'k');
set(gca,'Xlim',xlimit,'Ylim',ylimit,'Fontsize',7,'FontName','Arial','TickDir','out','box','off')

% Save
saveas(gcf, [fid '_burst_shape'])
export_fig([fid '_burst_shape'],'-pdf','-transparent',gcf)




