function static_bic( classify_filename, lambda, a_floor )

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

  if ~exist( 'lambda', 'var' )
    % Default value: follow the original definition of BIC
    lambda = 1;
  end

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

  % ( 0 ) Check for input/output

  [ p, n ] = fileparts( classify_filename );

  output_filename = fullfile( p, [ 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

    % Compute the BIC score for each cluster

    [ cluster_list( 1:ncluster ).BIC ] = deal( [ cluster_list.log_lkld ] - lambda * [ cluster_list.K ] / 2 .* log( [ cluster_list.n ] ) );
    
    % ( 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 (BIC with lambda=%5.3f,  starting with %d clusters)...' ], lambda, numel( cluster_list ) ) );

    % Global BIC score
    MBIC = sum( [ cluster_list.BIC ] );

    while 1

      ncluster = numel( cluster_list );

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

      clear merge_list

      for b = 1:ncluster
	
	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;
	
	BIC = log_lkld - lambda * K / 2 * log( n );

	% Store

	merge_list( b ).u     = u;
	merge_list( b ).v     = v;
	merge_list( b ).mu    = mu;
	merge_list( b ).sigma = sigma;
	merge_list( b ).n     = n;
	merge_list( b ).log_lkld = log_lkld;
	merge_list( b ).K     = K;
	merge_list( b ).BIC   = BIC;

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

	merge_list( b ).MBIC  = MBIC - sum( [ cluster_list( [ u v ] ).BIC ] ) + BIC;

      end
      
      % ( 3.2 ) Pick the best choice

      [ tmp_MBIC, max_ind ] = max( [ merge_list.MBIC ] );

      % ( 3.3 ) Termination test

      if tmp_MBIC < MBIC
	% Stop if global BIC does not improve
	break
      end

      % ( 3.4 ) Apply the best merge

      u = merge_list( max_ind ).u;
      v = merge_list( max_ind ).v;

      cluster_list( u ).tag      = [ cluster_list( [ u v ] ).tag ];
      cluster_list( u ).mu       = merge_list( max_ind ).mu;
      cluster_list( u ).sigma    = merge_list( max_ind ).sigma;
      cluster_list( u ).n        = merge_list( max_ind ).n;
      cluster_list( u ).log_lkld = merge_list( max_ind ).log_lkld;
      cluster_list( u ).K        = merge_list( max_ind ).K;
      cluster_list( u ).BIC      = merge_list( max_ind ).BIC;

      cluster_list = cluster_list( [ 1:v-1 v+1:end ] );

      MBIC = merge_list( max_ind ).MBIC;

      % Sanity check
      if abs( MBIC - sum( [ cluster_list.BIC ] ) ) > 1e-8 * abs( BIC )
	error( [ mfilename ' is insane!' ] );
      end

      % ( 3.5 ) Termination test

      if numel( cluster_list ) < 2
	break;
      end
      
    end  % Iterative pairwise BIC merging

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


  % ( 4 ) Write out

  static_bic.classify_filename = classify_filename;
  static_bic.lambda            = lambda;

  save( output_filename, 'static_bic' );
  
  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
  

