function Performance = ...
    dtPerformance(Filename, detected_tonals, ground_truth_tonals, varargin)
% dtPerformance (Filename, detected_tonals, ground_truth_tonals,...
%                OptionalArgs)
% Calculate performance measure: Recall and Precision
%
% Filename - Name of a audio file
%            Example - 'palmyra092007FS192-071011-230000.wav'
% detected_tonals - tonals detected by graph search algorithm
% ground_truth_tonals - ground truth tonals generated by using GUI tool
%                       dtPlotUIGroundtruth.m
% It is assumed that both sets of tonals are sorted by start time.
% Optional arguments:

%   'Framing', [Advance_ms, Length_ms] - frame advance and length in ms
%       See dtThresh for defaults.
%   'ParameterSet', String or struct
%       Default set of parameters.  May either be a string
%       which is passed to dtThresh or a parameter structure
%       that has been loaded from dtThresh and possibly modified.
%       This argument is processed before any other argument, and other
%       arguments may override these values.
%   'Noise', Method
%       Method for noise compensation in spectrogram plots.
%       It is recommended that the same noise compensation as
%       used for creating the tonal set be used.  See
%       dtSpectrogramNoiseComp for valid methods.
%   'GroundTruthCriteria', [dB, RatioAbove_SNR, MinLen_s] 
%       Criteria for determining whether or not to expect each ground
%       truth tonal to be detected.  When tonals are detected but we
%       did not expect them to be, they are not counted towards the
%       recall, but they will not be used to penalize the precision.
%       Default:  [-inf, .5, 0]
%       Question: Shall we include a minimum coverage parameter?
% 
% Output:
% Performance structure:
%    false_pos - file containing unmatched detections
%    detect_all - structure describing performance irregardless of SNR
%    detect_snr - structure describing performance meeting ground truth
%        criteria.
%    Both detect_all and detect_snr have the following structure:
%        detections - file containing detections matching ground truth
%        gt_idx - index of ground truth to which each detection was matched
%        det_N - # of detections
%        det_matched_N - # detections matching ground truth
%        det_falsepos_N - # detections that do not correspond to gt
%            Note:  det_matched_N + det_falsepos_N is not necessarily det_N
%            in the detect_snr structure as weak ground truth tonals may
%            be detected and we do not penalize them.
%        gt_N - # of ground truth tonals
%        gt_detected_N # of detected ground tonals
%        gt_missed_N - # of undetected ground truth tonals
%        fragment_rate - det_matched_N / gt_detected_N
%        covered_s - seconds of coverage for detected ground truth tonals
%        length_s - length of ground truth tonals that were detected
%    
% Example call:
% [tonals graphs] = dtTonalsTracking(File, t0, t1);
% dtPlotUIGroundtruth(File, tonals, t0, t1);
% dtPerformance(File, tonals, ground_truth_tonal);

stopwatch = tic;

SNR_dB = -inf;  % threshold
RatioAbove_SNR = .333;
NoiseSub = {'median'}; 
%NoiseSub = {'median', [3 3], 3};  % with 3 s MA

%warning off;
import tonals.*;

% defaults

% The threshold set is processed before any other argument as other
% arguments override the parameter set.
thr = dtParseParameterSet(varargin{:});

% Our ground truth annotation tools fit tonals with cubic splines.
% As a result, it is not uncommon for the tonal path to deviate
% slightly from actual tonal.  As a consequence, we search up/down
% for the strongest peak within a specified interval.
PeakTolerance_Hz = thr.PeakTolerance_Hz;

% The adjusted ground truth peak (see thr.PeakTolerance_Hz) must
% be within this distance of the detected frequency to be
% considered a match.
MatchTolerance_Hz = thr.MatchTolerance_Hz;

ms_per_s = 1000;
Length_s = thr.length_ms / ms_per_s;  % frame parameters
Advance_s = thr.advance_ms / ms_per_s;


