function r = do_stc( r, T_short, N_future, stc_type )

% function r = do_stc( r, T_short, N_future, stc_type )
%
% Main implementation of short-term clustering,
% ( Maximum Likelihood or confident ).
% 
% "stc_type" can take two possible values: 'ml' or 'conf'.
%
% -> location estimates in the input "r" are partitioned. 
% -> in the output "r", two location estimates belong to
%    the same cluster iff they have the same tag r( i ).tag == r( j ).tag. 
%
% "do_stc" is typically called by either "do_ml" or "do_conf".
%
% This code can be useful to demonstrate short-term clustering.
% See IDIAP RR 04-14 (2004) or Lathoud's PhD thesis (2006).
%
% Notes:
% - matrices for binary graphs are only defined in the upper triangular part.
% - no particular code optimization was made.
%
% By Guillaume Lathoud, 2006 - lathoud@idiap.ch

  if nargin < 4
    error( [ mfilename ' needs 4 input arguments.' ] );
  end
  
  switch stc_type
   case 'ml'
    DO_CONF = 0;
   case 'conf'
    DO_CONF = 1;
   otherwise
    error( [ mfilename ': "stc_type" must be either ''ml'' or ''conf''.' ] );
  end

  % Confident clustering (DO_CONF==1) rely on a STRICT inequality test.
  % --> We have to be careful with the numerical imprecision of the machine.
  % --> the "10" factor is not extremely important (I tried 1, 10, 100, 1000
  %     and obtained almost the same result, correct in all cases). 
  CONF_EPSILON = 10 * eps;
  
  % Total number of location estimates
  N_r = numel( r );

  % Target: define "tags"
  % (locations estimates with the same tag belong to the same cluster)

  [ r.tag ] = deal( NaN );

  last_tag = 0;

  % We implement Table 6.1 of Lathoud's PhD thesis
  % + (if DO_CONF=1) confident clustering
  %  
  % For simplicity we optained for offline fitting of the bi-Gaussian models.
  % For an implementation with online fitting, check
  % http://mmm.idiap.ch/Lathoud/2006-multidetloc/

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % We can define all possible partitions of 
  % up to N_future elements in an offline manner 
  % -> let's do it!
  
  all_graphs = [];
  
  for a = 1:N_future
  
    if a < 2
      all_graphs( a ).graph_list = construct_graph_list_pruned( a, [], [], [], +Inf );
    else
      all_graphs( a ).graph_list = construct_graph_list_pruned( a, all_graphs( a-1 ).graph_list, [], [], +Inf );
    end
    
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( Step 1 ) Offline fitting of bi-Gaussian models
  % We use a somewhat simplified EM (see comments inside "fit_big_model()")
  %
  % -> in practice the precise values of sigma_same and sigma_diff
  % don't matter a lot. Only the orders of magnitude of these values are important.

  disp( [ mfilename ': Step 1: EM fitting of the bi-Gaussian models...' ] );
  
  big_list = [];
  
  a_par = [];
  a_par.init_threshold = pi / 8;  % Necessary because we use a somewhat simplified EM (see comments inside "fit_big_model()").
  a_par.max_sigma_same = +Inf;

  for T = T_short:-1:1
    
    % Get the data

    delta_theta = [];
    
    for n = 1:N_r
  
      ind = find( [ r.T ] - r( n ).T == T );

      if ~isempty( ind )
	delta_theta = [ delta_theta    r( ind ).theta - r( n ).theta ];
      end

    end
    
    % Make sure it is between -pi and +pi

    delta_theta = -pi + mod( pi + delta_theta, 2 * pi );

    % Fit a bi-Gaussian model

    big_list( T+1 ).T = T;
    big_list( T+1 ).big_model = fit_big_model( delta_theta, a_par );

    % Additional constraint: useful only on pathologic synthetic data,
    % not really necessary on real data

    a_par.max_sigma_same = min( a_par.max_sigma_same, big_list( T+1 ).big_model.sigma_same );

  end

  % T = 0:  probably not used on most synthetic data
  % We put something here, for completion

  big_list( 1 ).T = 0;
  big_list( 1 ).big_model = big_list( 2 ).big_model;

  % For convenience, stack standard deviations (std devs) in two arrays

  big_list_sigma_same = repmat( NaN, 1, T_short + 1 );
  big_list_sigma_diff = repmat( NaN, 1, T_short + 1 );

  for T = 0:T_short
    
    big_list_sigma_same( T + 1 ) = big_list( T + 1 ).big_model.sigma_same;
    big_list_sigma_diff( T + 1 ) = big_list( T + 1 ).big_model.sigma_diff;
    
  end
  
  disp( [ mfilename ': Steps 2+: ''' stc_type ''' short-term clustering of location estimates...' ] );

  for n = 1:N_future:N_r
    
    % Extract past and future
    
    T_n = r( n ).T;
    P_ind = find( ismember( [ r.T ], T_n-T_short:T_n-1 ) );
    P = r( P_ind );
    N_P = numel( P );

    F_ind = n:min( n+N_future-1, N_r );
    F = r( F_ind );
    N_F = numel( F );

    % Prepare likelihoods of pairwise decisions:
    %   H0: two location estimates belong to two DIFFERENT sources
    %   H1: two location estimates belong to the SAME source

    PUF = cat( 2, P(:).', F(:).' );
    PUF_ind = [ P_ind  F_ind ];
    N_PUF = N_P + N_F;

    T     = [ PUF.T ];
    theta = [ PUF.theta ];
    tag   = [ PUF.tag ];

    abs_delta_T = abs( repmat( T(:), 1, N_PUF ) - repmat( T, N_PUF, 1 ) );
    
    delta_theta = repmat( theta(:), 1, N_PUF ) - repmat( theta, N_PUF, 1 );
    
    loglkld_diff = repmat( NaN, N_PUF, N_PUF );
    loglkld_same = repmat( NaN, N_PUF, N_PUF );

    ind = find( abs_delta_T <= T_short );

    % H0: different sources
    loglkld_diff( ind ) = log_normpdf( delta_theta( ind ), 0, big_list_sigma_diff( 1 + abs_delta_T( ind ) ).' );

    % H1: the same source
    loglkld_same( ind ) = log_normpdf( delta_theta( ind ), 0, big_list_sigma_same( 1 + abs_delta_T( ind ) ).' );

    % The current partition
    % -> this is what we will modify
    % -> storage form: 0 or 1 on each edge, corresponding
    %    to pairwise decisions H0 or H1 (different sources or same source).
    %
    % Remember: ( NaN == NaN ) yields 0 (i.e. "not true").

    full_graph = sparse( N_PUF, N_PUF );
    
    ind_1 = 1:N_P;
    ind_2 = N_P + (1:N_F);

    full_graph( ind_1, ind_1 ) = triu( repmat( tag( ind_1 ), N_P, 1 ) == repmat( tag( ind_1 ).', 1, N_P ), 1 );
    
        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ( Step 2 ) Future alone
    
    % Define all possible partitions of F

    % --> already done (see "all_graphs" above)

    graph_list = all_graphs( N_F ).graph_list;

    % Choose the most likely partition of F

    if ~DO_CONF

      % ML clustering -> don't need the posteriors

      [ best_ind, F_graph_list, graph_like ] = select_graph( loglkld_same( ind_2, ind_2 ), loglkld_diff( ind_2, ind_2 ), graph_list );

    else

      % Confident clustering -> need the posteriors

      [ best_ind, F_graph_list, graph_like, post_same, post_diff ] = select_graph( loglkld_same( ind_2, ind_2 ) , loglkld_diff( ind_2, ind_2 ), graph_list );

      % Confident clustering: based on the posteriors,
      % implement the strict inequality to detect and isolate
      % location estimates around which "same/different" decisions are dubious.

      F_graph_list{ best_ind } = split_confident( F_graph_list{ best_ind }, post_same, post_diff, CONF_EPSILON );
      
    end
    
    full_graph( N_P + (1:N_F), N_P + (1:N_F) ) = F_graph_list{ best_ind };
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ( Step 3 ) Merge future with past

    if N_P > 0
      
      % Define all possible merges between P and F

      graph_1 = full_graph( ind_1, ind_1 );
      graph_2 = full_graph( ind_2, ind_2 );

      % Limitation of the permitted "same" edges
      % permitted_1 is a N x N matrix of binary values,
      % where "1" on an edge indicates a permitted "same" decision for this edge,
      % and   "0" on an edge indicates a forbidden "same" decision for this edge,
      permitted_1 = ones( N_P, N_F );
      permitted_1( find( isnan( loglkld_same( ind_1, ind_2 ) ) ) ) = NaN;

      graph_list = construct_merge_two( graph_1, graph_2, permitted_1, +Inf );

      % Choose the most likely merge between P and F
      
      if ~DO_CONF
	
	% ML clustering -> don't need the posteriors
	
	[ best_ind, graph_list ] = select_graph( loglkld_same, loglkld_diff, graph_list );

      else
	
	% Confident clustering -> need the posteriors
	
	[ best_ind, graph_list, graph_like, post_same, post_diff ] = select_graph( loglkld_same, loglkld_diff, graph_list );
	
	% Confident clustering: based on the posteriors,
	% implement the strict inequality to detect and isolate
	% location estimates around which "same/different" decisions are dubious.
	
	graph_list{ best_ind } = split_confident( graph_list{ best_ind }, post_same, post_diff, CONF_EPSILON );

	% Do not modify previous decisions taken at Step 2 ("future" windows)

	graph_list{ best_ind }( ind_1, ind_1 ) = graph_1;
	graph_list{ best_ind }( ind_2, ind_2 ) = graph_2;

      end
      
      % Store the merge
      
      full_graph( ind_1, ind_2 ) = graph_list{ best_ind }( ind_1, ind_2 );
      
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Return the result: convert the graph back to a list of tags
    
    for a = 1:N_PUF
      
      % Make sure this sample has a tag
      if isnan( tag( a ) )
	
	% Create a new cluster
	last_tag = last_tag + 1;
	tag( a ) = last_tag;
	
      end
      
      % Propagate the tag to other samples in the same cluster
      tag( find( full_graph( a, : ) ) ) = tag( a );
      
      % Are we done?
      if all( ~isnan( tag ) )
	break;
      end

    end
    
    % Store the result
    
    for a = 1:N_PUF
      r( PUF_ind( a ) ).tag = tag( a );
    end
    
  end  % of "for n = 1:N_future:N_r"

  
  % Sanity check

  if any( isnan( [ r.tag ] ) )
    error( [ mfilename ' is NaN-insane.' ] );
  end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              FUNCTIONS BELOW ARE IN ALPHABETICAL ORDER             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find all possible connected graphs of N points,
% where each possible edge belongs to the  "same" or "different" category.
% This is a recursive implementation: it calls itself with N-1.
% The difficulty is to make sure that each graph is consistent
% ( if A = B and B = C then A = C ).
%
% Output: GRAPH_LIST_N contains a list of graphs with N points,
% stored as a cell array of N - by - N matrices of binary values
% In each matrix, only the upper triangular part is used,
% where a 1 indicates "same" and a 0 indicates "different".
%
% An optional pruning is possible, through matrices "permitted_0" and
% "permitted_1". Each is a N - by - N matrix of ones (possible link)
% and zeros (forbidden link). Hence zeroes form the prior knowledge
% (default: ones(N,N)). Only the upper triangular part is used.
%
% N_GRAPHS_MAX is the maximum number of graphs to construct.
% If reached, we issue a warning and return the N_GRAPHS_MAX graphs.
% default value: +Inf.
%
% Each combination is stored as a sparse matrix ( only strictly upper triangular
% part is used, lower triangular part is left as zeroes ).

function graph_list_n = construct_graph_list_pruned( n, graph_list_n_minus_1, permitted_0, permitted_1, n_graphs_max )
  
  if nargin < 1
    error( sprintf( [ mfilename '.construct_graph_list_pruned[%d] needs at least one parameter' ] ,n) );
  end
  
  if ~isnumeric( n )
    error( sprintf( [ mfilename '.construct_graph_list_pruned[%d] needs an integer N >= 0' ] ,n) );
  end
  
  if rem( n, 1 ) | ( n < 0 )
    error( sprintf( [ mfilename '.construct_graph_list_pruned[%d] needs an integer N >= 0' ] ,n) );
  end
  
  n = floor( n );
  
  % No point: empty matrix
  if n < 1
    graph_list_n = { sparse( [], [], [], 0, 0 ) };
    return;
  end
  
  % List of graphs on which we build on
  % Default: empty
  if ~exist( 'graph_list_n_minus_1', 'var' )
    graph_list_n_minus_1 = [];
  end
  
  % Limitation of the permitted "different" edges
  % permitted_0 is a N x N matrix of binary values,
  % where "1" on an edge indicates a permitted "different" decision for this edge,
  % and   "0" on an edge indicates a forbidden "different" decision for this edge,
  if ~exist( 'permitted_0', 'var' )
    permitted_0 = [];
  end
  
  % Limitation of the permitted "same" edges
  % permitted_1 is a N x N matrix of binary values,
  % where "1" on an edge indicates a permitted "same" decision for this edge,
  % and   "0" on an edge indicates a forbidden "same" decision for this edge,
  if ~exist( 'permitted_1', 'var' )
    permitted_1 = [];
  end
  
  if ~exist( 'n_graphs_max', 'var' )
    n_graphs_max = [];
  end
    
  %%
  
  if isempty( permitted_0 )
    permitted_0 = ones( n );
  end
  
  if isempty( permitted_1 )
    permitted_1 = ones( n );
  end
  
  if isempty( n_graphs_max )
    n_graphs_max = +Inf;
  end
  
  % 1 or more points:
  % first make sure we have the result from the previous order
  if isempty( graph_list_n_minus_1 )
    graph_list_n_minus_1 = construct_graph_list_pruned( n - 1, [], permitted_0(1:n-1, 1:n-1), permitted_1(1:n-1, 1:n-1), n_graphs_max );
  end
  
  % Quick sanity check
  if size( graph_list_n_minus_1{ 1 }, 1 ) ~= n-1
    error( [ mfilename '.construct_graph_list_pruned() is insane' ] );
  end
  
  % Make sure we have not reached the maximum number of graphs yet
  if length( graph_list_n_minus_1 ) >= n_graphs_max
    disp( sprintf( 'construct_graph_list_pruned[%d]: warning! graph_list_n_minus_1 already has reached maximum number of graphs.', n ) );
    graph_list_n = graph_list_n_minus_1;
    return;
  end
  
  % Algorithm to construct all possible graphs of N points:
  %
  % For each possible graph of N-1 points:
  %
  %   Find the G PERMITTED cliques of connected points.
  %   PERMITTED cliques means that all its vertices of the clique have ones 
  %   in PERMITTED_0_POINTS and in PERMITTED_1_POINTS.
  %   There are G possible ways to add the nth point to an existing clique.
  
  permitted_0_points = permitted_0( 1:n-1, n );
  permitted_1_points = permitted_1( 1:n-1, n );
  
  graph_list_n = {};
  
  for a = 1:length( graph_list_n_minus_1 )
    
    prev_combi  = graph_list_n_minus_1{ a };
    [ link_i, link_j ] = find( prev_combi );
    
    point_visited = zeros( 1, n - 1 );
    
    cliques = {};
    
    
    % 1) Build cliques of connected points
    % for this combination
    current_clique = [];
    while ~all( point_visited )

      current_clique_size = length( current_clique );
      
      if isempty( current_clique )
	
	% Create new clique
	aux = find( ~point_visited );
	current_clique = [ aux( 1 ) ];
	point_visited( aux( 1 ) ) = 1;
	
      else
	
	% Expand current clique
	v = [];
	v = [ v link_j( find( ismember( link_i, current_clique ) ) ).' ];
	v = [ v link_i( find( ismember( link_j, current_clique ) ) ).' ];
	if ~isempty( v )
	  current_clique = unique( [ current_clique, v ] );
	  point_visited( v ) = 1;
	end
	
      end
      
      % Current clique finished?
      if ( length( current_clique ) == current_clique_size ) | all( point_visited )
	% Yes, store it and start a new one
	cliques{ end+1 } = current_clique;
	current_clique = [];
      end
      
    end
    
    % 2) Which cliques are permitted?
    
    permitted_cliques = [];
    for b = 1:length( cliques )

      if all( permitted_1( cliques{ b } ) ) & all( permitted_0( ~cliques{ b } ) )
	permitted_cliques = [ permitted_cliques  b ];
      end
    
    end
    
    % 3) Based on these permitted cliques build new combinations
    
    G = length( permitted_cliques );

    % Set the right size
    prev_combi( n, n ) = 0;
    
    % Don't forget the case where the point is isolated
    graph_list_n{ end + 1 } = prev_combi;
    
    for b = 1:G
      
      the_clique = cliques{ permitted_cliques( b ) };
      p = length( the_clique );
      graph_list_n{ end + 1 } = prev_combi + sparse( the_clique, n * ones( 1, p ), ones( 1, p ), n, n );
    
      if length( graph_list_n ) > n_graphs_max
	disp( sprintf( ['construct_graph_list_pruned[%d]: warning! reached maximum number of graphs.'], n) );
	graph_list_n = { graph_list_n{ 1:n_graphs_max } };
	return;
      end
      
    end
    
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function graph_list = construct_merge_two( in_graph1, in_graph2, permitted_1, n_merges_max, n_elements_min )

% function graph_list = construct_merge_two( in_graph1, in_graph2, permitted_1, n_merges_max, n_elements_min )
%
% Construct all possible merges of two graphs.
%
% ( used to select the best merge of the past window and the future window ).
%
% in_graph1 and in_graph2 are N x N matrices of binary values.
%
% Optional limitation of the permitted "same" edges:
% permitted_1 is a N x N matrix of binary values,
% where "1" on an edge indicates a permitted NEW "same" decision for this edge,
% and   "0" on an edge indicates a forbidden NEW "same" decision for this edge.
%
% n_merges_max limits the number of constructed graphs (default: +Inf).

  
  % Deal with parameters
  
  if nargin < 3
    error( [ mfilename 'construct_merge_two() needs at least 3 input arguments' ] )
  end
  
  if ~exist( 'n_merges_max', 'var' )
    n_merges_max = +Inf;
  end
 
  if ~exist( 'n_elements_min', 'var' )
    n_elements_min = 1;
  end
  
  % Quick check
  siz1 = size( in_graph1 );
  siz2 = size( in_graph2 );
  siz3 = size( permitted_1 );
  
  if siz1(1)~=siz1(2)
    error( [ mfilename 'construct_merge_two(): "in_graph1" must be a square matrix !' ] );
  end
  n1 = siz1( 1 );
  
  if siz2(1)~=siz2(2)
    error( [ mfilename 'construct_merge_two(): "in_graph2" must be a square matrix !' ] );
  end
  n2 = siz2( 1 );
  
  if any( [ n1 n2 ] ~= siz3 )
    error( [ mfilename 'construct_merge_two(): wrong size for "permitted_1" !' ] );
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % If nothing to do, just return
  
  if ~any( permitted_1(:) == 1 )
    
    new_graph = in_graph1;
    n = 1:n2;
    new_graph( end+n, end+n ) = in_graph2;
    
    graph_list{ 1 } = new_graph;
    
    return;
  
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Otherwise get to work
  % First, convert graphs to tags
  
  tag1 = graph_to_tag( in_graph1 );
  tag2 = max( tag1 ) + graph_to_tag( in_graph2 );
  
  tag1_list = unique( tag1 );
  tag2_list = unique( tag2 );

  ntag1 = length( tag1_list );
  ntag2 = length( tag2_list );

  % List the clusters we refuse to merge with any other cluster
  % Elements of "exclusion1_list" are indices in "tag1_list"
  % Elements of "exclusion2_list" are indices in "tag2_list"
  exclusion1_list = [];
  exclusion2_list = [];
  
  if n_elements_min > 1
    
    for a = 1:length( tag1_list )
      a_tag = tag1_list( a );
      n_elements = sum( tag1 == a_tag );

      if n_elements < n_elements_min
	exclusion1_list( end+1 ) = a;
      end
    
    end
    
    for a = 1:length( tag2_list )
      a_tag = tag2_list( a );
      n_elements = sum( tag2 == a_tag );

      if n_elements < n_elements_min
	exclusion2_list( end+1 ) = a;
      end
    
    end
    
  end
  
  % Convert permitted_1 into a list of excluded merges
  %
  % forbidden_merge is a cell array of ntag2 vectors of indices,
  % each index is an index of an element the tag1 vector.
  
  forbidden_merge = cell( 1, ntag2 );
  
  m = zeros( ntag1, ntag2 );
  
  [u,v] = find( permitted_1 == 0 );
  
  if ~isempty( u );
    
    % Convert to tag 
    u = tag1( u );
    v = tag2( v );
    
    % Convert to index in each tag list
    tag1_to_index1( tag1_list ) = 1:ntag1;
    tag2_to_index2( tag2_list ) = 1:ntag2;
    
    u = tag1_to_index1( u );
    v = tag2_to_index2( v );
    
    % Mark the forbidden merges
    m( sub2ind( size( m ), u, v ) ) = 1;
    
  end
  
  for a = 1:ntag2
   
    % For each cluster in graph2
    % we list the forbidden clusters in graph1
    
    forbidden_merge{ a } = find( m( :, a ) );
    
  end
    
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Second, define all possible merges between clusters
  % of in_graph1 and clusters of in_graph2. 
  % ( No merge between clusters of the same graph. )
  % 
  % Recursive implementation:
  % Each element of "merge_list" is a ntag1 x ntag2 sparse matrix of zeroes and ones,
  % such that there is only one "one" on each row.
  % We recursively build a list of (n+1)-column matrices from each n-column matrix.
  % 
  % I have no time for a better explanation
  % If this is not clear: lathoud@idiap.ch
  
  merge_list = { sparse( [], [], [], ntag1, 0 ) };

  reached_max = 0;
  
  while ~reached_max  &  size( merge_list{ 1 }, 2 ) < ntag2     % Do we have the right number of columns?
    
    % If not, then add a column : store all possibilities
    % that respect the constraint of only one "one" on each row
    
    new_merge_list = {};
    
    for a = 1:length( merge_list )
      
      if length( new_merge_list ) >= n_merges_max
	s = warning;
	if ~strcmp( s( 1 ).state, 'off' )
	 disp( 'construct_merge_two: WARNING! Reached the maximum number of merges.' ); 
	end
	
	% Fall back on the existing list of merges
	reached_max = 1;
	break;
	
      end
      
      m = merge_list{ a };
      
      % We can always add a column of zeroes (means no additional merge)
      
      new_merge_list{ end+1 } = [ m  zeros( ntag1, 1 ) ];
      
      if ~ismember( size( m, 2 ) + 1, exclusion2_list )
      
	% All possible non-zero columns
	
	if size( m, 2 ) < 1
	  ind = 1:size( m, 1 );
	else
	  ind = find( ~any( m, 2 ) );
	end
	
	ind = setdiff( ind, [ exclusion1_list  forbidden_merge{ size( m, 2 ) + 1 } ] );
	
	for b = ind(:).'
	  new_merge_list{ end+1 } = [ m  (1:ntag1).'== b ];
	end
	
      end
      
    end

    if reached_max
      break;
    end
    
    % Update with the new list of matrices with n+1 columns
    
    merge_list = new_merge_list;
    
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Third convert back to graphs
  % Each vector of tags becomes a matrix of binary values.

  ind = [];
  
  graph_list = {};
  
  for a = 1:length( merge_list )
    
    [u,v] = find( merge_list{ a } );
    v = tag2_list( v );
    
    m = ( repmat( tag2(:), 1, length( v ) )  ==  repmat( v(:).', length( tag2 ), 1 ) );
    [ r, c ] = find( m );
    
    tag2bis = tag2;
    tag2bis( r ) = u( c );
    
    tag = [ tag1 tag2bis ];
    
    graph_list{ a } = tag_to_graph( tag );
    
    if isempty( ind )
      m = zeros( size( graph_list{ a } ) );
      m( 1:n1, n1+1:end ) = permitted_1;
      ind = find( isnan( m ) );
    end
    
    graph_list{ a }( ind ) = 0;
    
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Fourth make sure the constructed graphs
  % don't have the "extra ones"
  % ( e.g. for T > T_short )
  
  if ~all( permitted_1 == 1 )
    
    for a = 1:length( graph_list )
      graph_list{ a }(1:n1, n1+1:end ) = graph_list{ a }(1:n1, n1+1:end ) & ( permitted_1 == 1 );
    end
    
  end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function big_model = fit_big_model( delta_theta, in_p )

function big_model = fit_big_model( delta_theta, in_p )
  
  if nargin < 2
    error( [ mfilename '.fit_big_model() needs two input arggments!' ] );
  end
  
  p = in_p;
  
  p_default.std_floor       = 1e-10;   % To avoid "warning log zero" messages
  p_default.verbose         = 0;
  p_default.init_threshold  = +Inf;
  p_default.max_sigma_same  = +Inf;
  p = fill_default( p, p_default );
  
  
  if p.verbose
    disp( [ mfilename '.fit_big_model() parameters:' ] )
    disp( p );
  end
  
  
  last_sigma_same = [];
  last_max_abs_delta_theta = [];
  
  delta_theta = [ -pi +pi delta_theta ];

  % Initialization
  if isempty( last_max_abs_delta_theta )
      label = abs( delta_theta ) < p.init_threshold;
  else
    label = abs( delta_theta ) < last_max_abs_delta_theta;
  end

  % Expectation-Maximization (EM) fit of a bi-Gaussian model
  % Note: this is a very simplified version, as posteriors are forced 
  % to be either 0 or 1, through the "label" array.
  %
  % For a complete implementation of EM,
  % including automatic initialization with multiple thresholds,
  % please look at the "update_stc_model" function of:
  % http://mmm.idiap.ch/Lathoud/2006-multidetloc/FAST_detect_locate.m
  
  finished = 0;
  old_label = zeros( size( delta_theta ) );  % for convergence test
  while ~finished
    
    % Convergence?
    finished = 0.01 > ( length( find( old_label ~= label ) ) / length( label ) );
    
    aux = find( label );
    if isempty( aux )
      
      % No "connection" model available for this delay
      % This happens on some very specific synthetic data
      % This does not happen on real data
      
      mu_same    = [];
      sigma_same = [];
      mu_diff    = [];
      sigma_diff = [];
      finished   = 1;
      
    else
      
      % Retrain models
      ind = find( label );
      mu_same         = 0;
      sigma_same      = p.std_floor;
      if length( ind ) > 1
	sigma_same      = sqrt( sum( ( delta_theta( ind ) - mu_same ) .^ 2 ) / ( length( ind ) - 1 ) );
      end
      
      ind = find( ~label );
      mu_diff = 0;
      sigma_diff = sqrt( sum( ( delta_theta( ind ) - mu_diff ) .^ 2 ) / ( length( ind ) - 1 ) );
      
      % Enforce causality
      % This has never been observed to be needed on real data
      % However it is not costly to keep it!
      
      if ~isempty( p.max_sigma_same )
	sigma_same = min( sigma_same, p.max_sigma_same );
      end
      
      % Floor the std dev to avoid "warning log zero" messages
      % This is necessary to cope with VERY specific synthetic data examples
      % This is not necessary on real data
      
      sigma_same      = max( sigma_same, p.std_floor );
      sigma_diff   = max( sigma_diff, p.std_floor );
	
      % Relabel samples
      old_label = label;
      label = normpdf( delta_theta, mu_same, sigma_same ) > normpdf( delta_theta, mu_diff, sigma_diff );
      
      if p.verbose
	x = -179.5:179.5;
	y = delta_theta / pi * 180;
	n = hist( y, x );
	n2 = hist( y( find(label)),x);
	figure(1); plot(x,[n;n2].');
	title( sprintf( [ mfilename '.fit_big_model()' ] ) );
	pause
	
      end
      
    end
    
  end
  
  % Store the model, if any
  big_model.mu_same       = mu_same;
  big_model.sigma_same    = sigma_same;
  big_model.mu_diff    = mu_diff;
  big_model.sigma_diff = sigma_diff;
  
  % Store 2 values needed to enforce causality
  % (sigma is the strict necessary)
  if ~isempty( mu_same ) & ~isempty( sigma_same )
    last_sigma_same = sigma_same;
    last_max_abs_delta_theta = max( [ p.std_floor abs( delta_theta( find( label ) ) ) ] );
  end
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert a N x N graph (matrix of binary values)
% into a vector of N tag values.

function tag = graph_to_tag( graph )

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

  N = size( graph, 1 );
  n_tags = 0;

  tag = zeros( 1, N );

  for a = 1:N

    % Make sure this vertice has a tag
    if ~tag( a )
      n_tags = n_tags + 1;
      tag( a ) = n_tags;
    end

    % Propagate this tag to connected vertices
    tag( find( graph( a, : ) ) ) = tag( a );

  end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Select the ML graph out of list,
% given "same" likelihoods and "different" likelihoods.

function [ best_ind, graph_list, graph_like, post_same, post_diff ] = select_graph( loglkld_same, loglkld_diff, graph_list )

  n = size( loglkld_same, 1 );
  ng = length( graph_list );
  
  % ( 1 ) Compute graph likelihoods
  
  graph_like = repmat( NaN, 1, ng );
  
  % Likelihood computation is restricted to |t_i - t_j| <= T_{short}
  bin_mask = ~isnan( loglkld_same );
  
  for a = 1:ng

    g = graph_list{ a };
    graph_like( a ) = sum( loglkld_same( find( g & bin_mask ) ) ) + sum( loglkld_diff( find( ~g & bin_mask ) ) );
    
  end

  % Pick the ML one
  [tmp, best_ind] = max( graph_like );

  % Do we need posteriors?

  if nargout < 4

    % No
    return;

  end

  % Yes

  % ( 2 ) Compute posteriors for each possible local decision,
  % using all possible graphs.
  %
  % This can be useful for confident clustering
  % ...but definitely useless in audio meeting recordings (e.g. M4 corpus).
  
  graph_scaled = graph_like - my_logsum_fast( graph_like );
  
  post_diff = repmat( NaN, n, n );
  post_diff( find( bin_mask ) ) = -Inf;
  
  post_same = repmat( NaN, n, n );
  post_same( find( bin_mask ) ) = -Inf;
  
  for a = 1:ng
  
    g = graph_list{ a };
    
    ind = find( ~g & bin_mask );
    ind = ind(:).';
    m = post_diff( ind );
    m( 2,: ) = graph_scaled( a );
    post_diff( ind ) = my_logsum_fast( m );
    
    ind = find( g & bin_mask );
    ind = ind(:).';
    m = post_same( ind );
    m( 2,: ) = graph_scaled( a );
    post_same( ind ) = my_logsum_fast( m );
    
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function graph = split_confident( graph, post_same, post_diff, conf_epsilon )

% function graph = split_confident( graph, post_same, post_diff, conf_epsilon )
%
% Confident clustering: based on the posteriors,
% implement the strict inequality to detect and isolate
% the location estimates around which 
% the "same/different" decisions are dubious.
  
  if nargin < 4
    error( [ mfilename '.split_confident() needs 4 input arguments.' ] );
  end
  
  npoints = size( post_same, 1 );

  % ( 1 ) Determine which points we need to isolate,
  % based on a comparison between "post_same" and "post_diff"
  
  points_to_isolate = [];

  bin_mask = ~isnan( post_same );
  
  % List the "same source" pairwise decisions, in the graph
  bin_same =  graph & bin_mask;

  % List the "different sources" pairwise decisions, in the graph
  bin_diff = ~graph & bin_mask;

  % Now detect the points that we need to isolate
  
  if any( bin_same(:) ) & any( bin_diff(:) )
    
    max_post_same = max( post_same( find( bin_same ) ) );
    
    % Implement the STRICT inequality test
    
    [ u, v ] = find( bin_diff & ( post_diff < max_post_same - conf_epsilon ) );

    points_to_isolate = unique( [ u; v ] );

  end

  % ( 2 ) Isolate those points (which also breaks clusters)

  if ~isempty( points_to_isolate )

    tag = graph_to_tag( graph );
    
    for zz = points_to_isolate(:).'
      
      % Create a new 1-element cluster
      % with in it, ONLY the point to be isolated.

      old_tag   = tag( zz );
      tag( zz ) = max( tag ) + 1;
      
      % Break the rest of the cluster into two parts
      % ( before and after the point to isolate )
      % -> the part "after" gets a new tag.

      [ tmp, ind ] = find( tag( zz+1:end ) == old_tag );
      
      if ~isempty( ind )
      
	tag( zz + ind ) = max( tag ) + 1;
	
      end

    end

    graph = tag_to_graph( tag );

  end
  

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert a vector of N tag values
% into a N x N graph (matrix of binary values)

function graph = tag_to_graph( tag )

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

  n = length( tag );
  graph = sparse( triu( repmat( tag(:), 1, n ) == repmat( tag(:).', n, 1 ), 1 ) );
