
function [data, fdata, filtered_bursts, stats] = network_burst_analysis(data_string, active_units_percentile, min_num_spikes_in_burst, threshold_mode)

do_plot = 0;    
dt = 0.05; %minimal timespan
    dn = 2;  %number of minimal time intervals between two detected bursts s.t. they are considered one
    num_seconds_hour = 1200;
    
    l = (-5:5);
    gg = exp(-l.^2);
    gg = gg/sum(gg);
    lg = length(gg);
   % gg = [1 3 8 3 1]/16;
    data=new_loadData(data_string);
   
%% put all events on the same timeline

    N = size(data.unitSpikeTimes,2); 
    fdata = cell(N,1);
    spike_times = [];
    for i = 1:N,
       spike_times = [spike_times;  data.unitSpikeTimes{i}(:)]; 
       fdata{i} = data.unitSpikeTimes{i}(:)';
    end
    spike_times = sort(spike_times);
    Lt = length(spike_times);



    tmax = ceil(max(spike_times));
    tmin = floor(min(spike_times));
    Nsamp = (tmax - tmin)/dt;
    f = zeros(Nsamp,1);
     nsamp_hour = floor(num_seconds_hour/dt);
    Nhours = floor(Nsamp/nsamp_hour);

    thresh = calc_thresh(spike_times, Nhours, num_seconds_hour, gg, dt);

   
%% build a non-parametric histogram (Parzen) of the events 

     f = build_nonparam(spike_times, dt, Nsamp, gg);
    
%     imax = find(spike_times/dt < Nsamp-lg, 1, 'last');
%     imin = find(spike_times/dt > lg, 1);
%     
%     for i = imin:imax,
%        m0 = floor(spike_times(i)/dt-lg/2+1);
%        for k = 1:lg,
%            m = m0 + k;
%            f(m) = f(m) + gg(k);
%        end
%     end

 %% for each segment calculate the detection threshold; detect bursts
 %  calculate start and end of the bursts  
 
    burst_struct = cell(Nhours,1);
    burst_fin = [];
    for k = 1:Nhours,

        ff = f(((k-1)*nsamp_hour+1):k*nsamp_hour);
        %thresh = median(ff) + 1;

        ind_peak = find(ff > thresh(k));  
        tind = zeros(nsamp_hour,1);
        tind(ind_peak) = 1;
        tpeak = zeros(nsamp_hour,1);

        i = 1;
        idx_start = 1;
        idx_end = 1;
        burst_se = [];

        while (i < nsamp_hour-1)
           i = i+1;
           if (tind(i) == 1),
              idx_start = i;
              idx_end = i;
              while(tind(idx_end) == 1)
                  idx_end = idx_end+1;
                  if (idx_end == nsamp_hour)
                      break;
                  end
              end
              i = idx_end;
              burst_se = [burst_se; (k-1)*nsamp_hour + idx_start-1, (k-1)*nsamp_hour + idx_end];
           end
        end
        
%%  unite close bursts

        lb = size(burst_se,1);
        j = 1;
        while (j <= lb-1)
            js = burst_se(j,1);
            je = burst_se(j,2);
            while(j <= lb-1 && je >= burst_se(j+1,1) - dn)
                j = j+1;
                je = burst_se(j,2);
            end
            j = j+1;
            burst_fin = [burst_fin; js je];
        end
    end

    burst_fin = burst_fin*dt;
    filtered_bursts = [];

