function varargout = ampAtStim_funtest(file, varargin)
%% Function to visually validate the results of ampAtStim
% This function is identical to 'ampAtStim' in most respects, but it 
% additionally creates a plot for every trial showing the data that is used
% for the calculations. These plots can facilitate the verification of the
% results. The code sections that differ from the main function are marked
% as such. Only one file at a time is analysed with ffor all stimuli, but
% the trials can be specified.
% -------------------------------------------------------------------------
% Input     [additional arguments are positional, 
% -----      optional arguments are name-value pairs]
%
%   files:      ID-number(s) of the '.mat' data-file(s), 
%                   Type: integer or vector
%   trials:     number(s) of the trials to be analysed, type: integer or
%   (additional)    vector, default: 'all'
%   Filter:     select between three filter options: 
%   (optional)         'hn'    applys a hum-noise-remove filter
%                      'hnhfn' additonally applys a low pass filter to
%                              remove high frequency noise
%                      If not set, it will be set to "no" (default).
%   StimCell:   set whether to plot the data for a stimulated T-cell ('T')
%   (optional)    or a stimulated interneuron ('INT'). Default: 'T'
%   ErrorBars:  set whether error bars (standard deviation) should be 
%   (optional)    plotted (1 or 0). Default: 1 (yes)
% ------
% Output
% ------
%   visually output of the plots if no output variable is set, new figure 
%   windows are created
%
%   varargout{1}: structure, containing the calculated values, amplitudes 
%                   for the stimuli are sorted from the smallest to the 
%                   biggest stimulus
%                    - MeanAmplitudeT/MeanAmplitudeINT: averaged amplitudes
%                      for the trials where the 'stimCell' was stimulated
%                    - AmplitudeT/AmplitudeINT: amplitudes for the cells
%                      for all trials (not averaged)
%                    - sdT/sdINT: standard deviation for all amplitudes   
%                      over all trials for each stimulus
%                    - stimulatedT: trials with a stimulated T cell
% -------------------------------------------------------------------------
% Program by: Bjarne Schultze [last modified 07.03.2022]
% -------------------------------------------------------------------------

%% Preperations
% parse the user input
[file_num, trials, filter_choice, tol, ~, stimCell, error_bars, ~,~] = ...
    parse_input_peakfun(file, varargin{:});

% Access the data using the function 'extract_data'
dataIn = extract_data(file_num, filter_choice);
recData = dataIn{1,1};
timeVector = recData.timeVector;

%% Find and sort the stimulus peaks
% finds the points, where stimulus is positive
[stim_peaks, stim_locs] = findpeaks(recData.stimulus, 10000);

% find the points, where stimulus is negativ (therefor inverting the
% stimulus data)
[min_stim_peaks, min_stim_locs] = ...
    findpeaks(-recData.stimulus, 10000, 'MinPeakHeight', 0);
% invert the peak values back
min_stim_peaks = -min_stim_peaks;

% concatenating positive and nagative stimlus peak data
all_stim_locs = [min_stim_locs; stim_locs];
all_stim_peaks = [min_stim_peaks; stim_peaks];
% sort the values for the current in ascending order
[all_stim_peaks_sort, sort_order] = sort(all_stim_peaks);
all_stim_locs_sort = all_stim_locs(sort_order);