% Process arguments
error(nargchk(3, inf, nargin));
k = 1;
while k <= length(varargin)
    switch varargin{k}
        case 'Framing'
            if length(varargin{k+1}) ~= 2 || ~ isnumeri(varargin{k+1})
                error('%s must be [Advance_ms, Length_ms]', varargin{k});
            else
                Advance_s = varargin{k+1}(1) / ms_per_s;
                Length_s = varargin{k+1}(2) / ms_per_s;
            end
            k=k+2;
        case 'Noise'
            NoiseSub = varargin{k+1}; k=k+2;
        case 'GroundTruthCriteria'
            if ~ isnumeric(varargin{k+1}) || length(varargin{k+1}) > 3
                error('%s expecting [SNR_dB, RatioAbove_SNR, MinLen_ms]', ...
                    varargin{k});
            else
                params = varargin{k+1};
                if length(params)
                    SNR_dB = params(1);
                    if length(params) > 1
                        RatioAbove_SNR = params(2);
                        if abs(RatioAbove_SNR - .5) > .5
                            error('%s expecting RatioAbove_SNR in [0, 1]', ...
                                varargin{k});
                        end
                        if length(params) > 2
                            % Minimum length s --> ms
                            thr.minlen_ms = params(3)*1000.0;
                        end
                    end
                end
            end
            k=k+2;
        case 'ParameterSet'
            k=k+2;  % already processed
        otherwise
            try
                if isnumeric(varargin{k})
                    errstr = sprintf('Bad option %f', varargin{k});
                else
                    errstr = sprintf('Bad option %s', char(varargin{k}));
                end
            catch
                errstr = sprintf('Bad option in %d''optional argument', k);
            end
            error('%s: %s', mfilename, errstr);
    end
end

header = ioReadWavHeader(Filename);  % Get header information
handle = fopen(Filename, 'rb', 'l'); % open audio data

% end of file time in s
eof_s = sum(header.xhd.byte_length / ...
    (header.nch * header.samp.byte))/header.fs;

% Compute the start time of the last possible start frame given
% the current framing parameters.  Compute the last possible
% start position that will result in a full frame.
lastframe_s = floor((eof_s-Length_s)/Advance_s)*Advance_s;

% Frame parameters in samples
Length_samples = round(header.fs * Length_s);
Advance_samples = round(header.fs * Advance_s);

Nyquist_bin = floor(Length_samples/2);
bin_Hz = header.fs / Length_samples;    % Hz covered by each freq bin

thr.high_cutoff_bins = min(ceil(thr.high_cutoff_Hz/bin_Hz)+1, Nyquist_bin);
thr.low_cutoff_bins = ceil(thr.low_cutoff_Hz / bin_Hz)+1;

% save indices of freq bins that will be processed
range_bins = thr.low_cutoff_bins:thr.high_cutoff_bins;
range_binsN = length(range_bins);  % # freq bin count


Performance.falsePos = java.util.LinkedList();
gtN = ground_truth_tonals.size();
% performance as compared to ground truth without consideration of
% criteria
Performance.all = initPerformance();
% somewhat of a misnomer, this is the performance when we consider
% ground truth tonals that meet a certain criteria which includes
% SNR.
Performance.snr = initPerformance();  
Performance.file = Filename;

block_s = thr.blocklen_s;  % processing block length

% indicator showing which ground truth tonals have been processed
gt_processedI = zeros(gtN, 1);  
% access start/end times for each ground truth tonal
gt_times = zeros(gtN, 2);
for k=1:gtN
    t = ground_truth_tonals.get(k-1).get_time();
    % For search purposes, we consider between the 
    gt_times(k, :) = [t(1), min(t(end), lastframe_s)];
end

falsePos = detected_tonals.clone();
done = gtN == 0;

% Determine what padding is needed.
% padding will be added to each side of the current
% audio data block and removed after the spectrogram
% is computed.
[block_pad_s, block_pad_frames] = ...
    dtSpectrogramNoisePad(Length_s, Advance_s, NoiseSub{:});

% When checking for time frequency nodes that match a ground truth
% tonal, provide a half frame leeway on either side (we don't know
% what the ground truth framing was, so we give the benefit of
% the doubt).
quanterr_s = Advance_s / 2;

