function [pte, ve, ue, gpe, vpe] = pt_eval(gtp, sysp, thr)
% [pte, ve, ue, gpe, vpe] = pt_eval(gtp, sysp)
%    Calculate pitch tracking error.
%    <gtp> is a vector of ground-truth pitch values, in Hz.
%    Values of 0 indicate unvoiced frames; values < 0 indicate
%    missing ground truth (ambiguous frames).
%    <sysp> is a corresponding vector of system output pitch values
%    in Hz.  0 is unvoiced, and negative values are taken as
%    unvoiced too.
%    <thr> is a threshold performance for match.  Default is 0.1 
%    (values within 10% of ground truth are considered a match).
%    Return <ve> as voicing error (proportion of sys frames that
%    are correct where gtruth is > 0), <ue> as unvoicing error 
%    (proportion of sys frames that are correct where gtruth == 0), 
%    <pte> as pitch tracking error (simple average of <ve> and
%    <ue>), and <gpe> as traditional gross pitch error (proportion
%    of frames where both gt and sys report a pitch that are
%    correct). 
%    <vpe> is voicing pitch errors, i.e. the proportion of true
%    voiced frames have a pitch reported, but it's the wrong
%    pitch.  <ve>-<vpe> gives the voicing false reject rate.
% 2013-02-03 Dan Ellis dpwe@ee.columbia.edu

% example usage:
% [ptc, vc, uc, gpc] = pt_eval(addprefixsuffix(idlist, 'ptk/pgt/','-pgt.txt'), addprefixsuffix(idlist, 'ptk/sacc/','-sacc.txt'));


if nargin < 3;  thr = 0.1;  end

if iscell(sysp)

  % passed a cell array of sysps, score them all
  ncell = length(sysp);
  pte = zeros(1, ncell);
  ve = zeros(1, ncell);
  ue = zeros(1, ncell);
  gpe = zeros(1, ncell);
  vpe = zeros(1, ncell);
  for i = 1:ncell
    if iscell(gtp)
      gtp1 = gtp{i};
    else
      gtp1 = gtp;
    end
    sysp1 = sysp{i};
    [pte(i), ve(i), ue(i), gpe(i), vpe(i)] = pt_eval(gtp1, sysp{i}, thr);
  end

else
  
  % regular call for a single evaluation

  % if gtv or sysp are filenames, load them as ptk files
  if ischar(gtp)
    gtpname = gtp;
    [tt,gtp] = pt_read(gtpname);
  end
  if ischar(sysp)
    syspname = sysp;
    [tt,sysp] = pt_read(syspname);
  end
  
  % Now the actual core evaluation of two pitch sequences
  
  gtv = find(gtp > 0);
  gtuv = find(gtp == 0);
  if length(sysp) < length(gtp)
    % pad end of sysp with unvoiced values so it's long enough
    sysp(length(gtp)) = 0;
  end
  sysp = sysp(1:length(gtp));

  % score all frames
  % (works for pitch zero too, but not for -ve frames)
  syscorrect = (sysp >= ((1-thr)*gtp)) & (sysp <= ((1+thr)*gtp));

  ve = 1 - mean(syscorrect(gtv));
  ue = 1 - mean(syscorrect(gtuv));
  gpe = 1 - mean(syscorrect(find((gtp > 0) & (sysp > 0))));
  vpe = 1 - mean(syscorrect(gtv) | (sysp(gtv) == 0));
  
  pte = mean([ve,ue]);

end
