function dapmat = DAP(SS, varargin)
%DAP Finds dAP(s) on specified channels
%   DAP(SS) finds dAPs on every channel pair specified by
%   SS.analysispars.channelrels by generating a peristimulus histogram,
%   convolving it with a filter (default is a low-pass Gaussian) and 
%   applying an algorithm to determine if there are dAPS in this
%   histogram. The default algorithm is to find all peaks, find the highest
%   valley left or right of this peak, multiplying this valley by 2 and
%   adding .5, and if the peak if greater than this value, then it is
%   considered to be a dAP.
%   
%   DAP(..., cfilter, dap_alg, ploton, varargin) specifies 
%   CFILTER, filter to convolve peristimulus histogram
%   DAP_ALG, function handle to evaluate spike histogram generated by this 
%       function and returns index or indices where a dAP occurs
%   PLOTON, logical to specify whether function will plot convolved
%       peristimulus histogram
%   VARARGIN, either user-defined analysis function parameters or
%       analysispars struct 

    if nargin == 1
        [trange, ~, spktype, resprange, channels] = ...
            SS.ReturnAnalysisPars(SS.analysispars, 1);
        cfilter = normpdf(-3:.2:3);
        dap_alg = @getdap;
        ploton = 0;
    else
        cfilter = varargin{1};
        dap_alg = varargin{2};
        ploton = varargin{3};
        [trange, ~, spktype, resprange, channels] = ...
                SS.ReturnAnalysisPars(varargin(4:end), 1);
    end

    %get ranged spike and stim data
    [spike stim] = SS.ReturnRangedData(trange, spktype);
    
    %calculate bins
    binsize = 1/SS.fs;
    bins = ceil(resprange / binsize);

    %create dapmat, get unique stimulating channels
    dapmat = cell(size(channels, 1), 1);
    stchs = unique(channels(:, 1));
    dmind = 1; %dapmat index

    for x=stchs'
        %for each unique stimulating channel, get stimulations and
        %responding channels
        stims = stim.time(stim.channel == x);
        rechs = channels(channels(:, 1) == x, 2);
        for y=rechs'
            %for each responding channel, get spikes
            spikes = spike.time(spike.channel == y);
            spkhist = zeros(bins, 1);
            for z=1:length(stims)
                %for each stimulation, get responding spikes
                a = find(spikes >= stims(z), 1);
                b = find(spikes < stims(z) + resprange, 1, 'last');
                %convert spks into milliseconds, multiply to get times in
                %range of the size of histogram
                spks = (spikes(a:b)-stims(z))*1000*(bins/(resprange*1000));
                if isempty(spks)
                    continue;
                end
                %create new histogram and add to spike histogram
                spkhist = spkhist + hist(spks, 1:bins)';
            end
            %convolute
            spkhist = conv(spkhist, cfilter, 'same');
            %store results of getdap function into dapmat, increment index
            %since getdap returns index of where dAP(s) occur, must
            %multiply by the binsize to get the actual time of dAP
            dapmat{dmind} = dap_alg(spkhist);
            dapmat{dmind} = [dapmat{dmind}(:, 1) * binsize dapmat{dmind}(:, 2)];
            if ploton
                figure; hold on;
                plot(resprange/length(spkhist):resprange/length(spkhist):resprange, spkhist);
                %next three lines plot a red vertical line where a dAP was
                %detected
%                 for z=1:length(dapmat{dmind})
%                     plot([dapmat{dmind}(z, 1) dapmat{dmind}(z, 1)], [0 spkhist(floor(dapmat{dmind}(z, 1)/binsize))], 'r');
%                 end
                title(['Stim channel: ' num2str(x) ' resp channel: ' num2str(y)])
                xlabel('Time (sec)');
                ylabel('No. of Spikes');
            end
            dmind = dmind + 1;
        end
    end
end

function dap = getdap(spkhist)
%GETDAP Finds dAP(s) according to Bakkum et al. paper
%   DAP = GETDAP(SPKHIST) finds dAP(s) given a spike histogram SPKHIST. The
%   method finds peaks of the histogram, gets the highest valley left and
%   right of the peak, and returns the index of the peak if the peak is
%   2 * highest valley + .5
%   spike histograms can be made with the hist() MATLAB function

    dap = zeros(10, 2);
    dapind = 1;
    [maxima, ind] = findpeaks(spkhist);
    for x=1:length(maxima)
        %for each peak, find left valley (lv) and right valley (rv)
        %left valley
        z=ind(x)-1;
        lv = spkhist(z);
        while z > 1 && spkhist(z-1) < lv
            lv = spkhist(z-1);
            z = z - 1;
        end
        %right valley
        z=ind(x)+1;
        rv = spkhist(z);
        while z < length(spkhist) && spkhist(z+1) < rv
            rv = spkhist(z+1);
            z = z + 1;
        end
        %compare left and right valleys, pick larger value
        if lv > rv
            z = lv;
        else
            z = rv;
        end
        %if peak is 2 * highest valley + .5, then this is a dAP
        if maxima(x) > 2 * z + .5
            dap(dapind, 1) = ind(x);
            dap(dapind, 2) = maxima(x);
            dapind = dapind + 1;
        end
    end
    dap = dap(1:dapind-1, :);
end    