function sbmfcc_sns_model( par )

% EM fit of a speech/non-speech model on sector-based MFCCs.
% See IDIAP RR 06-26 for detailed documentation.
% Main idea:
% - diagonal covariance matrices for non-speech.
% - full covariance matrices for speech.
%
% An automatic selection of the size of the model is included.
%
% The final binary decision (speech/non-speech) is taken using
% automatic threshold selection, as in IDIAP RR 05-52.
%
% Note that we first reduce the data to a fixed number of frames,
% keeping approximately the same distribution of C0 coefficients.
%
% By Guillaume Lathoud, 2006.
% lathoud@idiap.ch
  
  
  if nargin < 1
    error( [ mfilename ' needs at least 1 input argument.' ] );
  end
  
  dbstop error
  dbstop warning

  % Check for required parameters

  required_fields = { 'in_filename', 'n_comp_list', 'far_t', 'out_dir' };
  check_param( required_fields, fieldnames( par ) );

  % Define default parameters

  par_default = [];
  
  % To select the size of the GMM, within the "par.n_comp_list"
  par_default.bic_lambda = 1;

  % Before reducing the data, we can sort it in increasing C0 value
  par_default.do_c0_sort = 1;
  
  % To limit the duration of EM fitting, we'll reduce the MFCC data
  % to a fixed number of vectors
  par_default.fit_max_data_size = 10000;
  
  % EM fitting termination parameters
  par_default.em_max_iter = 1000;
  par_default.em_cv_prop  = 1e-8;

  par_default.verbose = 1;
  
  par = fill_default( par, par_default );

  if par.verbose
    disp( ' ' );
    disp( [ mfilename ' parameters:' ] );
    disp( par );
  end

  
  % Load input and process each recording

  load( par.in_filename, 'result' );
  if par.verbose
    disp( [ mfilename ' loaded "' par.in_filename '".' ] );
  end
  n_rec = numel( result.list );
  
  for i_rec = 1:n_rec 
    
    % ( 1 ) Define output filename
    
    [ p, n ] = fileparts( par.in_filename );

    out_filename = fullfile( par.out_dir, [ mfilename '-' n '.mat' ] );
    
    if exist( out_filename, 'file' )

      if par.verbose

	disp( ' ' );
	disp( [ mfilename ' skipped recording ' sprintf( '%d/%d', i_rec, n_rec ) ': "' result.list( i_rec ).name '".' ] );
	disp( [ '( output file already exists: "' out_filename '" )' ] );

      end

    else
      
      % Lock the output file by creating it
      fid = fopen( out_filename, 'w' );
      fclose( fid );
      
      % ( 2 ) Fit the SNS model
      
      if par.verbose
	disp( ' ' );
	disp( [ mfilename ': processing recording ' sprintf( '%d/%d', i_rec, n_rec ) ': "' result.list( i_rec ).name '"...' ] );
      end

      % Extract and reduce sector-based MFCCs
      sbmfcc = result.list( i_rec ).sbmfcc_mat( 3:end, : );
      
      % Reduced version of the data,
      % used for GMM fitting only.
      fit_sbmfcc = sbmfcc;

      N = size( fit_sbmfcc, 2 );
      
      if N > par.fit_max_data_size
	
	if par.do_c0_sort
	  
	  % First, sort the sbMFCC vectors by increasing C0
	  
	  [ tmp, ind ] = sort( fit_sbmfcc( 1,: ) );
	  
	  fit_sbmfcc = fit_sbmfcc( :, ind );
	  
	  if par.verbose
	    disp( ' ' );
	    disp( [ mfilename ' sorted sector-based MFCCs by increasing C0.' ] );
	  end
	  
	end
	
	% Second, pick a reduced number of sbMFCC vectors, at regular intervals

	fit_sbmfcc = fit_sbmfcc( :, 1 + round( ( N - 1 ) * ( 0.5:par.fit_max_data_size ) / par.fit_max_data_size ) );
	
	if par.verbose
	  disp( sprintf( [ mfilename ' reduced sector-based MFCCs from %d frames to %d frames.' ], N, size( fit_sbmfcc, 2 ) ) );
	end
	
	N = size( fit_sbmfcc, 2 );
	
      end
      

      % Fit one GMM for each tentative size in "n_comp_list"
      
      model_list = [];
      for a = 1:numel( par.n_comp_list )
	
	% Full covariance GMM fitting
	
	model_list( a ).n_components    = par.n_comp_list( a );
	
	if par.verbose
	  disp( ' ' );
	  disp( sprintf( 'Fitting a GMM with %d components...', model_list( a ).n_components ) );
	end


	rand( 'state', 206619.17 );   % Same init conditions for all
	
	[ model_list( a ).mvn_model, model_list( a ).BIC ] = train_gmm_bic( fit_sbmfcc, model_list( a ).n_components, 0,  par.verbose, par.bic_lambda, par.em_max_iter, par.em_cv_prop );
	
	K = estimate_K_gmm( model_list( a ).mvn_model );
	
	% Compute misc. stuff related to BIC criterion
	% + (FAR,FRR)-related stuff
	
	if par.verbose
	  disp( [ mfilename ': taking final posterior-based SNS decision, for this model...' ] );
	end
	
	D = size( sbmfcc, 1 );
	N = size( sbmfcc, 2 );
	
	% Compute the posterior of each GMM component
	% Note that we do it on the full data, not just "fit_sbmfcc".
	[ log_post_gmm ] = compute_post_gmm( model_list( a ).mvn_model, sbmfcc );
	
	% Also extract the log-posterior probability of wideband activity
	% (defined in IDIAP RR 05-52).
	log_post_wa = result.list.sphloc_mat( 6,: );
	
	% Also compute the estimated (FAR,FRR) point,
	% for the target FAR_T
	
	ind_speech = find( ~[ model_list( a ).mvn_model.diag_Sigma ] );
	
	if isempty( ind_speech )
	  
	  % No speech model
	  
	  % Store it
	  
	  model_list( a ).par = par;
	  model_list( a ).D = D;
	  model_list( a ).N = N;
	  model_list( a ).M = numel( model_list( a ).mvn_model );
	  model_list( a ).n_diag = sum( [ model_list( a ).mvn_model.diag_Sigma ] );
	  model_list( a ).n_full = sum( ~[ model_list( a ).mvn_model.diag_Sigma ] );
	  model_list( a ).K = K;
	  
	  model_list( a ).log_post_speech = [];

	  model_list( a ).wa_est_far_point = NaN;
	  model_list( a ).wa_est_frr_point = NaN;
	  model_list( a ).wa_threshold     = NaN;
	  
	  model_list( a ).speech_est_far_point = NaN;
	  model_list( a ).speech_est_frr_point = NaN;
	  model_list( a ).speech_threshold     = NaN;
	  
	  model_list( a ).wa_speech_est_far_point = NaN;
	  model_list( a ).wa_speech_est_frr_point = NaN;
	  model_list( a ).wa_speech_threshold     = NaN;
	  
	  model_list( a ).speech_weight = 0;
	  
	  
	else
	  
	  % Usual case: at least one non-diagonal component
	  if numel( ind_speech ) < 2
	    log_post_speech = log_post_gmm( ind_speech,: );
	  else
	    log_post_speech = my_logsum_fast( log_post_gmm( ind_speech, : ) );
	  end

	  % ( 3 ) We will now determine an automatic threshold, as in IDIAP RR 05-52.
	  %
	  % Three types of cues are used:
	  % - wideband activeness alone (from SAM-SPARSE-MEAN, IDIAP RR 05-52).
	  % - speech activeness alone (from sector-based MFCCs, IDIAP RR 06-26).
	  % - both.
	  
	  % ( 3.1 ) wideband alone

	  post = exp( log_post_wa(:).' );

	  [ sorted_post, sorted_ind ] = sort( post );

	  est_frr = cumsum( sorted_post );
	  est_frr = est_frr / sum( sorted_post );

	  est_far = fliplr( cumsum( 1 - fliplr( sorted_post(:).' ) ) );
	  est_far = est_far / sum( 1 - sorted_post );

	  [ tmp, min_ind ] = min( abs( est_far - par.far_t ) );
	  
	  model_list( a ).wa_est_far_point = est_far( min_ind );
	  model_list( a ).wa_est_frr_point = est_frr( min_ind );
	  model_list( a ).wa_threshold     = sorted_post( min_ind );
	  
	  % Store final decision
	  model_list( a ).wa_decision = ( post > model_list( a ).wa_threshold );

	  % ( 3.2 ) MFCC alone
	  
	  post = exp( log_post_speech(:).' );

	  [ sorted_post, sorted_ind ] = sort( post );

	  est_frr = cumsum( sorted_post );
	  est_frr = est_frr / sum( sorted_post );

	  est_far = fliplr( cumsum( 1 - fliplr( sorted_post(:).' ) ) );
	  est_far = est_far / sum( 1 - sorted_post );

	  [ tmp, min_ind ] = min( abs( est_far - par.far_t ) );
	  
	  model_list( a ).speech_est_far_point = est_far( min_ind );
	  model_list( a ).speech_est_frr_point = est_frr( min_ind );
	  model_list( a ).speech_threshold     = sorted_post( min_ind );
	  
	  % Store final decision
	  model_list( a ).speech_decision      = ( post > model_list( a ).speech_threshold );

	  % ( 3.3 ) wideband + mfcc
	  % (independence assumption)

	  post = exp( log_post_wa(:).' + log_post_speech(:).' );
	  
	  [ sorted_post, sorted_ind ] = sort( post );
	  
	  est_frr = cumsum( sorted_post );
	  est_frr = est_frr / sum( sorted_post );
	  
	  est_far = fliplr( cumsum( 1 - fliplr( sorted_post(:).' ) ) );
	  est_far = est_far / sum( 1 - sorted_post );
	  
	  [ tmp, min_ind ] = min( abs( est_far - par.far_t ) );

	  model_list( a ).wa_speech_est_far_point = est_far( min_ind );
	  model_list( a ).wa_speech_est_frr_point = est_frr( min_ind );
	  model_list( a ).wa_speech_threshold     = sorted_post( min_ind );

	  % Store final decision
	  model_list( a ).wa_speech_decision      = ( post > model_list( a ).wa_speech_threshold );

	  % ( 3.4 ) Store misc. stuff
	  
	  model_list( a ).par = par;
	  model_list( a ).D = D;
	  model_list( a ).N = N;
	  model_list( a ).M = numel( model_list( a ).mvn_model );
	  model_list( a ).n_diag = sum( [ model_list( a ).mvn_model.diag_Sigma ] );
	  model_list( a ).n_full = sum( ~[ model_list( a ).mvn_model.diag_Sigma ] );
	  model_list( a ).K = K;

	  model_list( a ).log_post_speech = log_post_speech;
	  
	  model_list( a ).speech_weight = sum( [ model_list( a ).mvn_model( ind_speech ).w ] );
	  
	end % End of fit of full-covariance GMM
	
      end % End of loop through candidate number of GMM components
      
      % Write out 
      
      save( out_filename, 'model_list', 'par' );
      
      if par.verbose
	disp( [ mfilename ' wrote out "' out_filename '".' ] );
      end
      
    end % End of check for output filename
    
  end % End of loop through recordings
  
  if par.verbose
    disp( ' ' );
    disp( [ mfilename ': done.' ] );
  end
  