function eval_detloc( classify_cluster_filename )
  
  if nargin < 1
    error( [ mfilename ' needs 1 input argument.' ] );
  end
  
  % Some parameters
  PAR = [];

  % Where we get the 3D position of the mics,
  % relative to the cameras.
  % ( see documentation about AV16.3 corpus )
  PAR.GTMICARRAYPOS = 'DATA/calibration.mat';

  % Where we get the 3D position of the people,
  % relative to the cameras.
  % ( ibid. )
  % '%s' is to be replaced with the recording name.
  PAR.GTDIR_TEMPLATE = 'DATA/%s';

  % Position of the audio time 0
  % in terms of video timestamp
  %
  % In other words audio starts at video timestamp "PAR.AUDIO_ZERO_VIDEOTIMESTAMP"
  PAR.AUDIO_ZERO_VIDEOTIMESTAMP = 10;

  disp( ' ' );

  % Check for the output file

  [ p, n ] = fileparts( classify_cluster_filename );

  out_filename = fullfile( p, [ mfilename '-' n '.mat' ] );

  % Has this output file been created already?

  if exist( out_filename )

    % Yes -> skip it
    
    disp( [ mfilename ': output file already exists. Skipped.' ] );
    disp( [ '("' out_filename '")' ] );

    return;

  end

  % No -> create it (to lock it from other processes)

  fclose( fopen( out_filename, 'wt' ) );
 
  % Load result file (output of cluster-level speech/non-speech classification)

  global result  % it may be big -> pass it to subfunctions through global access
  load( classify_cluster_filename, 'result' );
  disp( [ mfilename ' loaded "' classify_cluster_filename '".' ] );
 
  % Process each recording in the result

  n_rec = numel( result.list );
  
  for i_rec = 1:n_rec 
    
    % ( 1 ) Load corresponding ground-truth file

    name  = result.list( i_rec ).name;
    gt_dir          = sprintf( PAR.GTDIR_TEMPLATE, name );

    if ~exist( gt_dir, 'dir' ) | ~exist( PAR.GTMICARRAYPOS, 'file' )

      disp( [ mfilename ': not all GT available for recording "' name '". Skipped' ] );

    else
    
      gt = load_gt( gt_dir, name, PAR.GTMICARRAYPOS, PAR.AUDIO_ZERO_VIDEOTIMESTAMP );

      % ( 2 ) Evaluate both detection and localization,
      % considering only frames where all persons have a ground-truth location.
      
      restricted = one_eval( i_rec, gt, 1);

      % ( 3 ) Evaluate considering all frames.

      allframes = one_eval( i_rec, gt, 0 );

      % ( 4 ) Store the evaluation result
      
      eval_result.list( i_rec ).name = name;
      eval_result.list( i_rec ).gt_dir         = gt_dir;
      eval_result.list( i_rec ).gt             = gt;
      eval_result.list( i_rec ).classify_cluster_filename = classify_cluster_filename;
      eval_result.list( i_rec ).PAR            = PAR;
      eval_result.list( i_rec ).restricted     = restricted;
      eval_result.list( i_rec ).allframes      = allframes;
      
    end

  end

  % ( 5 ) Write out
  
  save( out_filename, 'eval_result' );
  disp( [ mfilename ' wrote out "' out_filename '".' ] );

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load the video ground-truth 3D location of each person,
% and interpolate it to match the audio time scale.