%% Calculation of the amplitude and its standard deviation
% length of stimulus (here it's the same for all stimuli)
stim_len = 0.5;
stim_num = length(all_stim_locs_sort);                  % number of stimuli

% get the number of trials
if isequal(trials,'all')
    trials = 1:size(recData.recordingT,2);
end
trial_num = length(trials);

% pre-define vectors for the results of the calculations
ampT = zeros(stim_num,trial_num);
ampINT = ampT;

% ------------------------------------------------------------------------
% TEST PART:
if isequal(stimCell, "INT")
    stimulated = setdiff(trials, recData.stimulatedT(:), 'stable');
else
    stimulated = recData.stimulatedT(:);
    if ~isequal(trials,'all')
        stimulated = intersect(stimulated,trials);
    end
end
% ------------------------------------------------------------------------

for exp = trials
    % ----------------------------------------------------
    % TEST-PART: Plot the recording data of the current trial
    if any(ismember(stimulated,exp))
        figure();
        tiledlayout(2,1);
        plt1 = nexttile;
        plot(plt1,recData.timeVector, recData.recordingT(:,exp)); 
        xlabel("time [s]"); ylabel("membrane potential [mV]");
        plt2 = nexttile;
        plot(plt2, recData.timeVector, recData.recordingINT(:,exp));
        xlabel("time [s]"); ylabel("membrane potential [mV]");
    end
    % ----------------------------------------------------
    % iterate through the stimuli calculating amplitude and sd for both cells
    for stim = 1:length(all_stim_locs_sort)
        % select the time frame in which one cell was stimulated
        search_index_amp = timeVector >= all_stim_locs_sort(stim) + tol ...
            & timeVector <= all_stim_locs_sort(stim) + stim_len;
        % select a time frame before the stimulus with the same length as
        % the stimulus as a reference
        search_index_preamp = timeVector <= all_stim_locs_sort(stim) ...
            & timeVector >= all_stim_locs_sort(stim) - stim_len + tol;

        % ----------------------------------------------------
        % TEST-PART: Color the sections used for calculations and add the
        %            corresponding means
        if any(ismember(stimulated,exp))
            % select the stimulus time frame and in the reference time frame
            timeVector_amp = recData.timeVector(search_index_amp);
            timeVector_preamp = recData.timeVector(search_index_preamp);
            % calculate the mean of the recording data for the stimulus time
            % frame and the reference time frame for both cells
            mean_ampT = mean(recData.recordingT(search_index_amp,exp));
            mean_preampT = mean(recData.recordingT(search_index_preamp,exp));
            mean_ampINT = mean(recData.recordingINT(search_index_amp,exp));
            mean_preampINT = mean(recData.recordingINT(search_index_preamp,exp));
            % Color the sections respectively and add the means
            hold(plt1, 'on');
            plot(plt1, timeVector_amp, ...
                recData.recordingT(search_index_amp, exp),'-',...
                'Color', [1,0.38,0.11]);
            plot(plt1, timeVector_preamp,...
                recData.recordingT(search_index_preamp,exp),'-',...
                'Color', [0.2,0.76,0.1]);
            plot(plt1, timeVector_amp([1,end]),[mean_ampT, mean_ampT],'-',...
                'Color',[0.6,0.1,0.1]);
            plot(plt1, timeVector_preamp([1,end]),...
                [mean_preampT, mean_preampT],'-', 'Color', [0.1,0.6,0.1])
            legend(plt1,'recording', 'on stimulus', 'reference',...
                'Location', 'best', 'AutoUpdate', 'off');
            hold(plt1, 'off');
            hold(plt2, 'on');
            plot(plt2, timeVector_amp, ...
                recData.recordingINT(search_index_amp, exp),'-',...
                'Color', [1,0.38,0.11]);
            plot(plt2, timeVector_preamp,...
                recData.recordingINT(search_index_preamp,exp),'-',...
                'Color', [0.2,0.76,0.1]);
            plot(plt2, timeVector_amp([1,end]),[mean_ampINT, mean_ampINT],'-',...
                'Color',[0.6,0.1,0.1]);
            plot(plt2, timeVector_preamp([1,end]), ...
                [mean_preampINT, mean_preampINT],'-', 'Color', [0.1,0.6,0.1])
            legend(plt2,'recording', 'on stimulus', 'reference',...
                'Location', 'best', 'AutoUpdate', 'off');
            hold(plt2, 'off');
        end
        % ----------------------------------------------------
        
        % calculate the amplitude of the membrane potential in the T
        % cell in relation to the amplitude in the reference time frame
        ampT(stim,exp) = ...
            mean(recData.recordingT(search_index_amp,exp)) ...
            - mean(recData.recordingT(search_index_preamp,exp));
        % if the amplitude is smaller than 3 sd in the reference time
        % frame it is considered as zero
        if abs(ampT(stim,exp)) < ...
                3 * std(recData.recordingT...
                (search_index_preamp,exp))
            ampT(stim,exp) = 0;
        end

        % same calculations for the interneuron
        ampINT(stim,exp) = ...
            mean(recData.recordingINT(search_index_amp,exp)) ...
            - mean(recData.recordingINT(search_index_preamp,exp));
        if abs(ampINT(stim,exp)) < ...
                3 * std(recData.recordingINT...
                (search_index_preamp,exp))
            ampINT(stim,exp) = 0;
        end
        % ---------------------------------------------------
        % TEST PART: Add little scales for each amplitude
        if any(ismember(stimulated, exp))
            hold(plt1, 'on'); 
            plot(plt1, timeVector_preamp([end-50,end-50]),...
                [mean_preampT,mean_preampT + ampT(stim,exp)], '-_k'); 
            text(plt1, timeVector_preamp(end-100),...
                ((2 * mean_preampT + ampT(stim,exp))/2), ...
                num2str(round(ampT(stim,exp),1)), 'Color', 'black', ...
                'HorizontalAlignment', 'right','FontSize',8);
            hold(plt1, 'off');
            hold(plt2, 'on'); 
            plot(plt2, timeVector_preamp([end-50,end-50]),...
                [mean_preampINT,mean_preampINT + ampINT(stim,exp)], '-_k'); 
            text(plt2, timeVector_preamp(end-100),...
                ((2 * mean_preampINT + ampINT(stim,exp))/2), ...
                num2str(round(ampINT(stim,exp),1)), 'Color', 'black', ...
                'HorizontalAlignment', 'right', 'FontSize', 8);
            hold(plt2, 'off');
        end
        % ---------------------------------------------------
    end
end

%% Averaging of the amplitudes
% pre-define vectors for the mean and standard deviation of the spiking 
% frequency over all trials
m_ampT = zeros(stim_num,1); m_ampINT = zeros(stim_num,1);
sd_ampT = zeros(stim_num,1); sd_ampINT = zeros(stim_num,1); 

if isequal(stimCell, "INT")
    stimulated = setdiff(trials, recData.stimulatedT(:), 'stable');
    plot_title = sprintf("File: %02.f - Interneuron stimulated", file_num);
else
    stimulated = recData.stimulatedT(:);
    plot_title = sprintf("File: %02.f - T-cell stimulated", file_num);
    if ~isequal(trials,'all')
        stimulated = intersect(stimulated,trials);
    end
end

% find the average and standard deviation of the spiking frequency for each 
% stimulus
for stim = 1:stim_num
    % mean and sd of trials where the T-cell was stimulated
    m_ampT(stim) = mean(ampT(stim,stimulated));
    sd_ampT(stim) = std(ampT(stim,stimulated));
    
    % mean and sd for the interneurons
    m_ampINT(stim) = mean(ampINT(stim,stimulated));
    sd_ampINT(stim) = std(ampINT(stim,stimulated));
end

%% Return and plot results 
% if an output is required, the calculated values are returned (structure)
if isequal(nargout, 1)
    varargout{1} = struct('MeanAmplitudeT',m_ampT, ...
                          'MeanAmplitudeINT', m_ampINT, ...
                          'AmplitudeT', ampT, ...
                          'AmplitudeINT', ampINT, ...
                          'sdT', sd_ampT, 'sdINT', sd_ampINT, ...
                          'stimulatedT', recData.stimulatedT);
else
    % plot the results either with or without error bars
    figure();
    if isequal(error_bars, 1)
        errorbar(all_stim_peaks_sort, m_ampT, sd_ampT);
        axis padded
        hold on;
        errorbar(all_stim_peaks_sort, m_ampINT, sd_ampINT);
        hold off;
        yline(0, 'Color', [0.6,0.6,0.6]);
        title(plot_title);
        xlabel("current [nA]"); ylabel("amplitude [mV]");
        legend("T-cell", "Interneuron", 'Location', 'best');
        subtitle("- error bars representing standard deviation between trials -");
    else
        plot(all_stim_peaks_sort, m_ampT, 'o-',...
            all_stim_peaks_sort, m_ampINT, 'o-', 'MarkerSize', 4);
        axis padded
        title(plot_title);
        xlabel("stimulation current [nA]"); ylabel("amplitude [mV]");
        legend("T-cell", "Interneuron", 'Location', 'best');
        yline(0, 'Color', [0.6,0.6,0.6]);
    end
end

% end of function definition
end