%% File ID
% fid = '2012_03_05_Lollipop-10803';
r_idx = [1, 1, 1]; % Bursts to use in example rasters

fid = '2013_02_26_Yam-22565';
r_idx = [1, 1, 1]; % Bursts to use in example rasters

%% Load data
load([fid '_COMBO.SS'],'-mat')
SS.RemoveUnit(0);

%% Get burst shapes/times

% Timing parameters
dt = 0.01; % ASDR time bin
pre_sec = 0.5;% seconds
post_sec = 1;% seconds
minIBI = 10; % seconds

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

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

% Storage
bursts.dt = dt;
bursts.pre_sec = pre_sec;
bursts.post_sec = post_sec;
bursts.pre_drug.times = [];
bursts.pre_drug.waveforms = [];
bursts.drug.times = [];
bursts.drug.waveforms = [];
bursts.post_drug.times = [];
bursts.post_drug.waveforms = [];

% Pre-drug
steps = 0:3600:3*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
steps = 3*3600:3600:27*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.times = [bursts.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.drug.waveforms = [bursts.drug.waveforms; B];
    end
end

bursts.drug.waveforms = bursts.drug.waveforms/(numel(unique(SS.unit))-1);

% Post-drug
steps = 27*3600:3600:33*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'])

close all
figure('units','centimeter','position',[5 5 20 20]);
t = 1000*n*dt;

% Rasters
axes('units','centimeter','position',[5 10 6 3]);


% Burst profile
xlimit = [-250 500];
ylimit = [0 110];
axes('units','centimeter','position',[5 5 6 3]);
s = std(bursts.pre_drug.waveforms);
boundedline(t+20,mean(bursts.pre_drug.waveforms),s,'b','alpha');
s = std(bursts.drug.waveforms);
boundedline(t+20,mean(bursts.drug.waveforms),s,'k','alpha');
s = std(bursts.pre_drug.waveforms);
boundedline(t+20,mean(bursts.post_drug.waveforms),s,'r','alpha');


set(gca,'Xlim',xlimit,'Ylim',ylimit,'Fontsize',7,'FontName','Arial','TickDir','out')
saveas(gcf, [fid '_burst_shape'])
export_fig([fid '_burst_shape'],'-pdf','-transparent',gcf)