while ~ done
    % find first unprocessed ground truth
    gt_idx = 1;  % probably cheaper to search each time instead of 
                 % taking the memcpy hit we would have with deleting rows
                 % my kingdom for a good set of ADTs...
    while ~ done && gt_processedI(gt_idx)
        gt_idx = gt_idx + 1;
        done = gt_idx > gtN;
    end
    if ~ done
        % Set up processing window to span current ground truth
        start = gt_times(gt_idx, 1);
        stop = gt_times(gt_idx, 2);

        % Use the desired block size or the end of this tonal, whichever
        % is longer, making sure that we don't exceed the end of file.
        duration_s = min(eof_s, max(stop, start+block_s)) - start;
        
        % fix stop to reflect selected duration
        stop = start + duration_s;
        
        % Use preferred channel if known
        channel = channelmap(header, Filename);
        
        % Read in the block and compute spectra
        [Signal, spec_dB, t_info, dft, clickP] = ...
            dtProcessBlock( ...
            handle, header, channel, start, block_s, ...
            [Length_samples, Advance_samples], 'Pad', block_pad_s, ...
            'Range', range_bins, ...
            'ClickP', [thr.broadband * range_binsN, thr.click_dB], ...
            'Noise', NoiseSub);
        
        t_idx = t_info.timeidx;
        
        % find all ground truth tonals completely in this range that
        % have not yet been processed
        process = gt_times(:,1) >= start & gt_times(:,2) <= stop & ...
            gt_processedI == false;
        
        
        % for now:
        fprintf('Processing %d gt tonals in %.2f s:  %f to %f:', ...
            sum(process > 0), duration_s, start, stop);
        fprintf(' %d ', find(process));
        fprintf('\n');
        if 0
            for dbg_idx = find(process)'
                fprintf('%4d:  %3.3f - %3.3f\n', dbg_idx, gt_times(dbg_idx, :));
            end
        end
        
        % score unprocessed ground truth tonals in this range
        for gt_idx = find(process)'
            % find time bins corresponding to call
            % Note:  Daniel had +/- quanterr_s padding on groundtruth
            % times, but this causes problems for the mainline code.
            % TODO:  Find other place to include padding on check
            % of matchng tonals
            indices = find(t_idx >= gt_times(gt_idx, 1) & ...
                t_idx <= gt_times(gt_idx, 2))';
            gt_tonal = ground_truth_tonals.get(gt_idx-1);
            if gt_tonal.duration() < Advance_s
                % tonal duration is less than the distance of a single
                % frame.  We cannot process it with these parameters.
                % Either the parameters are incorrect or there was an
                % analyst error.  In any case, we will not try to match
                % anything up to it.
                gt_processedI(gt_idx) = true;  
                continue  % skip
            end

            t = gt_tonal.get_time();
            f = gt_tonal.get_freq();
            [t2, tindices] = unique(t, 'first');
            if length(t2) ~= length(t)
                fprintf('Sort problem - ground truth %d\n', gt_idx);
                % fix it
                t = t2;
                f = f(tindices);
            end
            % interpolate frequencies to current time resolution
            f_i = interp1(t, f, t_idx(indices));
            % quantize frequency
            f_q = round(f_i / bin_Hz) - thr.low_cutoff_bins + 1;
            
            % ground truth was smoothed and peak SNR may not be
            % directly on path. Create matrix to search for max within
            % +/- search_bins on frequency axis
            snr = zeros(1, length(indices));
            search = round(PeakTolerance_Hz / bin_Hz);
            lowbins = f_q - search;
            highbins = f_q + search;
            
            % prevent search from running off spectrogram
            % Overfit splines can exhibit large excursions that can
            % easily exceed the analaysis range.  Hopefully, the 
            % analyst was careful, but this will allows us to process
            % such excursions without an indexing error.
            lowbins(lowbins < 1) = 1;
            lowbins(lowbins > range_binsN) = range_binsN;
            highbins(highbins < 1) = 1;
            highbins(highbins > range_binsN) = range_binsN;
            
            % Find strongest peaks within search range.
            for k=1:length(indices)
                snr(k) = max(spec_dB(lowbins(k):highbins(k), indices(k)));
            end
            sorted_snr = sort(snr);
            gt_snr = sorted_snr(max(1,round((1-RatioAbove_SNR)*length(snr))));
            
            % Are we expecting that the detection algorithm should
            % be able to detect this tonal?
            valid = gt_snr > SNR_dB  && gt_tonal.duration() >= thr.minlen_ms / ms_per_s;
            
            % find overlapping detections
            deviations = [];
            matchedP = false;
            excess_t = 0;
            covered_t = 0;
            % We consider tonals that have already been matched as
            % two whistles may have been reported as a single detection.
            % While the detection will only be counted once with respect
            % to precision, each of the ground truth whistles will be
            % counted as having been detected.
            overlap = gt_tonal.overlapping_tonals(detected_tonals);
            
            % process overlapping detections
            overlap_it = overlap.iterator();
            while (overlap_it.hasNext())
                d = overlap_it.next();
                dt = d.get_time();
                df = d.get_freq();
                % find sections that overlap
                overlap_idcs = find(dt >= t(1) & dt <= t(end));
                t_ov = dt(overlap_idcs);
                f_ov = df(overlap_idcs);
                
                % interpolate ground truth frequencies for this detection
                f_interp = interp1(t,f, t_ov);
                % compute deviations
                deviation = abs(f_ov - f_interp);
                if mean(deviation) <= MatchTolerance_Hz
                    matchedP = true;
                    if falsePos.contains(d)
                        falsePos.remove(d);
                    end
                    Performance.all.detections.add(d);
                    if valid
                        Performance.snr.detections.add(d);
                    end
                    
                    % save statistics
                    % Note:  If multiple detections cover the same
                    % interval, they will both contribute to the coverage
                    % statistic.  This should be addressed in a future
                    % release.  It did not appear to be a major problem
                    % in Silbido (Roch et al. 2011).
                    deviations = [deviations; deviation];
                    excess_t = max(0, t(1)-dt(1)) + max(0, dt(end)-t(end));
                    covered_t = covered_t + t_ov(end) - t_ov(1);
                end
            end
            gt_processedI(gt_idx) = true;
            if matchedP
                Performance.all = updPerformance(...
                    Performance.all, gt_tonal, mean(deviations), ...
                    covered_t, excess_t);
                if valid
                    Performance.snr = updPerformance(...
                        Performance.snr, gt_tonal, mean(deviations), ...
                        covered_t, excess_t);
                end
            else
                Performance.all.gt_miss.add(gt_tonal);
                if valid
                    Performance.snr.gt_miss.add(gt_tonal);
                end
            end
        end
    end
end
fclose(handle);
% Write out information about tonals
Performance.falsePos = falsePos;

sperday = 3600*24;
fprintf('Elapsed time:  %s\n', datestr(toc(stopwatch)/sperday, 13))

function detect = updPerformance(detect, gt_tonal, deviation, ...
    covered_s, excess_s)
% detect = updPerformance(detect, gt_tonal, deviation, covered_s, excess_s)
% update detection statistics in detect structure

detect.gt_match.add(gt_tonal);
detect.deviations(end+1) = deviation;
detect.covered_s(end+1) = covered_s;
detect.excess_s(end+1) = excess_s;
detect.length_s(end+1) = gt_tonal.duration();

function detect = initPerformance()
% detect = initPerformance()
% Given the number of detections and ground truth tonals, initialize
% a structure which tracks statistics on detections.
detect.detections = java.util.LinkedList();
detect.gt_match = java.util.LinkedList();
detect.gt_miss = java.util.LinkedList();
detect.deviations = [];
detect.gt_missed_N = NaN;  % computed at end
detect.fragment_rate = NaN;
detect.covered_s = [];
detect.length_s = [];
detect.excess_s = [];
