function [outputTable, outputBursts] = Times2Stats(event_times, Max_Time)
%TIMES2STATS Convert the spike list into the statistical summary
%   event_times should be a horizontal list of the times at which events were detected on a single electrode.
outputTable = table();

% Information on Axions standard output.
% outputTable.('Measurement') = missing;
outputTable.('Number of Spikes') = NaN;
outputTable.('Mean Firing Rate (Hz)') = NaN;
outputTable.('Mean Firing Rate Alt (Hz)') = NaN;
outputTable.('ISI') = NaN;
outputTable.('ISI Coefficient of Variation') = NaN;
outputTable.('Number of Bursts') = NaN;
outputTable.('Burst Duration - Avg (s)') = NaN;
outputTable.('Burst Duration - Std (s)') = NaN;
outputTable.('Number of Spikes per Burst - Avg') = NaN;
outputTable.('Number of Spikes per Burst - Std') = NaN;
outputTable.('Mean ISI within Burst - Avg') = NaN;
outputTable.('Mean ISI within Burst - Std') = NaN;
outputTable.('Median ISI within Burst - Avg') = NaN;
outputTable.('Median ISI within Burst - Std') = NaN;
outputTable.('Inter-Burst Interval - Avg (s)') = NaN;
outputTable.('Inter-Burst Interval - Std (s)') = NaN;
outputTable.('Burst Frequency (Hz)') = NaN;
outputTable.('IBI Coefficient of Variation') = NaN;
outputTable.('Normalized Duration IQR') = NaN;
outputTable.('Burst Percentage') = NaN;

outputBursts = Intervals2Bursts([], []);

if ~isempty(event_times)
    % Number of Spikes, Average Frequency
    event_times(event_times>Max_Time) = []; % Limit the spike times to the maximum time passed to the function.

    outputTable.('Number of Spikes') = max(size(event_times));
    if exist('Max_Time', 'var')
        outputTable.('Mean Firing Rate (Hz)') = (outputTable.('Number of Spikes')) / Max_Time; 
    else
        outputTable.('Mean Firing Rate (Hz)') = (outputTable.('Number of Spikes')) / max(event_times); % Average frequency = Number of events -1 devided by the time between first and last event.
    end

    % Inter-spike-intervals (ISIs) (Periods)
    ISIs = Times2Periods(event_times);      % Find the inter-spike-intervals.
    ISI_Stats = Periods2Stats(ISIs);        % Get the basic stats from the ISIs.
    ISI_Stats = renamevars(ISI_Stats,...    % Rename the output variables.
        ["N",                   "f_avg",                        "i_avg",    "i_COV"],...
        ["Number of Spikes",    "Mean Firing Rate Alt (Hz)",    "ISI",      "ISI Coefficient of Variation"]);
    ISI_Stats = removevars(ISI_Stats,...    % Remove unnecessary variables.
        ["f_std", "f_COV", "f_med", "i_std", "i_med"]);
    outputTable(:, ISI_Stats.Properties.VariableNames) = ISI_Stats;
%     outputTable.('Mean Firing Rate (Hz)') = (outputTable.('Number of Spikes')) / max(event_times); % Average frequency = Number of events devided by the time of the last event.
    
    % Burst Stats
    outputBursts = Intervals2Bursts(ISIs, event_times); % Generate a table with all of the bursts found.
    BurstStats = Bursts2Stats(outputBursts);            % Generate summary stats of those bursts.
    
%     outputTable = [outputTable BurstStats]; % Add the burst stats to the output table.
    outputTable(:, BurstStats.Properties.VariableNames) = BurstStats; % Add the burst stats to the output table.
    outputTable.('Burst Percentage') = 100 * outputTable.('Number of Bursts') * outputTable.('Number of Spikes per Burst - Avg') / outputTable.('Number of Spikes'); % For consistency with the Axis Summary, the percentage is used here. (*100)
end
% pause; % For debugging purposes.

end


function [intervals] = Times2Periods(event_times) % Find the periods between events. 
    intervals = diff([NaN event_times]); % Treat the first event interval as NaN. Since it is unknown when the previous spike occured.
end