function gt = load_gt( gt_dir, name, micarrayposfile, audio_zero_videotimestamp )

  if nargin < 4
    error( [ mfilename '.load_gt() needs 4 input argument.' ] );
  end
  
  % We will need some information in order to do the time interpolation  
  
  global result 
  
  % Load the position of the microphone array

  load( micarrayposfile, 'gt' );

  % Position of array #1

  gt.array1_pos = mean( gt.mic_pos( :, 1:8 ), 2 );
    
  % Append the parameters given to this function
  
  gt.gt_dir = gt_dir;
  gt.name = name;
  gt.micarrayposfile = micarrayposfile;
  gt.audio_zero_videotimestamp = audio_zero_videotimestamp;

  % Append the motion GT of each person
  % that was derived from video
  % ( no speech/silence GT )

  a_template = fullfile( gt_dir, [ name '-person%d-interpolated.3dmouthgt' ] );

  a = 0;

  while 1

    % Move to the next speaker (if any)
    a = a + 1;

    a_filename = sprintf( a_template, a );

    if ~exist( a_filename, 'file' ),  break, end;

    % Load the matrix
    % From http://mmm.idiap.ch/Lathoud/av16.3_v6/FORMATS :
    %   col 1: video timestamp
    %   col 2: binary, 0 if interpolated GT, 1 if manual GT
    %   col 3: x
    %   col 4: y
    %   col 5: z
    %
    % We transpose it.

    gt.person( a ).video_cart_mat = load( a_filename, '-ASCII' ).';
    
    % Convert into a list of audio timestamps
    % and azimuth w.r.t. mic. array #1
    
    audio_time = gt.person( a ).video_cart_mat( 1,: ) - audio_zero_videotimestamp;
    is_manual  = gt.person( a ).video_cart_mat( 2,: );
    
    xyz        = gt.person( a ).video_cart_mat( 3:5,: );
    xyz        = xyz - repmat( gt.array1_pos(:), 1, size( xyz, 2 ) );
    
    sph        = cart2sph_mat( xyz );
    
    % Store the audio information
    % "audio_time" is expressed in terms of seconds (non-integer value).
    
    gt.person( a ).audio_array1_sph_mat = [ audio_time; is_manual; sph ];    
  
    % Interpolate it so that we can directly compare it with the results
    % The new "frame_ind" time is expressed in terms of audio frames
    % (integer value).

    frame_ind = result.list.start_frame:result.list.stop_frame;
    
    % Define the tolerance for linear interpolation, in numbers of audio frames
    % ( +1% to tolerate numerical imprecision errors )
    
    tmp = diff( sort( audio_time ) );
    vgtfd_sec = min( tmp( find( tmp > 0  ) ) );
    if isempty( vgtfd_sec )
      error( [ mfilename ' could not determine video ground-truth frame duration.' ] );
    end
    
    tolerance = 1.01 * vgtfd_sec / result.list.frameshift_sec;
    
    % Get the original video GT timestamps in terms of audio frame index
    % -> non-integer values.
    
    v_frame_ind = 1 + ( gt.person( a ).audio_array1_sph_mat( 1,: ) - result.list.framelength_sec / 2 ) / result.list.frameshift_sec;
    
    % First we need to find the times on which to interpolate it
    
    d_frame_ind = repmat( NaN, size( frame_ind ) );
    
    for b = 1:numel( frame_ind )
      
      % We can only interpolate BETWEEN two CONSECUTIVE video frames
      % -> measure the time distance to the closest two video GTs.
      
      [ d1, ind ] = min( abs( frame_ind( b ) - v_frame_ind ) );  % previous
      d2 = min( abs( frame_ind( b ) - v_frame_ind( [ 1:ind-1 ind+1:end ] ) ) );

      % BETWEEN two CONSECUTIVE (...) iff MAX( d1, d2 ) < 1 video frame
      d_frame_ind( b ) = max( d1, d2 );
      
    end
    
    interp_frame_ind = frame_ind( find( d_frame_ind < tolerance ) );
    
    % Now that we now what we want, let's do linear interpolation
    
    % timestamp in terms of frames
    
    interp_sph = interp_frame_ind(:).';
    
    % is_manual
    
    interp_sph( 2,: ) = 0;
    
    % azimuth
    
    old_az = gt.person( a ).audio_array1_sph_mat( 3,: );
    
    interp_sph( 3,: ) = interp1( v_frame_ind, unwrap( old_az ), interp_frame_ind, 'linear', 'extrap' );
    
    % elevation
    
    old_el = gt.person( a ).audio_array1_sph_mat( 4,: );
    
    interp_sph( 4,: ) = interp1( v_frame_ind, unwrap( old_el ), interp_frame_ind, 'linear', 'extrap' );
    
    % radius
    
    old_r  = gt.person( a ).audio_array1_sph_mat( 5,: );
    
    interp_sph( 5,: ) = interp1( v_frame_ind, old_r, interp_frame_ind, 'linear', 'extrap' );
    
    % Store it
    
    gt.person( a ).interp_sph = interp_sph;
    
  end

  % Determine on which audio frames we have GT...

  % ... for all persons
  gt.intersect_frame_ind = gt.person( 1 ).interp_sph( 1, : );

  % ... for at least one person
  gt.union_frame_ind     = gt.person( 1 ).interp_sph( 1, : );
  
  
  for a = 2:numel( gt.person )

    gt.intersect_frame_ind = intersect( gt.intersect_frame_ind, gt.person( a ).interp_sph( 1, : ) );

    gt.union_frame_ind     = union( gt.union_frame_ind, gt.person( a ).interp_sph( 1, : ) );
    
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% One evaluation
% - optionally restrict to frames where all people have a GT.
% - compute the azimuth errors.
% - fit a Gaussian + Uniform model on azimuth errors.
% - derive detection statistics.
% - return it all!