%% filter the detected bursts by the number of participating units

    Nb = size(burst_fin,1);
    active_struct = zeros(Nb, N);
    for n = 1:N,
        ts = fdata{n};
        KK = length(ts);
        k = 1;
        m = 1;
        while (k < KK & m < Nb)

            curr_sp = ts(k);
            curr_bst =  burst_fin(m,1);
            curr_be =  burst_fin(m,2);
            nsp = 0;

            while (k < KK & curr_sp < curr_bst)
                k = k+1;
                curr_sp = ts(k);
            end

            while (k < KK & curr_sp < curr_be)
                nsp = nsp+1;
                k = k+1;  
                curr_sp = ts(k);
            end

            active_struct(m,n) = (nsp >= min_num_spikes_in_burst);
            m = m+1;
        end

    end
    Num_active = sum(active_struct,2);
    burst_filter = zeros(size(Num_active));
    
    active_units = cell(Nhours,1);
    for k = 1:Nhours,      
        nstart = (k-1)*nsamp_hour+1;
        nend = k*nsamp_hour;
        burst_unfilt_idx = floor(find(burst_fin(:,1) > nstart*dt & burst_fin(:,2) < nend*dt));       
        active_units{k} = Num_active(burst_unfilt_idx);
                  
        [a,b] = hist(active_units{k}, 1:N);
        N_th = b(floor(find(cumsum(a/sum(a)) > active_units_percentile, 1)));
        burst_filter(burst_unfilt_idx) = (active_units{k} >  N_th); 
    end
    
    
    if (strcmp(threshold_mode, 'global_th') == 1)    
        [a,b] = hist(Num_active, 10);
        N_th = b(floor(find(cumsum(a/sum(a)) > active_units_percentile, 1)));
        burst_filter = (Num_active >  N_th);  
    end


    Num_active_filtered = [];
    for m = 1:Nb,
        if (burst_filter(m))
            filtered_bursts = [filtered_bursts; burst_fin(m,:)];
            Num_active_filtered  = [Num_active_filtered; Num_active(m)];
        end
    end

    filtered_bursts = refine_bursts(spike_times, filtered_bursts, 0.001);
    
%% plot the results
if do_plot == 1;
% figure(1)
    plot((1:length(f))*dt, f);
    hold on
    plot(filtered_bursts(:,1), f(floor(filtered_bursts(:,1)/dt)), 'go');
    plot(filtered_bursts(:,2), f(floor(filtered_bursts(:,2)/dt)), 'ro');
    hold off
end
    
 %% Statistics - to be completed 
 
    %stats = struct([]);
 
    g = cumsum(f);
    stats.mean_rate = zeros(Nhours,1);                  %mean rate -  for reference
    stats.norm_burst_freq = zeros(Nhours,1);            %burst frequency normalized to mean rate
    stats.fraction_spikes_in_bursts = zeros(Nhours,1);  %fraction of spikes in burst
    stats.spikes_per_burst = zeros(Nhours,1);           %number of spikes per burst
    stats.burst_duration = zeros(Nhours,1);             %duration of burst in secs
    stats.spikes_in_bursts_freq = zeros(Nhours,1);      %frequency of spikes in bursts
    stats.mean_frac_active_units = zeros(Nhours,1); 
    stats.mean_rate_units = zeros(Nhours, N);
    stats.active_units = active_units;

    
    for k = 1:Nhours,
               
        nstart = (k-1)*nsamp_hour+1;
        nend = k*nsamp_hour;
        gg = g(nstart:nend);
        for m = 1:N
            nspikes_unit = length(find(data.unitSpikeTimes{m}(:) > nstart*dt  & data.unitSpikeTimes{m}(:) < nend*dt));
            stats.mean_rate_units(k,m) = nspikes_unit/num_seconds_hour;
        end 
        burst_idx = floor(find(filtered_bursts(:,1) > nstart*dt & filtered_bursts(:,2) < nend*dt));       
        stats.mean_frac_active_units(k) = mean(active_units{k})/N;
             
        nbursts = length(burst_idx);
        stats.norm_burst_freq(k) = nbursts/(gg(end) - gg(1));
        
        stats.mean_rate(k) = (gg(end) - gg(1))/(N*num_seconds_hour);
        bursts_hour = filtered_bursts(burst_idx,:);
        delta_tb =  bursts_hour(:,2) - bursts_hour(:,1);
        stats.burst_duration(k) = mean(delta_tb);
        total_bursts_duration = sum(delta_tb);
        spikes_num_b = 0;
        stats.spikes_in_bursts_freq(k) = 0;
        for i = 1:nbursts,
           idx_start = floor(filtered_bursts(burst_idx(i),1)/dt);
           idx_end = floor(filtered_bursts(burst_idx(i),2)/dt);
           nspikes_b = sum(f(idx_start:idx_end)); 
           spikes_num_b = spikes_num_b + nspikes_b;
           stats.spikes_in_bursts_freq(k) = stats.spikes_in_bursts_freq(k) + (nspikes_b/Num_active_filtered(burst_idx(i)))/delta_tb(i);
        end
        stats.spikes_in_bursts_freq(k) = stats.spikes_in_bursts_freq(k)/nbursts;
        stats.spikes_per_burst(k) = spikes_num_b/nbursts;
        stats.fraction_spikes_in_bursts(k) = spikes_num_b/(gg(end) - gg(1));
        
    end

end