function [outputTable] = Periods2Stats(intervals) % Find the basic stats from a set of intervals.
% intervals is a list of the intervals between events in either orientation. 1*n or n*1
    frequencies = intervals.^-1;
    outputTable = table();

    outputTable.N     = max(size(intervals));                   % Number of events being analysed.
    outputTable.f_avg = mean(frequencies, 'omitnan');           % Find the mean frequency.
    outputTable.f_std = std(frequencies, 0, 'omitnan');         % Find the related std.
    outputTable.f_COV = outputTable.f_std / outputTable.f_avg;  % Find the frequency Coefficient of Variation. std / avg.
    outputTable.f_med = median(frequencies, 'omitnan');         % Find the median frequency.
    outputTable.i_avg = mean(intervals, 'omitnan');             % Find the mean interval.
    outputTable.i_std = std(intervals, 0, 'omitnan');           % Find the related std.
    outputTable.i_COV = outputTable.i_std / outputTable.i_avg;  % Find the interval Coefficient of Variation. std / avg.
    outputTable.i_med = median(intervals, 'omitnan');           % Find the median interval.
end


function [burst_table] = Intervals2Bursts(ISIs, event_times)    % Detect bursts and store them in a table.
% ISIs is the unflitered interspike interval list, and event times are the times that each interval was completed.
    burst_table = table();
    
    Burst_Threshold = 0.1;  % Time duration (s) for two events to happen to be counted as a possible burst.
    Burst_MinEvents = 5;    % Number of Burst threshold events required to be counted as a burst.
    
    BurstFlags = ISIs <= Burst_Threshold;               % Create an index of events that qualify as bursts.
    BurstFlags(BurstFlags ~= 1) = 0;                    % Force non-true flags to 0. In case of NaN.
%     BurstFlags = BurstFlags & [BurstFlags(2:end) 1];    % Include the preceding spike as a burst in all cases.
    
    BurstStarts = strfind([0 BurstFlags 0], [0 1]) - 1; % Find where the potential bursts Start. strfind technically finds the second spike, not the first in the burst.
    BurstEnds   = strfind([0 BurstFlags 0], [1 0]) - 1; % Find where the potential bursts End. The strfind is actually shifted by 1 position, hence the -1 adjustment.
    
    BurstLengths = BurstEnds - BurstStarts + 1;         % Find number of spikes per burst.
    
    BurstStarts( BurstLengths < Burst_MinEvents) = [];  % Remove potential bursts that do not have enough events.
    BurstEnds(   BurstLengths < Burst_MinEvents) = [];  
    BurstLengths(BurstLengths < Burst_MinEvents) = [];
    
    burst_table.B_Start = BurstStarts';                 % Store the qualifying Burst Starts.
    burst_table.B_End   = BurstEnds';                   % Store the qualifying Burst Ends.
    burst_table.B_Len   = BurstLengths';                % Store the qualifying Burst Lengths.
    
    burst_count = size(BurstStarts, 2);                 % Find number of bursts.

    burst_table.B_Time      = NaN(burst_count, 1);      % Default/placeholder for Burst time.
    burst_table.B_Duration  = NaN(burst_count, 1);      % Default/placeholder for Burst duration.
    burst_table.B_N_Spikes  = NaN(burst_count, 1);      % Default/placeholder for Burst number of spikes.
    burst_table.B_iISI_Avg  = NaN(burst_count, 1);      % Default/placeholder for Burst interspike interval average.
    burst_table.B_iISI_Med  = NaN(burst_count, 1);      % Default/placeholder for Burst interspike interval median.

    if burst_count                                                      % Process stats if a burst has been detected.
        for i_b = 1:burst_count                                         % Loop through each burst to determine stats.

            B_Time(i_b)     = event_times(BurstStarts(i_b));            % Find the time the burst starts.
            B_Index(i_b)    = {BurstStarts(i_b):BurstEnds(i_b)};        % Find all spikes associated with the burst.
            B_Duration(i_b) = event_times(BurstEnds(i_b)) - event_times(BurstStarts(i_b)); % Find  the duration of the burst.

            B_Stats = Periods2Stats([NaN ISIs(B_Index{i_b}(2:end))]);   % Replace the first ISI with a NaN, as the first spike in the burst does not have a suitable ISI.

            B_N_Spikes(i_b) = B_Stats.N;                % Extract the Number of Spikes in the burst.
            B_iISI_Avg(i_b) = B_Stats.i_avg;            % Extract the interspike interval average.
            B_iISI_Med(i_b) = B_Stats.i_med;            % Extract the interspike interval median.
        end

        burst_table.B_Time          = B_Time';      % Store the burst time in the output burst table.
        burst_table.B_Duration      = B_Duration';  % Store the burst duration in the output burst table.
        burst_table.B_N_Spikes      = B_N_Spikes';  % Store the burst number of spikes in the output burst table.
        burst_table.B_iISI_Avg      = B_iISI_Avg';  % Store the burst interspike interval average in the output burst table.
        burst_table.B_iISI_Med      = B_iISI_Med';  % Store the burst interspike interval median in the output burst table.
    end