function one_eval_result = one_eval( i_rec, gt, IS_RESTRICTED )

  if nargin < 3
    error( [ mfilename '.one_eval() needs 3 input arguments.' ] );
  end

  % To access the result

  global result 
  
  % Note: for some recordings (e.g. "seq37"), we KNOW that the ground-truth is always available.

  if IS_RESTRICTED & ~ismember( { result.list( i_rec ).name }, { 'seq37-3p-0001' } );

    % - optionally restrict to frames where all people have a GT.
    which_loc_est = find( ismember( result.list( i_rec ).sphloc_mat( 1, : ), gt.intersect_frame_ind ) );
    
  else

    % - or evaluate on all frames where there is at least one GT.
    which_loc_est =  find( ismember( result.list( i_rec ).sphloc_mat( 1, : ), gt.union_frame_ind ) );

  end
  
  % - compute the azimuth errors 

  % audio frame index
  frame_ind    = result.list( i_rec ).sphloc_mat( 1, which_loc_est );
  
  % signed values in radians
  az_error     = repmat( NaN, 1, numel( which_loc_est ) );

  % unsigned values in radians
  abs_az_error = repmat( +Inf, 1, numel( which_loc_est ) );
  
  for a = 1:numel( frame_ind )

    f = frame_ind( a );

    for b = 1:numel( gt.person )

      c = find( gt.person( b ).interp_sph( 1, : ) == f );

      if ~isempty( c )
	
	new_az_error = -pi + mod( pi + result.list( i_rec ).sphloc_mat( 3, which_loc_est( a ) ) - gt.person( b ).interp_sph( 3, c ), 2 * pi );
	new_abs_az_error = abs( new_az_error );

	if new_abs_az_error < abs_az_error( a )

	  az_error( a ) = new_az_error;
	  abs_az_error( a ) = new_abs_az_error;

	end

      end

    end

  end

  % Sanity check
  if any( isnan( [ az_error abs_az_error ] ) ) | any( isinf( [ az_error abs_az_error ] ) )
    error( [ mfilename '.one_eval() found an undesirable NaN or Inf!' ] );
  end
  
  % - fit a Gaussian + Uniform model on azimuth errors.
  
  gu = em_train_gu( az_error, -pi, +pi );
  
  % - derive detection statistics.
  
  log_post = compute_post_gu( gu, az_error(:) );

  unique_frame_ind = unique( frame_ind );

  frame_n_correct   = zeros( size( unique_frame_ind ) );
  frame_n_incorrect = zeros( size( unique_frame_ind ) );

  for a = 1:numel( unique_frame_ind )
    
    ind = find( frame_ind == unique_frame_ind( a ) );
    
    frame_n_correct( a )   = exp( my_logsum_fast( log_post( ind, 1, 2 ) ) );

    frame_n_incorrect( a ) = exp( my_logsum_fast( log_post( ind, 1, 1 ) ) );
    
  end

  % "frame_n_correct" is most often an almost discrete variable
  % ( number of sources correctly found )
  % -> discretize it to read in a compact manner its histogram
  
  discrete_correct.n_max = ceil( max( frame_n_correct ) );
  
  discrete_correct.hcx   = 0:discrete_correct.n_max;
  
  discrete_correct.edges = 0.5 + ( [ -1  discrete_correct.hcx ] );

  discrete_correct.hcw   = histc( frame_n_correct, discrete_correct.edges );
  discrete_correct.hcw   = discrete_correct.hcw( 1:end-1 );

  %


  if ~IS_RESTRICTED

    % if we also evaluate on frames where not all true locations are known,
    % then the number of incorrect measurements loses its meaning.

    discrete_correct.hcw( 1 ) = 0;
    
    discrete_correct.hcw_norm = discrete_correct.hcw / sum( discrete_correct.hcw );
    
    discrete_incorrect = [];



  else
    
    % Case where we evaluate only on frames where ALL true location are known.
    % In such a case, the number of incorrect measurements makes sense.

    discrete_correct.hcw_norm = discrete_correct.hcw / sum( discrete_correct.hcw );
    
    %
    
    discrete_incorrect.n_max = ceil( max( frame_n_incorrect ) );
    
    discrete_incorrect.hcx   = 0:discrete_incorrect.n_max;
    
    discrete_incorrect.edges = 0.5 + ( [ -1 discrete_incorrect.hcx  ] );
    
    discrete_incorrect.hcw   = histc( frame_n_incorrect, discrete_incorrect.edges );
    discrete_incorrect.hcw   = discrete_incorrect.hcw( 1:end-1 );
    
    discrete_incorrect.hcw_norm = discrete_incorrect.hcw / sum( discrete_incorrect.hcw );

  end

  % - return it all!

  one_eval_result.which_loc_est       = which_loc_est;
  one_eval_result.az_error            = az_error;
  one_eval_result.gu                  = gu;
  one_eval_result.unique_frame_ind    = unique_frame_ind;
  one_eval_result.frame_n_correct     = frame_n_correct;
  one_eval_result.frame_n_incorrect   = frame_n_incorrect;
  one_eval_result.discrete_correct    = discrete_correct;
  one_eval_result.discrete_incorrect  = discrete_incorrect;

  % Also copy some  duration information,
  % it may well turn out to be useful
  
  one_eval_result.framelength_sec     = result.list( i_rec ).framelength_sec;
  one_eval_result.frameshift_sec      = result.list( i_rec ).frameshift_sec;
  one_eval_result.duration            = result.list( i_rec ).duration;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Train Gaussian + Uniform on bounded data
