% FUNCTION [ THE_FAR_FRR, EER, X_N, X_EER ] = FAR_FRR_CURVE( X, GT, NPOINTS, VERBOSE );
%
% For each of NPOINTS threshold we compute FAR and FRR using the FAR_FRR function.
% Result is stored in the corresponding column of THE_FAR_FRR.
% Threshold are defined from the data X.
%
% X:  vector of values. May contain Inf or NaN values. NaN values are ignored.
% GT: vector of binary values. 1 means "to be detected". 0 means "not to be detected".
%     GT( k ) indicates whether X( k ) should be above or below a detection threshold.
% NPOINTS: integer. Number of thresholds = number of points in the FAR/FRR curve.
% VERBOSE: binary value. Default 0.
%
% THE_FAR_FRR: 2 x NPOINTS matrix of values in the [ 0 1 ] interval.
%              Row 1 contains FAR values.
%              Row 2 contains FRR values.
% EER: scalar value. Equal Error Rate point estimated from THE_FAR_FRR.
% X_N: vector of NPOINTS values. List of the thresholds.
%
% By Guillaume Lathoud 2004 - lathoud@idiap.ch

function [ the_far_frr, eer, x_n, x_eer ] = far_frr_curve( x, gt, NPOINTS, verbose );
  
  if nargin < 3
    error( 'far_frr_curve needs at least 3 input arguments!' );
  end

  if ~exist( 'verbose', 'var' )
    verbose = 0;
  end
  
  if sum( size( x ) > 1 ) > 1
    error( 'far_frr_curve needs a vector for input argument "x" !' );
  end
  
  if sum( size( gt ) > 1 ) > 1
    error( 'far_frr_curve needs a vector for input argument "gt" !' );
  end
  
  if length( x ) ~= length( gt )
    error( 'far_frr_curve: "x" and "gt" must have same length!' );
  end
  
  
  
  if NPOINTS < 3
    error( 'far_frr_curve needs NPOINTS >= 3 !' );
  end

  usable_data = find( ~isnan( x ) );
  if isempty( usable_data )
    error( 'far_frr_curve could not find any usable data!' );
  end

  % We use the data itself to define the list of thresholds

  x_us = sort( x( usable_data ) );
  n_us = length( x_us );
  k = 1 + round( ( 0:NPOINTS-1 ) / ( NPOINTS - 1 ) * ( n_us - 1 ) );
  x_n = x_us( k );
  
  % Perturbate a bit so that thresholds don't fall exactly on data points
  
  k = round( NPOINTS / 2 );
  
  x_n( 1:k ) = x_n( 1:k ) - 1e-10;
  x_n( k+1:end ) = x_n( k+1:end ) + 1e-10;
  
  % Loop through the list of thresholds
  
  the_far_frr = repmat( NaN, 2, NPOINTS );
  
  if verbose
    chrono = chrono_start( NPOINTS );
  end
  
  for a = 1:NPOINTS
    
    threshold = x_n( a );
    
    result = x > threshold;
    
    tp = sum( result & gt );
    fp = sum( result & (~gt) );
    tn = sum( (~result) & (~gt) );
    fn = sum( (~result) & gt );
    
    [ far, frr ] = far_frr( tp, fp, tn, fn );
    
    the_far_frr( 1, a ) = far;
    the_far_frr( 2, a ) = frr;
    
    if verbose
      chrono = chrono_check( chrono, a );
    end
    
  end
  
  if verbose
    chrono_stop( chrono );
  end
  
  % Find equal error rate (EER)
  
  tmp = diff( the_far_frr( 1:2, : ), [], 1 );
  ind = find( abs( diff( sign( tmp ) ) ) ~=0 );
  
  if isempty( ind )
    error( 'far_frr_curve could not define EER!' );
  end
  
  if length( ind ) > 1
    if ~strcmp( warning, 'off' )
      disp( 'far_frr_curve: WARNING! Multiple definitions for EER!' );
    end
  end
  
  tmp = the_far_frr( 1:2, ind:ind+1 );
  eer = mean( tmp(:), 1 );
  x_eer = mean( x_n( ind:ind+1 ) );