end

function outputTable = Bursts2Stats(BurstTable) % Generates the statistics from sets of bursts.
% BurstTable is expected to be a list of the individual statistics of each burst.
    outputTable = table();
    
    N_Bursts = size(BurstTable, 1);
    outputTable.('Number of Bursts')    = N_Bursts;
% Default/placeholders for Burst Stats.
    outputTable.('Burst Duration - Avg (s)')            = NaN(1);
    outputTable.('Burst Duration - Std (s)')            = NaN(1);
    outputTable.('Number of Spikes per Burst - Avg')    = NaN(1);
    outputTable.('Number of Spikes per Burst - Std')    = NaN(1);
    outputTable.('Mean ISI within Burst - Avg')         = NaN(1);
    outputTable.('Mean ISI within Burst - Std')         = NaN(1);
    outputTable.('Median ISI within Burst - Avg')       = NaN(1);
    outputTable.('Median ISI within Burst - Std')       = NaN(1);
    outputTable.('Inter-Burst Interval - Avg (s)')      = NaN(1);
    outputTable.('Inter-Burst Interval - Std (s)')      = NaN(1);
    outputTable.('Burst Frequency (Hz)')                = NaN(1);
    outputTable.('IBI Coefficient of Variation')        = NaN(1);
    outputTable.('Normalized Duration IQR')             = NaN(1);
    
    if N_Bursts                                                                             % Produce stats if bursts have been found.
        outputTable.('Burst Duration - Avg (s)')            = mean(BurstTable.B_Duration);  %   Average burst duration.
        outputTable.('Burst Duration - Std (s)')            = std(BurstTable.B_Duration);   %   Standard deviation of Burst duration.
        outputTable.('Number of Spikes per Burst - Avg')    = mean(BurstTable.B_N_Spikes);  %   Average number of spikes per burst.
        outputTable.('Number of Spikes per Burst - Std')    = std(BurstTable.B_N_Spikes);   %   Standard deviation of spike number per burst.
        outputTable.('Mean ISI within Burst - Avg')         = mean(BurstTable.B_iISI_Avg);  %   Average (Average (interspike interval)), within each burst.
        outputTable.('Mean ISI within Burst - Std')         = std(BurstTable.B_iISI_Avg);   %       Associated Standard deviation.
        outputTable.('Median ISI within Burst - Avg')       = mean(BurstTable.B_iISI_Med);  %   Average (Median (interspike interval)), within each burst.
        outputTable.('Median ISI within Burst - Std')       = std(BurstTable.B_iISI_Med);   %       Associated Standard deviation.
        
        IBIs    = Times2Periods(BurstTable.B_Time');                                        % Find the interburst periods.
        B_Stats = Periods2Stats(IBIs);                                                      % Produce stats on the interburst periods.
        
        outputTable.('Inter-Burst Interval - Avg (s)')      = B_Stats.i_avg;                % Store the Interburst interval average.
        outputTable.('Inter-Burst Interval - Std (s)')      = B_Stats.i_std;                % Store the Interburst interval standard deviation.
        outputTable.('Burst Frequency (Hz)')                = B_Stats.f_avg;                % Store the bursting frequency average.
        outputTable.('IBI Coefficient of Variation')        = B_Stats.i_COV;                % Store the Interburst interval coefficient of variation.
        outputTable.('Normalized Duration IQR')             = iqr(BurstTable.B_Duration) / median(BurstTable.B_Duration); % After consultation with the Axis developer, the IQR is normalised to the median burst duration.
    end
end