%
% ( e.g. -180 <= angle <= 180 )

function gu = em_train_gu( data, u_min, u_max )
  
  if nargin < 3
    error( [ mfilename '.em_train_gu() needs 3 input arguments.' ] );
  end
  
  data = data(:);
  
  % Define the uniform pdf
  
  gu.u_min = u_min;
  gu.u_max = u_max;
  gu.u_pdf = 1 ./ ( u_max - u_min );
  gu.u_prior = 0.5;
  
  % Initialize the Gaussian pdf
  
  gu.g_mu  = mean( data );
  gu.g_sigma = std( data );
  gu.g_prior = 1 - gu.u_prior;
  
  % Run a few iterations of EM
  
  for a = 1:20
    
    % E: compute posteriors
    
    log_post = compute_post_gu( gu, data );
  
    % M: re-evaluate priors
    
    tmp = [  my_logsum_fast( reshape( log_post( :,:,1 ), [], 1 ) );
	     my_logsum_fast( reshape( log_post( :,:,2 ), [], 1 ) ) ];
    
    priors = exp( tmp - my_logsum_fast( tmp ) );
    priors = priors / sum( priors );
    
    gu.u_prior = priors( 1 );
    gu.g_prior = priors( 2 );
    
    
    % M: re-evaluate parameters of the Gaussian
    % ( the uniform pdf is fixed )
    
    w     = exp( log_post( :,:, 2 ) );
    w     = w(:) / sum( w(:) );
    
    gu.g_mu  = data(:).' * w(:);
    gu.g_sigma = sqrt( ( ( data(:).' - gu.g_mu ) .^ 2 ) * w(:) );
  
  end
  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ log_post, log_joint, log_lkld ] = compute_post_gu( gu, data )
  
  if nargin < 2
    error( [ mfilename '.compute_post_gu() needs 2 input arguments.' ] );
  end
  
  if numel( size( data ) ) > 2
    error( [ mfilename '.compute_post_gu() at most deals with 2-D data.' ] );
  end
  
  log_joint = repmat( -Inf, [ size( data ), 2 ] );
  
  % Joint for the uniform pdf
  
  tmp = log_joint( :,:,1 );
  
  ind = find( ( gu.u_min <= data ) & ( data <= gu.u_max ) );
  tmp( ind ) = log( gu.u_prior ) + log( gu.u_pdf );
  
  log_joint( :,:,1 ) = tmp;
  
  % Joint for the Gaussian pdf
  
  log_joint( :,:,2 ) = log( gu.g_prior ) + log_normpdf( data, gu.g_mu, gu.g_sigma );
  
  % Posteriors
  
  log_lkld = my_logsum_fast( log_joint( :,:,1 ), log_joint( :,:,2 ) );
  
  log_post = log_joint - repmat( log_lkld, [ 1 1 2 ] );
  
