function sector_perf = eval_sector_perf( result_mat, gt_mat, verbose )
  
% FUNCTION SECTOR_PERF = EVAL_SECTOR_PERF( RESULT_MAT, GT_MAT, [ VERBOSE ] )
%
% (1) On active frames only: for each active source instance compute TP,FN and recall.
% (2) On all frames: sum over all instances AND over silence -> TP,FP,TN,FN and derived measures.
% 
% See also STATIC_GT_2_SECTOR_GT (includes more explanations).
%
% RESULT_MAT: N_sectors - by - N_frames matrix of zeros and ones.
%
% GT_MAT: N_instances - by - 1 cell array of N_sectors - by - N_frames matrices of integer values (either zero or object id). Non-zero values locate in both space and time the activity of a given active source instance. The value itself gives the identity of the object.
%
% Note: we use a cell array of 2D matrices instead of a 3D matrix 
% because sparse matrices cannot be more than 2-dimensional.
%
% VERBOSE:  optional parameter, 0 or 1 (default value is 1).
%
% "GT" means Ground-Truth
% 
% Note that a given object may be represented by several "active source instances" in the GT,
% e.g. one instance for each GT speech segment of a given speaker.
%
% OUTPUT: a structure containing various evaluation results.
%
% Typical use: some black box search method outputs a binary matrix (e.g. thresholding on a matrix of SAM values). We want to evaluate the result 1) overall 2) for each active source instance. We also want to determine how well multiple concurrent sources can be localized (histograms).

  if nargin < 2
    error( 'eval_sector_perf needs at least two input parameters' );
  end
  
  % ------------------------------------------------------------
  % ( 1 ) Check input parameters
  
  siz1 = size( result_mat );
  if length( siz1 ) ~= 2
    error( 'eval_sector_perf: "result_mat" must be a 2-dimensional matrix' );
  end
  
  if ~iscell( gt_mat )
    error( 'eval_sector_perf: "gt_mat" must be a CELL ARRAY of 2-dimensional matrices.' );
  end

  % How many active source instances ?
  
  N_instances = length( gt_mat );
  
  for a = 1:N_instances
    
    siz2 = size( gt_mat{ a } );
    
    if length( siz2 ) ~= 2
      error( sprintf( 'eval_sector_perf: "gt_mat{%d}" must be a 2-dimensional matrix', a ) );
    end
    
    if ~all( siz1 == siz2 )
      error( sprintf( 'eval_sector_perf: "result_mat" and "gt_mat{%d}" must have same size', a ) );
    end
    
  end
  
  if ~exist( 'verbose', 'var' )
    verbose = 1;
  end
  
  if ~ismember( verbose, [0 1] )
    error( 'eval_sector_perf: "verbose" must be 0 or 1' );
  end
  
  % How many sectors in space ? 
    
  N_sectors = siz2( 1 );
  
  % How many time frames ?
  
  N_frames  = siz2( 2 );
  
  % Reduce memory usage
  
  if ~issparse( result_mat )
    result_mat = sparse( result_mat );
  end
  
  for a = 1:N_instances
    if ~issparse( gt_mat{ a } )
      gt_mat{ a } = sparse( gt_mat{ a } );
    end
  end
    
  % Verbosity
  
  if verbose
    disp( sprintf( 'eval_sector_perf:  N_sectors = %d,  N_frames = %d,  N_instances = %d', N_sectors, N_frames, N_instances ) );
  end
  
  % ------------------------------------------------------------
  % ( 2 ) Store input parameters
  
  sector_perf.result_mat = result_mat;
  
  % We don't store the GT because it is usually quite big
  
  sector_perf.N_sectors   = N_sectors;
  sector_perf.N_frames    = N_frames;
  sector_perf.N_instances = N_instances;
  
  % -------------------------------------------------------------
  % ( 3 ) Evaluate and store results 
  % for each active source instance separately

  for a = 1:sector_perf.N_instances

    % Where we'll store everything about this instance
    an_instance = [];
    
    tmp = gt_mat{a};
    an_instance.object_id = unique( tmp( find( tmp ) ) );

    if length( an_instance.object_id ) < 1
      error( sprintf( 'eval_sector_perf: empty instance #%d', a ) );
    end
    
    if length( an_instance.object_id ) > 1
      error( sprintf( 'eval_sector_perf: multiple object ids for instance #%d', a ) );
    end
    
    % ( 3.1 ) First look at the GT itself
    
    [ ind_sector, ind_frame ] = find( gt_mat{ a } );
     
    % In the GT, list the frames occupied by this active source instance
    an_instance.frame_list  = sort( unique( ind_frame ) );
    an_instance.frame_list  = an_instance.frame_list(:).';
    an_instance.n_frames    = length( an_instance.frame_list );
    
    % In the GT, list the (sector,time frame) occupied by this active source instance
    % -> note that there may be more than one sector per time frame
    %    typically when the true location of an active source 
    %    is at the border between two sectors
    an_instance.sector_list = sub2ind( size( gt_mat{ a } ), ind_sector, ind_frame );
    an_instance.sector_list = an_instance.sector_list(:).';
    
    % ( 3.2 ) Second, compute results
    
    % "result_per_frame": vector of "an_instance.n_frames" binary values
    % 0 means active source not detected
    % 1 means active source detected
    %
    % Note that these 3 lines of code automatically take care of 
    % the borderline cases (equivalent to an OR operator),
    % i.e. when the true location of an active source instance is
    % in more than one sector.
    
    m = reshape( gt_mat{ a }( an_instance.sector_list ) & result_mat( an_instance.sector_list ), length( unique( ind_sector ) ), [] );
    
    an_instance.result_per_frame = any( m, 1 );
    
    % ( 3.3 ) Finally, some statistics
    % true positives, false negatives, and recall
    
    an_instance.TP = sum( an_instance.result_per_frame );
    an_instance.FN = sum( ~an_instance.result_per_frame );
    an_instance.recall = an_instance.TP / ( an_instance.TP + an_instance.FN );
    
    % ( 3.4 ) Store the new object

    sector_perf.instance_list( a ) = an_instance;
    
  end
  
  if verbose
    disp( 'eval_sector_perf: finished evaluating each active source instance' );
  end
  
  % ------------------------------------------------------------
  % ( 4 ) Overall results
  
  sector_perf.TP = sum( [ sector_perf.instance_list.TP ] );
  sector_perf.FN = sum( [ sector_perf.instance_list.FN ] );
  
  % Look at (sector,time frame) pairs NOT occupied by any active source instance
  
  occupied = sparse( N_sectors, N_frames );
  for a = 1:N_instances
    occupied = occupied | gt_mat{ a };
  end
  
  sector_perf.TN = sum( reshape( ~occupied & ~result_mat, 1, [] ) );
  sector_perf.FP = sum( reshape( ~occupied &  result_mat, 1, [] ) );
  
  [ far,  frr,  hter ] = far_frr( sector_perf.TP, sector_perf.FP, sector_perf.TN, sector_perf.FN );
  
  sector_perf.far = full( far );
  sector_perf.frr = full( frr );
  sector_perf.hter = full( hter );

  [ prc, rcl, F ] = prc_rcl( sector_perf.TP, sector_perf.FP, sector_perf.TN, sector_perf.FN );
  
  sector_perf.prc = full( prc );
  sector_perf.rcl = full( rcl );
  sector_perf.F   = full( F );
  
  % Average instance recall
  
  x = [ sector_perf.instance_list.recall ];
  w = [ sector_perf.instance_list.n_frames ];
  sector_perf.mean_instance_recall = sum( x .* w ) / sum( w );
  
  % Minimum instance recall
  
  x = [ sector_perf.instance_list.recall ];
  sector_perf.min_instance_recall = min( x );
  
  
  if verbose
    disp( 'eval_sector_perf: fininshed evaluating overall results' );
  end

  % ------------------------------------------------------------
  % ( 5 ) Results broken down per number of active sources
  %       -> histograms.
  
  % ( 5.1 ) For each frame, compute the number of active sources
  % in the GT and in the result.
  % "sources in the results" means sources correctly localized.
  
  sector_perf.gt_n_active = zeros( 1, N_frames );
  sector_perf.result_n_active = zeros( 1, N_frames );
  
  for a = 1:sector_perf.N_instances
    
    % Look at the GT
    frame_list = sector_perf.instance_list( a ).frame_list;
    sector_perf.gt_n_active( frame_list ) = sector_perf.gt_n_active( frame_list ) + 1;
    
    % Look at the result
    sector_perf.result_n_active( frame_list ) = sector_perf.result_n_active( frame_list ) + sector_perf.instance_list( a ).result_per_frame;
    
  end
  
  % ( 5.2 ) Compute the histogram of number of sources found 
  % for each TRUE number of active sources
  
  sector_perf.gt_n_active_list = sort( setdiff( unique( sector_perf.gt_n_active ), [0] ) );
  
  sector_perf.histo_multisource = cell( max( sector_perf.gt_n_active_list ), 1 );
  
  for a = sector_perf.gt_n_active_list
    
    x = 0:a;
    
    frame_list = find( sector_perf.gt_n_active == a );
    
    n = hist( sector_perf.result_n_active( frame_list ), x );
    
    sector_perf.histo_multisource{ a } = n;
    
  end
  
  if verbose
    disp( 'eval_sector_perf: finished computing histograms.' );
    disp( 'eval_sector_perf: done.' );
  end