function   all_su_bursts=get_unit_bursts(data, min_ISI, min_spikes_per_burst)
%data - basic data struct including information regarding recording time, number of units, spike times etc.
%min_ISI - min inter-spike interval in ms for spikes to be considered a burst
%min_spikes_per_burst - min number of spikes for a burst to be accepted
%% preallocations
min_ISI=min_ISI/1000;
sep_unit_spikes = cell(data.maxUnitFullHours, data.nUnits);
all_burst_data = cell(data.maxUnitFullHours, data.nUnits);
mean_spikes_in_burst = zeros(data.maxUnitFullHours, data.nUnits);
mean_burst_duration = zeros(data.maxUnitFullHours, data.nUnits);
mean_in_burst_freq = zeros(data.maxUnitFullHours, data.nUnits);
mean_percent_spikes_in_bursts = zeros(data.maxUnitFullHours, data.nUnits);
mean_IBI = zeros(data.maxUnitFullHours, data.nUnits);
mean_short_bursts = zeros(data.maxUnitFullHours, data.nUnits);
mean_burst_freq = zeros(data.maxUnitFullHours, data.nUnits);
mean_norm_burst_freq= zeros(data.maxUnitFullHours, data.nUnits);
%% make cell array of spikes per unit per time
for t=1:data.maxUnitFullHours
    for u=1:data.nUnits
        [r] = find(data.unitSpikeTimes{u}(:)>1200*(t-1) & data.unitSpikeTimes{u}(:)<1200*t);
        sep_unit_spikes{t,u} = data.unitSpikeTimes{u}(r);
    end
end
%%
for t=1:data.maxUnitFullHours
    for u=1:data.nUnits
        spikes = sep_unit_spikes{t,u}(:);
        [burst_data] = single_burst_detect(spikes, min_ISI, min_spikes_per_burst);
        all_burst_data{t,u} = burst_data;
        mean_spikes_in_burst(t,u) = mean(burst_data.spikes_in_burst);
        mean_burst_duration(t,u) = mean(burst_data.burst_duration);
        mean_in_burst_freq(t,u) = mean(burst_data.in_burst_freq);
        mean_percent_spikes_in_bursts(t,u) = mean(burst_data.percent_spikes_in_bursts);
        mean_IBI(t,u) = mean(burst_data.IBI);
        mean_short_bursts(t,u) = mean(burst_data.short_bursts);
        mean_burst_freq(t,u) = mean(burst_data.burst_freq);
        mean_norm_burst_freq(t,u)= mean(burst_data.norm_burst_freq);
        
    end
end
%% remove non-bursting units
for u=1:data.nUnits
    if all(mean_burst_freq(:,u)<0.004)
        mean_spikes_in_burst(:,u) = NaN;
        mean_burst_duration(:,u) = NaN;
        mean_in_burst_freq(:,u) = NaN;
        mean_percent_spikes_in_bursts(:,u) = NaN;
        mean_IBI(:,u) = NaN;
        mean_short_bursts(:,u) = mean(burst_data.short_bursts);
        mean_burst_freq(:,u) = NaN;
        mean_norm_burst_freq(:,u)= NaN;
        
    end
end
all_su_bursts.mean_spikes_in_burst = mean_spikes_in_burst;
all_su_bursts.mean_burst_duration = mean_burst_duration;
all_su_bursts.mean_in_burst_freq = mean_in_burst_freq;
all_su_bursts.mean_percent_spikes_in_bursts = mean_percent_spikes_in_bursts;
all_su_bursts.mean_IBI = mean_IBI;
all_su_bursts.mean_burst_freq = mean_burst_freq;
all_su_bursts.mean_norm_burst_freq = mean_norm_burst_freq;
end


%%
function [burst_data]=single_burst_detect(spikes, min_ISI, min_spikes_per_burst)
ISI=diff(spikes);
[r]= find(ISI<=min_ISI);
bs=[];
for i=1:numel(r)
    bs(i)=r(i);
    if ~isequal(i,1)
        while isequal(r(i-1), (r(i)-1))
            bs(i)=0;
            
            i=i+1;
            if i>=numel(r)
                break
            end
        end
    end
end
if ~isempty(bs)
    bs_r=find(bs);
    bs=bs(bs_r);
    be=zeros(numel(bs),1);
    for i=1:numel(bs)
        if isempty(find(ISI(bs(i):numel(ISI))>min_ISI))
            be(i)=numel(ISI);
        else
            be(i)=(find(ISI(bs(i):numel(ISI))>min_ISI, 1, 'first'))+(bs(i))-1;
        end
    end
end
%%
%%single_burst_stats;
if ~isempty(bs)
    burst_spike_num = [bs' be];
    spikes_in_burst = diff(burst_spike_num,1,2)+1;
    [short_bursts, filt_burst_spike_num, filt_spikes_in_burst]= find_short_bursts(min_spikes_per_burst, burst_spike_num,  spikes_in_burst);
    burst_times(:,1) = spikes(filt_burst_spike_num(:,1));
    burst_times(:,2) = spikes(filt_burst_spike_num(:,2));
    burst_duration = diff(burst_times,1,2)*1000; %in ms when spikes are in seconds
    in_burst_freq = filt_spikes_in_burst./(burst_duration/1000); %in Hz
    percent_spikes_in_bursts = sum(filt_spikes_in_burst)/numel(spikes);
    burst_freq = numel(burst_duration)/1200;
    if isempty(burst_times) || size(burst_times,1)==1
        IBI=[];
    else
        for b=1:size(burst_times,1)-1
            IBI(b,1) = 1000*(burst_times(b+1,1)- burst_times(b,2));
        end
        
    end
    burst_data.burst_spike_num = filt_burst_spike_num;
    burst_data.burst_times = burst_times;
    burst_data.spikes_in_burst = filt_spikes_in_burst;
    burst_data.burst_duration = burst_duration;
    burst_data.in_burst_freq = in_burst_freq;
    burst_data.percent_spikes_in_bursts = percent_spikes_in_bursts;
    burst_data.IBI = IBI;
    burst_data.short_bursts=short_bursts;
    burst_data.burst_freq = burst_freq;
    burst_data.norm_burst_freq= burst_freq/(numel(spikes)/1200);
    
else
    burst_data.burst_spike_num = [];
    burst_data.burst_times = [];
    burst_data.spikes_in_burst = [];
    burst_data.burst_duration = [];
    burst_data.in_burst_freq = [];
    burst_data.percent_spikes_in_bursts = [];
    burst_data.IBI = [];
    burst_data.short_bursts = [];
    burst_data.burst_freq = 0;
    burst_data.norm_burst_freq= 0;
end
%%
%%
end
