function static_sigma( classify_filename, output_dir, alpha, a_floor )

% Difference with static_bic.m:  instead of maximizing BIC,
% we minimize sum( [ cluster_list.sigma ] .^ alpha )
%
% (criterion already used in Ajmera, Lathoud & McCowan, ICASSP 2004).

  if nargin < 2
    error( [ mfilename ' needs at least 2 input arguments.' ] );
  end

  if ~exist( 'alpha', 'var' )
    % Power defining the sigma criterion
    alpha = 2;
  end

  if ~exist( 'a_floor', 'var' )
    % Floor on the standard deviation, mainly for 1-element clusters
    a_floor = 1e-4;
  end  

  % For a speedup on large meetings

  DO_FAST_NBEST = 20;
  
  % ( 0 ) Check for input/output

  % Check the input file

  disp( ' ' );
  
  if ~exist( classify_filename, 'file' )

    disp( [ mfilename ': input filename does not exists. Skipped.' ] );
    disp( [ '"' classify_filename '"' ] );
    return;

  end

  d = dir( classify_filename );

  if ( numel( d ) ~= 1 )
    
    disp( [ mfilename ': cannot determine a single input filename. Please avoid "*" in the filename.' ] );
    disp( [ mfilename ': skipped.' ] );
    disp( [ '"' classify_filename '"' ] );
    return;

  end

  if ( d.bytes < 2 )

    disp( [ mfilename ': input file contains less than 2 bytes. Skipped.' ] );
    disp( [ '"' classify_filename '"' ] );
    return;
    
  end

  % Ok let's process the input file
  
  [ p, n ] = fileparts( classify_filename );

  output_filename = fullfile( output_dir, [ mfilename '-' n '.mat' ] );

  if exist( output_filename, 'file' )
    disp( [ mfilename ': output filename already exists. Skipped.' ] );
    return;
  end

  % Lock the output file by creating it
  
  fid = fopen( output_filename, 'w' );
  fclose( fid );

  % Load the input file
  
  load( classify_filename, 'result' );

  disp( [ mfilename ' loaded "' classify_filename '".' ] );
  
  % Process several recordings in the list
  % ( I don't think this is used --> only one usually )

  for a = 1:numel( result.list )
    
    % ( 1 ) Fit a Gaussian on each short-term cluster, in azimuth space
    
    tag_list = unique( result.list( a ).sphloc_mat( 7, : ) );
    
    ncluster = numel( tag_list );

    clear cluster_list

    for b = 1:ncluster

      az = result.list( a ).sphloc_mat( 3, find( result.list( a ).sphloc_mat( 7, : ) == tag_list( b ) ) );

      % Avoid border effects in the modulo 2*pi space

      a_center = angle( mean( exp( j * az ) ) );
      az       = -pi + a_center + mod( pi + az - a_center, 2 * pi );
      
      % Evaluate the Gaussian parameters
      
      mu    = mean( az );
      sigma = std( az );
      n     = numel( az );

      % Flooring
      
      if sigma < a_floor
	disp( sprintf( [ mfilename ': warning! Flooring standard deviation for the cluster with tag #%d' ], tag_list( b ) ) );
	sigma = a_floor;
      end
      
      % Store
      
      cluster_list( b ).tag   = tag_list( b );
      cluster_list( b ).mu    = mu;
      cluster_list( b ).sigma = sigma;
      cluster_list( b ).n     = n;  % number of samples
      cluster_list( b ).log_lkld = log_lkld_single_gaussian( mu, sigma, n );
      cluster_list( b ).K     = 2;  % number of free parameters

    end
    
    % ( 2 ) Order the cluster in azimuth space

    [ tmp, ind ] = sort( [ cluster_list.mu ] );

    cluster_list = cluster_list( ind );

    % ( 3 ) Start iterative merging of the clusters.
    % We only try pairs of consecutive clusters in azimuth space,
    % in order to make it faster (square cost instead of cubic cost).

    disp( sprintf( [ mfilename ': starting iterative merging (sigma_criterion with alpha=%5.3f, starting with %d clusters)...' ], alpha, numel( cluster_list ) ) );

    % For a speedup on large meetings
    % - do nbest search as long as possible.
    % - then finish with full search.
    do_fast = 1;

    % Global criterion

    sigma_criterion = sum( [ cluster_list.sigma ] .^ alpha );

    % Verbosity

    chrono = chrono_start( numel( cluster_list ), 30 );

    while 1

      ncluster = numel( cluster_list );

      % ( 3.1 ) Try pairs of consecutive clusters in azimuth space

      clear merge_list

      merge_list = [];

      % Problem: in a long meeting (1 hour) there are many clusters
      %
      % Solution: 
      % - Start investigating ONLY pairs of consecutive clusters that are THE CLOSEST (do_fast = 1).
      % - Do this as long as possible.
      % - When it is not possible anymore, try ALL pairs of consecutive clusters (do_fast = 0).
      %
      % We mark this with a binary-valued flag "do_fast" (initially 1).

      if do_fast 

	% DO_FAST_NBEST pairs of consecutive clusters
	% ( the closest ones )
	
	tmp = [ cluster_list.mu ];

	delta_az = -pi + mod( pi + tmp - tmp( [ 2:end 1 ] ), 2 * pi );

	[ tmp, cluster_ind ] = sort( abs( delta_az ) );

	cluster_ind = cluster_ind( 1:min( numel( cluster_ind ), DO_FAST_NBEST ) );
	
      else

	% ALL pairs of consecutive clusters

	cluster_ind = 1:ncluster;

      end
      

      for b = cluster_ind(:).'
	
	u = b;
	v = 1 + mod( b, ncluster );
	
	mu_u = cluster_list( u ).mu;
	mu_v = cluster_list( v ).mu;

	a_center = angle( mean( exp( j * [ mu_u mu_v ] ) ) );

	mu_u = -pi + a_center + mod( pi - a_center + mu_u, 2 * pi );
	mu_v = -pi + a_center + mod( pi - a_center + mu_v, 2 * pi );
	
	[ mu, sigma, n ] = merge_single_gaussian( mu_u, cluster_list( u ).sigma, cluster_list( u ).n, mu_v, cluster_list( v ).sigma, cluster_list( v ).n );

	log_lkld = log_lkld_single_gaussian( mu, sigma, n );
	
	K = 2;
	
	% Store

	new_merge.u     = u;
	new_merge.v     = v;
	new_merge.mu    = mu;
	new_merge.sigma = sigma;
	new_merge.n     = n;
	new_merge.log_lkld = log_lkld;
	new_merge.K     = K;

	% Estimate what the global score would become if we apply this merge

	new_merge.sigma_criterion  = sigma_criterion - sum( [ cluster_list( [ u v ] ).sigma ] .^ alpha ) + sigma .^ alpha;

	% Store the new merge

	merge_list = cat( 1, merge_list, new_merge );

      end
      
      % ( 3.2 ) Pick the best choice

      [ tmp_sigma_criterion, min_ind ] = min( [ merge_list.sigma_criterion ] );

      % ( 3.3 ) Termination test

      if tmp_sigma_criterion > sigma_criterion

	if do_fast
	  
	  % DO_FAST_NBEST search finished,
	  % let's continue with full search.

	  do_fast = 0;

	  disp( [ mfilename ': finishing with full search.' ] );

	else

	  % Both DO_FAST_NBEST and FULL search are finished:
	  % stop, because the global sigma criterion did not diminish.
	  
	  break

	end
	
      else

	% ( 3.4 ) We did get an improvement (reduction of sigma_criterion),
	% so let's apply the best merge.
	
	u = merge_list( min_ind ).u;
	v = merge_list( min_ind ).v;
	
	cluster_list( u ).tag      = [ cluster_list( [ u v ] ).tag ];
	cluster_list( u ).mu       = merge_list( min_ind ).mu;
	cluster_list( u ).sigma    = merge_list( min_ind ).sigma;
	cluster_list( u ).n        = merge_list( min_ind ).n;
	cluster_list( u ).log_lkld = merge_list( min_ind ).log_lkld;
	cluster_list( u ).K        = merge_list( min_ind ).K;
	
	cluster_list = cluster_list( [ 1:v-1 v+1:end ] );
	
	sigma_criterion = merge_list( min_ind ).sigma_criterion;
	
	% Sanity check
	if abs( sigma_criterion - sum( [ cluster_list.sigma ] .^ alpha ) ) > 1e-8 * abs( sigma_criterion )
	  error( [ mfilename ' is insane!' ] );
	end
	
      end
      
      % ( 3.5 ) Termination test

      if numel( cluster_list ) < 2

	break;

      end

      % Verbosity

      %disp( sigma_criterion)

      chrono = chrono_check( chrono, chrono.n - numel( cluster_list ) );
      
    end  % Iterative pairwise merging, trying to minimize sum( [ cluster_list.sigma ] .^ alpha );
    
    % Verbosity

    %disp( tmp_sigma_criterion )

    chrono.n = chrono.n - numel( cluster_list );

    chrono_stop( chrono );

    % Store the result
    
    static_sigma.list( a ).cluster_list = cluster_list;
    static_sigma.list( a ).name         = result.list( a ).name;
    
  end % Loop through recordings
  

  % ( 4 ) Write out

  static_sigma.classify_filename = classify_filename;
  static_sigma.alpha             = alpha;

  save( output_filename, 'static_sigma' );
  
  disp( [ mfilename ' wrote "' output_filename '".' ] );
  
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analytical merge of 2 single-Gaussian processes
% mu: mean
% sigma: std dev
% n: number of samples

function [mu12, sigma12, n12, diag_sigma12] = merge_single_gaussian( mu1, sigma1, n1, mu2, sigma2, n2, unbiased )
  
  if nargin < 6
    error( [ mfilename '.merge_single_gaussian() needs at least 6 input arguments.' ] );
  end
  
  if ~exist( 'unbiased', 'var' )
    unbiased = 1;
  end

  mu12 = ( n1 * mu1 + n2 * mu2 ) / ( n1 + n2 );
  
  n12 = n1 + n2;

  if unbiased
    
    sigma12 = ( ( n1 - 1 ) * sigma1 + n1 * mu1(:) * mu1(:).' + ( n2 - 1 ) * sigma2 + n2 * mu2(:) * mu2(:).' - ( n1 + n2 ) * mu12(:) * mu12(:).' ) / ( n1 + n2 - 1 );
      
  else

    sigma12 = ( n1 * sigma1 + n1 * mu1(:) * mu1(:).' + n2 * sigma2 + n2 * mu2(:) * mu2(:).' - ( n1 + n2 ) * mu12(:) * mu12(:).' ) / ( n1 + n2 );
    
  end
  
  % Make sure that it is positive definite ( for the inversion in log_mvnpdf.m )
  [ R, p ] = chol( sigma12 );

  diag_sigma12 = ( p~= 0 );

  if diag_sigma12
    % If not positive definite, fall back on the diagonal covariance matrix
    sigma12 = diag( diag( sigma12 ) );
  end

  mu12     = full( mu12 );
  sigma12  = full( sigma12 );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analytical expression of the log likelihood
% of the data on which a single Gaussian is fitted
%
% mu: mean
% sigma: std dev
% n is the number of samples

function log_lkld = log_lkld_single_gaussian( mu, sigma, n, unbiased )

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

  if ~exist( 'unbiased', 'var' )
    unbiased = 1;
  end

  d = numel( mu );

  log_lkld = - n * d / 2 * log( 2 * pi ) - n / 2 * log( abs( det( sigma ) ) );

  if unbiased
    log_lkld = log_lkld - ( n - 1 ) * d / 2;
  else
    log_lkld = log_lkld - n * d / 2;
  end
  

