% function [APD,RestPotBeforeStim] = DetermineAPDx(FileName, X)
% Version 2.1, 08.08.2018
% Loads an abf file and determines the APDx
% Input:
%   FileName: .abf File - must contain measured trace and optionally a stimulation trace
%   APDx: APDx value – i.e. 90 for APD90
% optional:
%   'SampleFq' – Sample Frequency (in Hz)
%   'AmpThreshold' – min Threshold of an AP do be detected
%   'StimI' – if there is no Stimulus recorded, the Stimulation time  can be
%             set here (in seconds)
%   'SetWeirdAPtoNaN' - if an AP has its max before stimulated, the
%                       respective output will be set to NaN, set to false to not do that
%
% Output:
%    1dArray containing APDx values for each sweep
%    1dArray containing Resting potential before stimulation
% Example:
%     [APD,RestPotBeforeStim] = DetermineAPDx('MyTrace.abf',90,'StimI', 0.800,'SampleFq', 100, 'AmpThreshold', 0)
%
% Note: 
% This function needs another function to load the respective abffiles. 
% We recommend using the one from https://github.com/fcollman/abfload

function [APD,RestPotBeforeStim] = DetermineAPDx(FileName,APDx,varargin)
parser = inputParser;
DefaultAmpThreshold = 0;
DefaultSampleFq = 10000;
DefaultStimI = [];
DefaultSetWeirdAPtoNan = true;

addRequired(parser, 'FileName', @ischar)
addOptional(parser, 'APDx',APDx, @isnumeric)
addParameter(parser, 'SampleFq', DefaultSampleFq, @isnumeric)
addParameter(parser, 'AmpThreshold',DefaultAmpThreshold, @isnumeric)
addParameter(parser, 'StimI',DefaultStimI, @isnumeric)
addParameter(parser,'SetWeirdAPtoNaN',DefaultSetWeirdAPtoNan, @islogical)

parse(parser,FileName,APDx,varargin{:})

%Load the data
if (exist('abfload') > 0)
    Data = abfload(FileName);
else
    error(['Note: This function needs another function to load the respective abffiles.', ...
           'We recommend using the one from https://github.com/fcollman/abfload']);
end

%Determine Stim point
if isempty(parser.Results.StimI)
    [~, StimI] = max(Data(:,2,:));
    StimI = squeeze(StimI);
else
    StimI =parser.Results.StimI*parser.Results.SampleFq;
end

%Determine avg Potential before Stim
RestPot = mean(Data(1:floor(StimI*0.8),1,:));
StdRestPot = std(Data(1:floor(StimI*0.8),1,:));
[~, MaxI] = max(Data(:,1,:));

WeirdAP = squeeze(MaxI) < 0.98*StimI;
if ((sum(WeirdAP) > 0) && DefaultSetWeirdAPtoNan)
    disp('The following APs have their maximum before being stimulated');
    disp('-> setting respective output to NaN')
    disp(find(WeirdAP == 1));
end

%Find 95% highest Amplitude
for i = 1:size(Data,3)
    try
        [V,~] = sort(Data(Data(:,1,i) > RestPot(:,1,i)+5*StdRestPot(:,1,i),:,i),'ascend');
        MaxinTimeVal(1,1,i) = V((floor(size(V,1)*0.95)));
        %Calc APDx
        VarWinSize = 2*1/(parser.Results.APDx/100);
        [~,APD(i)]=min(abs(Data(StimI+floor(size(V,1)/VarWinSize):end,1,i)-(MaxinTimeVal(1,1,i)-(parser.Results.APDx/100)*abs(MaxinTimeVal(1,1,i)-RestPot(:,1,i)))));
        APD(i) = APD(i)+floor(size(V,1)/VarWinSize);
    catch
        MaxinTimeVal(1,1,i) = -100;
        APD(i) = 0;
    end
end

%Convert to sec and remove all AP < AmpThreshold as well as 'weird' APs
APD(MaxinTimeVal(1,1,:) < parser.Results.AmpThreshold) = 0;
APD=APD'/parser.Results.SampleFq;
RestPotBeforeStim = squeeze(RestPot);

if parser.Results.SetWeirdAPtoNaN
    APD(WeirdAP) = NaN;
    RestPotBeforeStim(WeirdAP) = NaN;
end

end