function varargout = findIntSpikes(data, varargin)
%% Find interneuron action potentials 
% This function is designed to find interneuron spikes and seperate them
% from noise and other fast changes in the membrane potential. The function
% is capable of finding most of the interneuon spikes with amplitudes down
% to approximately 1 mV. If requested, a plot can be created showing all
% spikes that were detected by the function for validation purposes. 
% -------------------------------------------------------------------------
% Input [optional input arguments are name-value-pairs]
% -----
%   data:           data to be analysed, type: vector
%   PlotFlag:       select if a plot should be created for each trial
%   (optional)          showing the detected peaks. 0 (no plots, default), 
%                       1 (plots)
% ------
% Output
% ------
%   varargout{1}:   peak values [mV]
%   varargout{2}:   peak locations [datapoints]
% -------------------------------------------------------------------------
% Program by:   Bjarne Schultze     [last modified: 01.03.2022]
% -------------------------------------------------------------------------

%% Preparations
% parse the user input
p = inputParser;
addRequired(p,'data', @(x) isvector(x))
addParameter(p,'PlotFlag', [], @(x) isequal(x,1)||isequal(x,0))
parse(p, data, varargin{:})

dataIN = p.Results.data;
plotflag = p.Results.PlotFlag;

% load additional file
load("../AdditionalFiles/standard_stimulus.mat", "standard_stimulus")

%% Main calculations 
[peak_values, locations, ~, prom] = findpeaks(dataIN,...
    'Annotate','extents', 'WidthReference','halfprom',...
    'MinPeakProminence', 1,'MaxPeakWidth',200);

% preassign a vector to gather the quater prominence width
qprom_width = zeros(length(peak_values),1);
% measure the peak width at one quarter of the prominence for each
% spike (default would be half prom)
for spike = 1:length(peak_values)
    % set the quarter prom as reference value for width measurement 
    measure_level = peak_values(spike) - 0.25*prom(spike);
    % select a window of 100 datapoints around the current spike if
    % possible, otherwise a one-sided smaller window
    if locations(spike) < 51
        spikeData = dataIN(1:locations(spike)+50);
        spikeLoc = locations(spike);
    elseif locations(spike)+50 > locations(end)
        spikeData = dataIN(locations(spike)-50:end);
        spikeLoc = 51;
    else
        spikeData = dataIN(locations(spike)-50:locations(spike)+50);
        spikeLoc = 51;
    end
    % devide the points in left and right of the peak
    dataL = spikeData(1:spikeLoc);
    dataR = spikeData(spikeLoc+1:end);
    % take the points (left and right) that are closest to the
    % measure level of a quarter of the prominence
    [~,minL] = min((dataL - measure_level), [], ...
        'ComparisonMethod','abs');
    [~,minR] = min((dataR - measure_level), [], ...
        'ComparisonMethod','abs');
    % calculate the final width at quarter peak prominence
    qprom_width(spike) = minR + spikeLoc - minL;
end
% sort and select the all found peaks in a 3-step algorithm 
% 1st step: exclude every peak beneath the overall mean +1 mV 
select_peaks = peak_values > mean(dataIN)+1;
peaks_subset1 = peak_values(select_peaks);
% select the corresponding values/properties
locs_subset1 = locations(select_peaks);
qprom_width_subset1 = qprom_width(select_peaks);
prom_subset1 = prom(select_peaks);

% 2nd step: select only peaks with a quarter peak prom over 4 and
% under 45
select_peaks = 4 < qprom_width_subset1 & qprom_width_subset1 < 45;
peaks_subset2 = peaks_subset1(select_peaks);
% select the coorresponding values/properties
locs_subset2 = locs_subset1(select_peaks);
prom_subset2 = prom_subset1(select_peaks);

% 3rd step: exclude the peaks smaller than 20% of the largest peak
select_peaks = prom_subset2 > max(prom_subset2)*0.2;
peaks_subset3 = peaks_subset2(select_peaks);
% select the corresponding values/properties
locs_subset3 = locs_subset2(select_peaks);

% using findpeaks again, but without a maximum peak width to find all
% spikes even those with incorrectly measured peak width
[peaks,locs,~,prom] = findpeaks(dataIN,...
        'Annotate','extents', 'WidthReference','halfprom',...
        'MinPeakProminence', 1);
% find onset locations of positive stimuli
[~, stim_locs] = findpeaks(standard_stimulus);
% seperate the one spike per stimulus that was ignored before
for i = 1:length(stim_locs)
    sel_peaks = locs >= stim_locs(i) & locs <= stim_locs(i)+5000;
    [~,maxIndex] = max(prom(sel_peaks));
    if ~isempty(maxIndex)
        locs_subset = locs(sel_peaks); peaks_subset = peaks(sel_peaks);
        extra_locs(i,1) = locs_subset(maxIndex);
        extra_peaks(i,1) = peaks_subset(maxIndex);
    end
end

% collect all spikes and sort them according to their location
[locs_final,sortOrder] = sort([locs_subset3;extra_locs]);
peaks_final = [peaks_subset3;extra_peaks];
peaks_final = peaks_final(sortOrder);

% plot the peaks that are found plus their prominence and width
if plotflag
    figure
    findpeaks(dataIN, 10000,...
        'Annotate','extents', 'WidthReference','halfprom',...
        'MinPeakProminence', 1,'MaxPeakWidth',200);
    locs_plot = locs_final / 10000;
    hold on; plot(locs_plot, peaks_final, 'og'); hold off;
    xlabel("time [s]"); ylabel("membrane potential [mV]");
end

% return peak and location values
if isequal(nargout,1)
    varargout{1} = peaks_final;
elseif isequal(nargout,2)
    varargout{1} = peaks_final;
    varargout{2} = locs_final;
end

% end of function definition
end